Tomography of the subducting Pacific slab and the 2015 Bonin deepest earthquake (Mw 7.9)

On 30 May 2015 an isolated deep earthquake (~670 km, Mw 7.9) occurred to the west of the Bonin Islands. To clarify its causal mechanism and its relationship to the subducting Pacific slab, we determined a detailed P-wave tomography of the deep earthquake source zone using a large number of arrival-time data. Our results show that this large deep event occurred within the subducting Pacific slab which is penetrating into the lower mantle. In the Izu-Bonin region, the Pacific slab is split at ~28° north latitude, i.e., slightly north of the 2015 deep event hypocenter. In the north the slab becomes stagnant in the mantle transition zone, whereas in the south the slab is directly penetrating into the lower mantle. This deep earthquake was caused by joint effects of several factors, including the Pacific slab’s fast deep subduction, slab tearing, slab thermal variation, stress changes and phase transformations in the slab, and complex interactions between the slab and the ambient mantle.

the subducting Pacific slab clearly, which shed new light on the causal mechanism of the 2015 Bonin deep earthquake and subduction dynamics in the study region.

Data
We used arrival-time data of earthquakes which are selected from the repressed ISC (International Seismological Center) data base 7 , the JMA Unified Earthquake Catalogue 1 , and the Annual Bulletin of Chinese Earthquakes 8 . For selecting a best set of earthquakes, the crust and mantle (0 to 700 km depth) are divided into cubic blocks. The block size is 1° × 1° × 20 km for the whole globe but 0.5° × 0.5° × 10 km for the target Izu-Bonin region (see Supplementary Figure S2). Among the many events located within each block, only one earthquake is selected which was recorded by the largest number of seismic stations and has the smallest error in its hypocenter location. Our data set thus collected contains 5,126,696 arrival times of P, pP, PP, PcP and Pdiff waves from 39,323 earthquakes recorded by 9141 seismic stations in the world (see Supplementary Figure S1). Hence the target Bonin region ( Fig. 2 Figure S2) is well sampled by the up-going and down-going rays of both the direct P waves and later phases 6,9 . (Figs 3 and 4). In Fig. 3d,e, the Pacific slab with deep seismicity shows a lower velocity around 410 km depth, which is an artifact due to a lower resolution of the tomographic image there. The active arc volcanoes are located above low-velocity (low-V) zones in the mantle wedge above the subducting slabs. The low-V zones represent the source areas of arc magmatism and volcanism caused by a combination of fluids from the slab dehydration and corner flow in the mantle wedge driven by the plate subduction e.g. refs 11-14. The subducting Philippine Sea slab is also imaged clearly to west of the Ryukyu Trench (Figs 3 and 5a and Supplementary Figure S3), indicating that the Philippine Sea slab has also reached the MTZ depth, being consistent with previous local and regional tomographic models e.g. refs 13, 15-17. In the vertical cross-sections north of ~28°N latitude ( Figure S4), which are discussed in the next section.

Results
The 2015 Bonin deep event is relocated using our 3-D Vp model. The relocated focal depth is 667.2 ± 0.5 km (see Supplementary Table S1), and the relocated hypocenter is located within the high-V slab but close to the eastern boundary of the near-vertical slab (Fig. 5). Note that the hypocenter is just the initial point of the earthquake rupture. Waveform inversions show that the average source depth of this great deep event is ~680 km and its overall source dimension is ~40 km in a shallowly-dipping fault plane 2 .
Our whole-mantle tomographic images beneath the Izu-Bonin region are shown in Supplementary Figure S5. Below the Pacific slab in the MTZ, intermittent high-V zones appear in the lower mantle, and broad high-V anomalies exist above the core-mantle boundary, which reflect old pieces of the Pacific slab that have collapsed down to the lower mantle and reached the core-mantle boundary due to gravitational instability caused by phase transformations e.g. refs 6 and 21.
Detailed resolution analyses are made to confirm the main features of the tomographic results (see Supplemen tary Information Figures S6-S11). The results of these resolution tests show that main features of the tomographic results (Figs 3-5), in particular, those in and around the source zone of the 2015 Bonin deep earthquake, are quite robust.

