Segmentation of the Izu-Bonin and Mariana slabs based on the analysis of the Benioff seismicity distribution and regional tomography results

We present a new model of P and S velocity anomalies in the mantle down to a depth of 1300 km beneath the Izu-Bonin and Mariana (IBM) arcs. This model is derived based on tomographic inversion of global travel time data from the revised ISC catalogue. The results of inversion are thoroughly verified using a series of different tests. The obtained model is generally consistent with previous studies by different authors. We also present the distribution of relocated deep events projected to the vertical surface along the IBM arc system. Unexpectedly, the seismicity forms elongated vertical clusters instead of horizontal zones indicating phase transitions in the slab. We propose that these vertical seismicity zones mark zones of intense deformation and boundaries between semi-autonomous segments of the subducting plate. The P and S seismic tomography models consistently display the slab as prominent high-velocity anomalies coinciding with the distribution of deep seismicity. We can distinguish at least four segments which subduct differently. The northernmost segment of the Izu-Bonin arc has the gentlest angle of dipping which is explained by backward displacement of the trench. In the second segment, the trench stayed at the same location, and we observe the accumulation of the slab material in the transition zone and its further descending to the lower mantle. In the third segment, the trench is moving forward causing the steepening of the slab. Finally, for the Mariana segment, despite the backward displacement of the arc, the subducting slab is nearly vertical. Between the Izu-Bonin and Mariana arcs we clearly observe a gap which can be traced down to about 400 km in depth. Based on joint consideration of the tomography results and the seismicity distribution, we propose two different scenarios of the subduction evolution in the IBM zone during the recent time, depending on the reference frame of plate displacements. In the first case, we consider the movements in respect to the Philippine Plate, and explain the different styles of the subduction by the relative backward and forward migrations of the trench. In the second case, all the elements of the subduction system move westward in respect to the stable Asia. Different subduction styles are explained by the “anchoring” of selected segments of the slab, different physical properties of the subducting plate and the existence of buoyant rigid blocks related to sea mount and igneous provinces.


Introduction
Subduction zones are the most spectacular geological objects where intensive geodynamic processes occur. Strong mechanical deformations and complex geochemical processes in subduction zones are responsible for the diversity of various geological structures observed on the surface. Investigating the links between surface geology and processes in the subducting slab is an actual problem that involves many specialists from different fields of geology and geophysics (e.g., Dobretsov et al., 2012). The general concept of subduction is confirmed by a number of different data and accepted by most specialists in the world. However, the details of this process are the subject of intensive debates. In particular, representation of the subduction zones as a two-dimensional conveyer-type process appears to be in most cases oversimplified and it does not explain many important phenomena. For example, the distribution of deep seismicity, arc  Peive, (1980): Pre-Cretaceous (green), Pre-Mesozoic (red), Pre-Miocene (yellow), Pre-Quaternary age (white), sea floor of Pacific Ocean, Pre-Triassic (blue). White dotted lines depict arcs (ancient and present) (Peive, 1980). Names and ages of arcs are indicated. Arrows with numbers present the displacement rates in respect to Philippine plate (DeMets et al., 2010). The red triangles indicate active volcanoes (including submarine and hydrothermal volcanoes) (Taylor and Nesbitt, 1998;Baker et al., 2008). volcanism properties and dip angles of slabs may abruptly change along the strike of the subduction zone. There are many open questions related to subduction zones and, to answer them, new robust and reliable data on the deep structure are still required.