Discussion
Based on the present results and also previous seismic studies of the Izu-Bonin region e.g. refs 16-23, we propose a model on the subducting Pacific slab and the source location of the 2015 Bonin deep earthquake in the slab (Fig. 6). The Pacific slab is split at ~28° N latitude, i.e., slightly north of the hypocenter of the 2015 deep event which occurred at 27.7° N latitude (see Supplementary Table S1). In the north, the slab becomes stagnant in the MTZ, whereas in the south the slab is directly penetrating into the lower mantle. Miller et al. 18,19 showed that the Pacific slab has an anomaly of physical properties in the upper mantle, which is located in an area north of the 2015 Bonin hypocenter. They found that most of the intermediate-depth and deep earthquakes occurring within this anomaly in the slab had lateral tension mechanisms, potentially related to the slab tearing. The 2015 Bonin deep event also had a lateral tension mechanism, which may be related to the slab tear 3 .
Our model (Fig. 6) is similar to the second model of Ye et al. 2 . In other words, our tomographic results do not support their first model but their second model, because a simple vertical high-V zone is clearly visible in the 2015 source area (Fig. 5b), which reflects the Pacific slab penetrating into the lower mantle. We agree with Takemura et al. 3 on their point that the 2015 Bonin deep event occurred within the slab and close to the lower slab boundary, but we do not prefer their model which shows a flat slab in the MTZ at the 2015 Bonin hypocenter. The slab pileup model 4 is similar to the first model of Ye et al. 2 and is even more complex, which is quite different from  our model (Fig. 6) and our tomographic images (Fig. 5) which show a simple slab geometry without buckling. The receiver-function image 4 might be over-interpreted, which was obtained beneath only one seismic station. The conversion points in and below the MTZ in the receiver-function image 4 may not represent boundaries of a buckled slab but complex phase changes in the slab at the MTZ depths.
To date, several studies of receiver-functions and waveform modeling have shown that the 670-km seismic discontinuity is depressed down to ~690 km depth beneath the Bonin region 4,[24][25][26] . The focal depth of the 2015 Bonin deep event was estimated to be in a range of 664-682 km (see Supplementary Table S1). Hence, it is almost certain that this deep event occurred within the MTZ. The slab-related high-V zone looks thicker at depths of 600-800 (Fig. 5b), which may suggest that the slab has thickened near the MTZ bottom due to resistance at the upper-lower mantle boundary e.g. refs 6, 21 and 27.
Previous 3-D velocity models of the Izu-Bonin-Mariana region are obtained by using global tomography or large-scale regional tomography e.g. refs 6, 16-22, which generally have a lower resolution because large blocks or grid intervals were adopted in the 3-D model parameterization. Because there are few seismic stations in this region, the tomographic models are sensitive to the data set, model parameterization, inversion algorithm, and damping and smoothing parameters adopted in the inversion. Hence the relationship between the 2015 Bonin deep event and the subducting slab is not clear in the previous tomographic models (e.g., Supplementary Figure S4). In contrast, our present study focuses on the source zone of the 2015 Bonin deep event, and we use a much better data set and adopt a dense 3-D grid in the target region, hence our present tomographic model shows a clear relation between the slab and the deep event (Fig. 5). To better resolve the detailed structure of the subducting slab in the MTZ, P and S waves reflected and/or converted at the MTZ discontinuities e.g. refs 26, 28-30 should be collected and used in tomographic imaging, in addition to deploying a dense network of ocean bottom seismometers in the study region.
Many previous studies have revealed tears in subducting slabs in the upper mantle and slab tear-related intermediate-depth seismicity in many regions e.g. refs 31 and 32. In addition to the Izu-Bonin region, a slab tear in the MTZ was detected beneath western Japan 33,34 . Slab buckling and folding in the MTZ have been observed in several regions, which provide an important control on the distributions and focal mechanisms of deep-focus earthquakes e.g. ref. 35.
The physical processes that permit the occurrence of deep earthquakes are still not well understood 36 . Several mechanisms for deep earthquakes have been proposed, including transformational faulting triggered by metastable olivine transforming to spinel in the cold, stressed core of the subducting slab e.g. refs 27, 37-40, thermal instability and run-away shear melting e.g. refs 41-43, and dehydration embrittlement e.g. refs 44-46. All of these mechanisms depend on the thermal structure of deep slabs and the deviatoric stress conditions associated with the slabs impinging on the upper-lower mantle boundary e.g. refs 2 and 43. It is challenging to distinguish between these possible mechanisms for deep earthquakes, because of the difficulty to clarify their source dimensions and rupture processes, as well as the fine slab structure in and around the source zones of deep earthquakes. Among these mechanisms, dehydration embrittlement is considered to operate for the intermediate-depth earthquakes but may not for the deep earthquakes e.g. refs 27 and 36. The viability of transformational faulting as a mechanism for deep earthquakes hinges on the presence of a sufficient amount of metastable phase 36 . So far, a metastable olivine wedge (MOW) has been detected within the subducting Pacific slab at the MTZ depths beneath Southwest Japan 47,48 ,Mariana 49,50 , and the Japan Sea 51 , and it is considered that the generation of deep earthquakes in these regions is related to the existence of MOW in the slab. Hence, transformational faulting is a conventional and popular mechanism for deep earthquakes. However, our present tomography could not image the MOW in the Pacific slab, due to the limited resolution of the tomographic model and because the MOW, if it exists, should be very thin at the great depth. Waveform inversions suggest that localized stress concentration associated with the pronounced deformation of the Izu-Bonin slab and proximity to the 670-km phase transition may have played a dominant role in generating this significant earthquake 2 .
Meng et al. 52 suggested that the 2013 Okhotsk deep earthquake (Mw 8.3, 610 km depth) was affected by thermal thinning of the Pacific slab because it occurred near the northern end of the slab as revealed by seismic tomography 53 . We think that the generation of the 2015 Bonin deep earthquake could be also affected by thermal variation of the slab, because this event occurred near the northern edge of the near-vertical segment of the Pacific slab (Fig. 6) where the slab must be heated by the hot ambient mantle.
Based on the above discussions, we deem that the generation of the 2015 Bonin deep earthquake was caused by joint effects of several factors, including the fast deep subduction of the Pacific slab, slab tearing and its related thermal variation, stress changes and phase transformations in the slab near the upper-lower mantle boundary, as well as complex interactions between the subducting slab and the ambient mantle. It is not clear why the isolated deep earthquake is so large, but its large size indicates the presence of a large volume of material still in seismogenic conditions 54,55 . There may be very infrequent mineral transformation or volatile release processes that occur only under particularly high deviatoric stress conditions allowing large dynamic stress relaxations to take place 2 .

Methods
We conducted tomographic inversions for 3-D Vp structure using a modified version of the multiscale tomographic method 5,6 . For expressing the 3-D Vp structure, a denser 3-D grid with a lateral grid interval of ~50 km is arranged in a depth range of 0-1000 km beneath the target Izu-Bonin region including the 2015 Bonin deep earthquake, whereas a coarse grid with a lateral grid interval of ~220 km is arranged in the whole crust and mantle of the Earth (see Supplementary Figure S2). Vp perturbations at every grid nodes from the one-dimensional iasp91 velocity model 56 are taken to be unknown parameters. The Vp perturbation at any location in the crust and mantle is computed by a linear interpolation of the Vp perturbations at the 8 grid nodes adjacent to that location. A 3-D ray tracing technique 11,21 is adopted to calculate theoretical travel times and ray paths. The LSQR algorithm 57 with smoothing and damping regularizations is applied to solve the large and sparse system of observational equations 5,11 . We conducted many tomographic inversions of our data set to search for the optimal smoothing and damping parameters considering the balance between the travel-time residual reduction and the norm of the obtained 3-D Vp model 5,6,9 .