Geological overview
Izu-Bonin-Mariana (IBM) arc system is located in the western Pacific (Fig. 1). The arc length is more than 2800 km and it extends from north to south, from Tokyo Bay to Guam. This intraoceanic convergent zone is the result of a multistage subduction of the Pacific plate beneath the Philippine plate. There are the remnants of ancient arcs, such as the Kyushu-Palau Ridge with the age of volcanism 37-25 million years (Peive, 1980;Berckhemer et al., 1982), the West Mariana Ridge with the age of volcanism 23-10 million years and Bonin Ridge. There are also inactive remnants of rift systems, such as the Shikoku Basin, Parece Vela Basin (extension age 25-15 Ma) and the West Philippine Basin (with the age of oceanic plate 100-150 Ma). There is also a presently active back-arc spreading basin of the Mariana Trough, where the extension started about 3 million years ago (Peive, 1980, Fig. 1) and occur at a rate of 15-30 mm yr −1 (Kato et al., 2003). The age of the Pacific Plate along the IBM arc system is generally old and it varies from 146 Ma in the northern segments to 155 Ma at the Mariana trench (e.g., Müller et al., 1997).
There are more than 20 volcanic islands along the Izu-Bonin arc and more than 70 along the Mariana arc (Taylor and Nesbitt, 1998;Baker et al., 2008). Most of the volcanoes along the IBM arc system are submarine. Figure 1 shows the location of the volcanoes and tectonic units of different ages. The volcanoes of the IBM system are generally characterised by low dacite and andesite content (Kelemen et al., 2004) and a frequent occurrence of the boninite series, which were first discovered and described in this region (Petersen, 1891;Crawford et al., 1981). These composition features make the IBM volcanoes particular compared to other subduction zones of the world.

Previous slab structure investigations
The structure of the subducting slabs beneath the IBM system have been investigated in a number of different studies. The first information on the shape of the subduction zone was derived from the analysis of depth distribution of seismicity in the Benioff zone (e.g., Katsumata and Sykes 1969;Engdahl et al., 1998). More details on the slab structure beneath the IBM area were found in a number of seismic models on the global (e.g., Bijwaard et al., 1998) and regional scales (Widiyantoro, 1999;1996, Van der Hilst and Seno;, Gorbatov et al., 2003Miller et al., 2005Miller et al., , 2006aMiller et al., , b, 2011 mostly constructed by a tomographic inversion of body wave travel times. These studies reveal generally consistent information on the shape of the subducting slab, but in details they appear to be considerably different. For example, the penetration of the slab into the lower mantle is still actively debated, and the behaviour of the slab in the contact zone between Izu-Bonin to the Mariana segments is not clear yet. The inconsistencies between different models are partly due to the fact that seismic tomography based on noisy earthquake data is not a perfect tool; often processing of same data by different authors leads to considerably different results. Therefore, any attempt to compute a new model using independent approaches is an important step to corroborate the existing information.

Kinematics of IBM area
The main purpose of this paper is to compute a new 3-D seismic model of the mantle beneath the IBM system and to use it to reconstruct the plate interactions and the subduction development in this area. This reconstruction needs robust information about the plate motions in the region of interest including the surrounding areas in a large extent. We found that the data on the rates of plate motions in and around the Philippine plate given by different authors are often not consistent and even contradictory. One of the most important problems is defining an appropriate reference frame which is necessary for accessing the subduction motion with respect to the lower mantle. For example, when considering the Philippine plate as a reference point, De Mets et al. (2010) show that the extension along the backarc Mariana Trough causes eastward displacement of the Mariana trench at about 2 cm yr −1 that forms a loop shaped arc (Fig. 1). In this reference frame, the rate of the Pacific plate westward motion near the trench with respect to the Philippine Plate is 4.9 and 4.0 cm yr −1 in the northern and southern segments of the Izu-Bonin arcs; in the Mariana arc it changes from 2.8 to 1.6 cm yr −1 when moving southward (DeMets et al., 1994(DeMets et al., , 2010. On the other hand, there are several frames based on different reference points which provide considerably different absolute rates of the plate motions in the IBM region. For example, Schellart et al. (2008) presents an overview of plate motions related to three reference frames which are related to different geological structures. Based on all these frames, they report a westward absolute displacement for the Mariana subduction of both subducting and overriding plates at the rates of 11.8-5 cm yr −1 , 8.1 -1.8 cm yr −1 and 7.1 -1.5 cm yr −1 , respectively.
The behaviour of the slab beneath the IBM area has been actively discussed in a series of papers by Miller et al. (2005Miller et al. ( , 2006aMiller et al. ( , 2006b. They reconstructed the position of the trench in the last 30 Ma and explained the variable shapes of the slab along the IBM arc system by the mutual displacements of the trench and the oceanic plate. The bend around the contact zone between Izu-Bonin and Mariana arc is explained by presence of highly buoyant segment in the incoming plate with large igneous Ogasawara Plateau. They explain the vertical shape of the Mariana slab by forward displacement of the trench (figure 10 in Miller et al., 2006b). There are many studies showing the evidence of eastward displacement of the Mariana segment of the trench (Seno and Maruyama, 1984;Seno et al., 1993;Hall et al., 1995aHall et al., , 1995bHall, 2002). Alternatively, the westward displacement of the trench is shown by (Carlson and Mortera-Gutierrez, 1990;Schellart et al., 2008;O'Neill et al., 2005;Kreemer et al., 2003;Gripp and Gordon, 2002). It was shown, by several authors, that the trench had several stages of migration at different phases of developing this subduction zone. For example, Sdrolias and Muller (2006) show three main stages of Mariana trench with a stable trench at 50-25 Ma, a backward displacement of the trench at 25-5 Ma and a forward displacement of the trench at 5-0 Ma. Faccena et al. (2009) show that the trench migration depends on the pull force of the slab, which increases when the lithosphere is older. The trench advancing occurs in the case of subduction of older and stronger oceanic lithosphere. Taking into account this effect (also considering other parameters including width and rheology of the lithosphere) the authors suggest that recent advancing of the IBM trench system occurred after a long-term history of the trench retreat. In turn, the trench displacement (advancing and retreating) controls the subduction style and the slab shape. It was shown in several papers (Van der Hilst and Heuret et al., 2007;Faccena et al., 2009) that in the case of the trench retreating, the subduction becomes shallower, whereas the trench advancing causes steeping of the slab of even backward turning if the bottom of the slab is anchored in the transition zone.
In this study, we present a new seismic model of the mantle seismic structure down to 1200 km depth beneath the IBM arc system. For this model, we use a more recent version of the ISC catalogue compared to the previous studies which contains considerably larger data amount. We use our own algorithm which has some particularities compared to the previous studies. We present several tests which demonstrate the limitations of the resolution and support the robustness of the derived structures. In this paper, we compare our model with the existing ones and discuss possible geodynamical interpretation of the results.

Data processing and algorithm
In this study, we use the travel times of P and S waves from the International Seismological Centre (ISC) for the time period from 1964 to 2007. The advantage of the ISC dataset is a long registration period and the global data coverage which cannot be achieved by any regional network. At the same time, the data quality in the initial ISC catalogue is low. The reported coordinates of events are strongly biased due to outdated location algorithm and a big number of outliers in the data. Therefore, before using for tomographic inversion, this catalogue should be reprocessed. Most of scientists use a revised EHB catalogue which was created from careful refining of the ISC data by Engdahl et al. (1998). For our studies, we perform our own revision in which the rejection of outliers and relocation of sources were performed using the algorithm described in (Koulakov and Sobolev, 2006). This revised catalogue has been constructed based on the 1D reference model AK 135 (Kennett at al., 1995) and takes into account the ellipticity of the Earth and variable crustal thickness according to the global model CRUST2.0 (Bassin et al., 2000). Our approach is less conservative compared to the one used by Engdahl et al. (1998), and it keeps more data, which is favourable for tomographic inversion.

K. Jaxybulatov et al.: Segmentation of the Izu-Bonin and Mariana slabs
Regional tomographic inversion in this study is based on the approach proposed by (Koulakov et al., 2002). In later works, this algorithm was modified and used to study the upper mantle structure beneath some regions (see, for example, Koulakov et al., 2009a;Koulakov, 2011). This tomographic inversion algorithm is based on the linearised approach which implies that all calculations are performed on the basis of rays constructed in one-dimensional model AK135 (Kennett at al., 1995). It should be noted that there is no theoretical difficulty in using a nonlinear algorithm with a three-dimensional ray tracing (for example, as used in the LOTOS code by Koulakov, 2009b). However, we do not believe that implying the iterative nonlinear inversions would bring principal changes in the results. First, the recovered amplitudes of anomalies in the upper mantle are relatively low (2-3 %) and they do not really lead to a considerable nonlinear effect. Second, the amplitudes derived from the inversion of the ISC data are usually lower than expected in nature because of the large damping which should be applied to reduce the noise-related instability. In this case, even if we iteratively trace the rays in 3-D, they might appear not significantly closer to the true paths than the 1D based rays. Third, the linear approach is fast, and it gives us a possibility for performing many trials based on real and synthetic data. This allows us to perform thorough selection of free inversion parameters which would be hard in the case of a single nonlinear inversion.
In this study, the tomographic inversion is performed separately in three overlapping circular windows with the radius of 10 • , covering the study area. For the IBM area we use the travel times from the revised ISC catalogue which mostly correspond to events located in the study region which were recorded by worldwide stations at all possible epicentral distances. The opposite scheme, with stations in the target area recording the teleseismic events, is almost not used in this case, because of the lack of seismic networks in the study region. The total amount of used rays for the entire region is more than 2 million. The amount of P data is about ten times larger than that of the S-data. Note that using this ray configuration does not allow obtaining high resolution in tomography reconstruction for the shallowmost sections (above ∼50 km depth) where the ray coverage is poor and there is a strong trade-off between velocity and source parameters.
For the inversion, the 3-D velocity distribution is parameterised with a set of nodes which are distributed in 15 depth levels (50, 100, 150, 220, 290, 360, 430, 500, 570, 650, 710, 800, 900, 1000, 1300 km). In each level, the nodes are installed according to the ray coverage; no nodes are set in areas with insufficient data amount. To minimise any artifacts related to grid configuration, the velocity models are calculated independently for four differently oriented grids with basic orientations at 0 • , 22 • , 45 • and 67 • and then averaged. The inversion was performed simultaneously for Pand S-velocity anomalies and for the source parameters (co-ordinates and origin times). The damping coefficients and weighting of different parameters were estimated based on results of synthetic modelling. The inversion of the matrix is performed using the LSQR algorithm (Page and Saunders, 1982; Van der Sluis and Van der Vorst, 1987). The main model was obtained by combining 12 separate inversion results (4 models in each of the three circular areas).

Slab-related seismicity
We show the distribution of seismicity from the revised ISC catalogue projected to the vertical profile along the IBM system in Fig. 2. The map view of the profile location and the band, from which the events were considered, is shown in Fig. 2a. The projection of the events to the profile was performed in each zone along the lines with proportional angles in respect to the orientations of opposite sides of the trapezium. We present both the individual events with magnitude indications (Fig. 2b) and the integral released energy in cells 50 × 50 km (Fig. 2c).
Plots in Fig. 2 display some important features which can be used for geodynamic interpretation. The paradox conclusion which can be drawn from this representation of the seismicity distribution is that the events located below 200 km depth beneath the IBM arc system form rather vertically oriented clusters than horizontal ones. This seems to contradict a point of view that seismicity in the slab reflects phase transitions which occur at certain depth. So we would rather expect horizontal clusters around 100 km, 200 km and other depth levels where the material of the slab is transformed. They can be seen locally, but less prominently than thin vertical zones of seismicity. The most clear are vertical clusters of seismicity around the distances of 700, 1380 and 1910 km. The later cluster corresponds to the Mariana segment with the seismicity along a relatively narrow vertical "column" of about 700 km high and only 100-200 km width.
Another striking feature is a large seismicity cluster related to the lower part of the northern segment of the Izu-Bonin plate. Between 0 and 900 km along the profile, it tends to dip from 200 km to 300 km depth for the upper boundary of the cluster, and from 400 km to 500 km for the lower boundary. The southern limit of this body roughly fit to the vertical cluster at ∼700 km distance.
These observations possibly indicate that along the profile there are at least four separate segments of the subducting plate which may behave differently. The first northern segment of the Izu-Bonin plate is represented by gently dipping plate with very active seismicity in the depth interval between 300 and 500 km. The next two segments of the Izu-Bonin arc are limited by vertical clusters at 700 km, 1379 km and 1912 km. Unlike the first segment, for these two segments, no significant deep seismicity is observed. The cluster in the middle at 1379 km is observed down to about 350 km in depth, which may imply that these two segments are only

Computed P-and S-velocity models
The derived models of seismic P-and S-anomalies are shown in seven horizontal sections (Fig. 3) and in 9 vertical sections for Izu-Bonin-Mariana arc (Fig. 4). The areas, where the distance to the nearest grid node is more than 80 km, are shaded; in these areas the data coverage is not sufficient for the tomographic inversion. The variance reduction after the inversion was more than 50 % which shows relatively high level of coherent signal.
The major feature which is clearly seen in horizontal and vertical sections is a slab-related high-velocity anomaly observed beneath the entire arc, except for the shallowmost part in the Mariana region where the resolution is poor. This anomaly is generally consistent in P and S model, although the amount of the S data was about ten times smaller than that of P. The correlation between the P and S models is an informal argument for the reliability of the model; on a regional scale, the major geological structures, such as subduction, usually affect similarly the P and S velocities. Note that this anomaly fits to the deep seismicity distribution in the Benioff zone (see Fig. 4).
By observing the vertical and horizontal sections, it helps to determine the shape of the subducting slab. In vertical Sect. 1, we observe the southern part of the Japan arc. The amount and the quality of data in this area are maximal, therefore, the image seems to be well resolved. Section 2 roughly crosses the transition zone of Japan arc to Izu-Bonin. In these sections, the shape of the slab is well repeated in Pand S-anomalies, excluding some features in deepest parts. The recovered slab-related anomaly fits to the distribution of deep seismicity. In these sections, the slab gently dips with the angle of about 35 • . For Sect. 3, the slab becomes steeper, but in the transition zone between 400 and 600 km depth it turns to be horizontal. In this section, an enigmatic feature is a prominent positive anomaly in the lower mantle which is clearly seen in both P and S models, though with different shapes. Note that in the neighbouring Sects. 2 and 4, this anomaly is not observed. This anomaly is not distinctly detected and discussed in other early tomography models, but in certain papers (e.g., Widiyantoro et al., 1999;Miller et al., 2006) anomalies penetrating into the lower mantle and sometimes separate anomalies lower than 660 km are observed. So presence of this feature can be taken into account.
In Sects. 4 and 5, the dipping part of the slab becomes steeper and reaches the angle of 60 • -70 • . The horizontal part of the slab above the transition zone is still observed here. In section 5, the dipping part of the slab disappears and in the map view, we see that the section passes through the gap area between the Izu-Bonin and Mariana plates. In Sects. 8 and 9, we observe the bright images of the Mariana slab which deepens almost vertically. Based on these results, we cannot say definitively if the slab penetrates into the lower mantle or not. In Sect. 9, the slab seems to stop at 700 km depth, whereas in Sect. 8 for P-and S-anomalies, we observe highvelocity anomaly below 700 km depth.

Comparison to previous tomography models
The derived model generally agrees with other previously published regional models (Bijwaard et al., 1998;Widiyantoro, 1999;Van der Hilst and Seno, 1993;Widiyantoro and Van der Hilst, 1996;Gorbatov et al., 2003;Miller et al., 2005Miller et al., , 2006. In Fig. 5, we show an example of comparing of our tomographic images with the model by Miller et al. (2005) which is one of the most recent results for this region. Same as in most of previously published models, our model reveals almost vertical dipping slab beneath the Mariana arc which penetrates into lower mantle. Generally consistent results are also obtained for the Izu-Bonin arc where shallower subduction is observed. However, in details our model has some particular features. First, our tomographic images provide sharper images of the slab than the previous models. Second, we observe more clearly the penetration of a part of the slab beneath Izu-Bonin arc into the lower mantle (Fig. 4,  Sect. 3). Third, tearing between Izu-Bonin and Mariana slab segments down to depths of 500 km is clearly observed in horizontal sections. In contrast to the results by Miller et al. (2005Miller et al. ( , 2006a who observed the gap in the subduction segments beneath Izu-Bonin and northern Mariana arcs, we can see the continuous transition of the subduction style (see the shapes of the fast anomaly between 400-450 km in    5). Finally, we present more different tests to prove the reliability of our results than in any of previous studies for this region.

Verification
As mentioned previously, the data of ISC catalogue are very noisy. Despite the efforts on pre-processing the data and rejection of outliers, we should confess that noise is still important and it can significantly affect the inversion. To estimate the influence of random noise upon the inversion results, we performed the "odd/even" test which consists in independent inversions of two equal subsets of data (e.g., with odd and even numbers of events) using the same inversion parameters as in the case of the entire dataset. If noise effect is significant, it would generate random anomalies which would be different in two independent results; they should be ignored or smoothed by increasing the damping. The results of this test shown in Fig. 6 demonstrate fairly high consistency of the main patterns which confirm their independence of random noise in the data.
To evaluate the spatial resolution of the model, we have performed the synthetic checkerboard test. The synthetic model was defined as alternating positive and negative rectangular anomalies. For the P model the horizontal size of the cell was 1 • × 1 • with an empty space of 0.5 • ; for the S model the corresponding parameters were 2 • × 1.5 • and 1 • . With the depth, the sign of anomalies alternated every 200 km. To produce the synthetic dataset, we used the same rays as in the case of real data inversion. The synthetic times were perturbed with random noise with the statistical distribution which was estimated after analysis of real residuals in the revised ISC catalogue. The magnitude of the noise was 0.5 s and 0.8 s for the P-and S-data, respectively, which lead to the same variance reduction after inversion (about 50 %) as in the real data case. When reconstructing the synthetic model, the values of the inversion parameters were the same as in the case of real data inversion. Figure 7 displays the reconstruction of the checkerboard model at different depths. It shows that in the central part of the area, where the slab is located, the resolution is fair. However, below 700 km in depth the resolution becomes poorer and does not allow the recovering of such small-scale features. An important test consists in reconstruction of a synthetic model with realistic configuration of anomalies. In Fig. 8 we present a model which reproduces the configuration of the slab in horizontal and vertical sections. The synthetic anomalies in this case are free-shaped anomalies along the vertical sections. They have fixed shape in the direction across the section within a certain band. This method allows defining the realistic shapes of different segments of the slab. In this model, the Izu-Bonin slab is gently dipping and becomes horizontal in the transition zone. On the contrary, The Mariana slab consists of three segments and almost vertically dips down to 1200 km. Between the Mariana and Izu-Bonin arcs we define a small gap. As in the case of the checkerboard test, the synthetic data were perturbed with random noise of 0.5 and 0.8 s magnitude for P and S data, respectively. The results of the reconstruction in horizontal and vertical sections show that the main patterns in the mantle are correctly reconstructed. However, the uppermost structures in the Mariana area are strongly biased because of lack of data and tradeoff between velocity structures and source parameters. In the lower mantle, the slab-related anomaly is strongly smeared, especially for the S-model, however it is retrieved at a correct place.

Discussion
The main feature of the obtained model is the slab-related high-velocity anomaly which is clearly observed in all horizontal sections of P and S model, except for the uppermost part where it is not resolved for the Mariana area. In vertical sections we observe complex geometry of the subducting plate which changes the shape along the Izu-Bonin and Mariana arcs. Based on these results in Figs. 9 and 10 we propose two scenarios based on different reference frames which explain the present shape of the slab.

Scenario #1
The first scenario shown in Fig. 9 corresponds to the reference point located in the middle of the Philippine Sea Plate. According to (Seno and Maruyama, 1984;Seno et al., 1993;Hall et al., 1995a, b;Hall, 2002), presently curved shape of the IBM trench (violet line in Fig. 9) was much straighter in the recent past (e.g., 12-25 Ma) as schematically shown with black line in left plot in Fig. 9. An indirect argument for this hypothesis is the nearly linear distribution of high-velocity anomalies at 975 km depth shown in Fig. 9 which probably represent the curvature of the trench in the past (dotted line).
Taking the present rate of subduction (∼5 cm yr −1 ) this depth corresponds to the location of the trench about 20-25 Ma. It can be seen that since that time in different segments of the trench there were both backward and forward displacements. Note that a segment with the forward displacement corresponds to the presence of the Ogasawara large igneous plateau. The change of the subduction shape roughly corresponds to transitions between different displacement styles of the trench. Based on the results of tomographic inversion we propose our back-time reconstruction to about 10-12 Ma.
In sections to the right of Fig. 9, which roughly correspond to vertical Sects. 2, 3, 4 and 9 in Fig. 4, we present the hypothetical evolution of subduction in different segments of IBM arc system. We propose that initially the subduction occurred along a weakly curved trench and the shape of the slab was generally similar throughout the IBM system (namely with the dipping angles of 40-50 degrees). For Section A1-A2 we expect a considerable backward displacement of the trench due to eastward movement of Asia. In this section the oceanic lithosphere seems to be united with a large segment beneath Japan where the slab is relatively shallow and turns to be horizontal in the transition zone to a large extent (Zhao et al., 2009). The reason why this old and dense lithosphere beneath Japan behaves in this way is actively discussed in the literature (e.g., Zhao et al., 2009). It might be due to "roof" shape of the slab which gives additional straight and prevents steep bending of the plate. Moving back of the trench in this area causes decreasing the dipping angle of the slab, as shown in Section A1-A2 in Fig. 9. For Section B1-B2, the situation seems to be significantly different. Here the shallow part of the slab seems to be detached from the horizontal body in the transition zone. In this part of the arc the trench remains at the same position. The continued subduction at a rate of about 4 cm yr −1 (DeMets et al., 2010) causes accumulation of the slab material in the transition zone (stage 2). When reaching a critical volume, the accumulated denser material starts to descend to the lower mantle.
In Section C1-C2 we observe steepening of the slab. In this segment of the arc the trench is though to migrate forward (Schellart et al., 2008, Miller et al., 2006 which is probably caused by the presence of the large igneous Ogasawara Plateau. It is possible that this plateau increased the strength of the lithosphere and made it more buoyant. That is why, when this plateau reached the trench, it stopped the subduction and leaded to displacement of the trench along the incoming plate movement. In this case, the anchored lower part was fixed, and the displacement of the upper part caused the steepening of the slab. The last Section D1-D2 corresponds to the vertical Sect. 9 (Fig. 4) which passes across the Mariana segment. According to Miller et al. (2006b) and authors referred therein, the Mariana trench is moving backward, and initially the situation was similar to Section A1-A2. The Mariana plate was detached from the Isu-Bonin plate by the Ogasawara plateau which served as a "knife" which cut the lithosphere. We propose that this detachment was facilitated by remnant heat which may be kept in plume traces beneath the plateau. The isolated heavy Mariana plate cannot behave as the northernmost part of the Izu-Bonin plate which is supported by the Japanese shallow slab. Even if in stage 1 the subduction had a normal dipping angle (e.g., 45 • ), in stage 2, the shallow part of the slab had to be detached from the deep one, and a new subduction zone originated. The newly formed subduction is mostly driven by high density of the subducted lithosphere which pulls it down. In this case, it behaves not as a normal subduction: the descending of the plate occurs due to both ocean plate movement and trench displacement, as shown in section D1-D2. This is supported by relatively slow rate of the subduction (∼2 cm yr −1 ) and an approximately same rate of the slab retreat in the Mariana trench (DeMets, 2010). The remnant part of the previous subduction zone could descend to the transition zone and then further could continued sinking to the lower mantle. A possible location of sinking of the older subducted slab is in the area, where the resolution of the tomographic model is poor. However, a strongly smeared anomaly of high P-velocity model in the lower mantle can be interpreted as a trace of the previous subduction and it would fit to this hypothesis.

Scenario #2
A weak point of the scenario shown in Fig. 9 is that fixing the reference point in the Philippine Sea Plate needs a very fast eastward movement of Asia. To keep the divergence and convergence balances in subduction and spreading zones, the rate of eastward displacement of Asia should be about 5-6 cm yr −1 , which is similar to the displacement rate of the Pacific Plate. At the same time, much slower relative displacements in continents compared to the fast rates of subduction and spreading processes in oceans suggest that the absolute movements of the entire continents are much slower that those of oceanic plates. Based of this statement, we propose an alternative scheme of plate movement corresponding to the fixed location of Asia (Fig. 10). To estimate the total displacement balance between the Pacific Plate and Asia, we need to consider the entire Philippine Plate including the Ryukyu subduction zone and the Okinawa backarc spreading in the west. The indicated velocities are compiled from different studies including (Suppe J., 1984;Sdrolias et al., 2006;Schellart et al., 2008;Hsu et al., 2012). For the northern part of the IBM system (Section A1-A2 in Fig. 9) we see that in the vertical section the slab penetrates at a high rate (more than 10 cm yr −1 ) between a fixed Asian lithosphere and lower mantle having zero velocity in this reference frame. This movement is only possible in the case of a sufficiently thick layer in the transition zone between the slab and lower mantle playing the role of lubricants. This explains the shallow angle of dipping the slab and its nearly horizontal shape in the left side of the section.
The second section includes two converging (Mariana and Ryukyu subductions) and two diverging zones (Okinawa and Mariana troughs). The rates of shortening and extension in these zones are given as average estimates of several results given in different studies (Suppe J., 1984;Kato et al., 2003;Hsu et al., 2012). The rates of the plates between these zones are presumed to be constant. If we assume that the Asia is fixed in respect to the lower mantle, the entire system including Philippine Sea and Pacific Plates moves westward, as predicted by Schellart et al. (2008). In this case, the fastest Pacific Plate moves at a rate of 12 cm yr −1 . The lower part of the subducting Pacific slab seems to be anchored in the lower mantle. Westward migration of the trench causes steepening of the slab which may lead to vertical and even reverse orientation. A similar scenario has been proposed by Schellart (2011) based on numerical modeling.
The variability of the subduction styles leaded to clear segmentation of the Izu-Bonin and Mariana subduction zones with the boundaries marked with vertical seismicity clusters (Fig. 2). Similarly as in Miller et al. (2006a), between the Mariana and Izu-Bonin segments we can see the gap corresponding to vertical Sect. 6 in Fig. 4 which is observed down to 400 km depth. Initially it was formed by the effect of the Ogasawara plateau. Further opening of this gap occurred due oblique orientation of the Pacific Plate movement in respect to the northern part of the Mariana loop. The strike-slip displacement along the trench in this area causes the further detachment of the upper part of the slab. The lower part of the slab seems to be united and it probably corresponds to the stage when the trench was not perturbed by relative forward and backward displacements. A problem of this model is the existence of deep seismicity in the gap area where the highvelocity slab is not observed. We propose that the gap area is formed by plastic extension of the subducting plate. In this case the resolution of the tomography model might be not sufficient to resolve the thinned part of the slab in the extension zone. Furthermore, the presence of remnant heat beneath the Ogasawara plateau might lead to degrading some physical properties of the slab and to decreasing seismic velocities.

Conclusions
The main result of this study is the new model of seismic velocity heterogeneities of P-and S-waves in the mantle beneath the Izu-Bonin and Mariana arcs which was constructed from tomographic inversion of global travel time data (ISC). This model was thoroughly verified using a series of different tests which showed generally high robustness of the results. The main feature of this result is a slab-related anomaly which is similarly expressed in the P and S models. There is a general correlation with previous results derived by other authors; however in details we have detected some new features. We have found that the slab beneath IBM arc system has strongly variable dipping angle which generally correspond to the loops of the trench shape. We also analyze the distribution of the deep seismicity which was projected to the profile along the trench. The striking feature of the seismicity distribution is that deep seismic events form clear vertical clusters which generally fit to the limits of the IBM arc system segments. We propose that these clusters mark the boundaries between autonomous parts of the slab which may subduct independently with different rates.
We propose the two scenarios based on different reference frames which explain the present shape of subducting plates in the IBM system. The first scenario explains the variable dipping angle of the slab beneath the IBM arc system by relative forward and backward displacements of the trench. In the second scenario, the entire system is moving westward, and the different angles of the slab can be explained by properties of the subducting lithosphere and the existence of rigid and buoyant blocks related to seamounts and intrusive plateaus. These two scenarios present two extreme cases and the truth might appear to be in the middle. Also, we aware about the uncertainties related to determination of relative rates in divergent and convergent zones which may be significantly different in different studies. In addition, we consider the balances of displacements along 2-D profiles, whereas the displacements are clearly 3-D which may also affect our estimates. All these reasons show the importance of implying quantitative estimates based on 3-D numerical modeling. In this case, the tomography results might serve for defining the boundary conditions.