Slab morphology and deformation beneath Izu-Bonin

Seismic tomography provides unique constraints on the morphology, the deformation, and (indirectly) the rheology of subducting slabs. We use teleseismic double-difference P-wave tomography to image with unprecedented clarity the structural complexity of the Izu-Bonin slab. We resolve a tear in the slab in the mantle transition zone (MTZ) between 26.5° N and 28° N. North of the tear, the slab is folded in the MTZ. Immediately above the fold hinge, a zone of reduced P-wavespeed may result from viscous dissipation within an incipient shear zone. To the south of the tear, the slab overturns and lies flat at the base of the MTZ. The ~680 km deep 2015 Bonin earthquake (Mw~7.9) is located at the northernmost edge of the overturning part of the slab. The localised tearing, shearing and buckling of the Izu-Bonin slab indicates that it remains highly viscous throughout the upper mantle and transition zone.

The regional model in the tomographic inversion is bounded by the 8º N and 36º N lines of latitude, 134º E and 149º E lines of longitude, and extends from the surface to 900 km depth (Supplementary Figure 3). We adopted a two-step strategy for tomographic inversion, in which a coarser inversion grid is first used and the resulting velocity model is used as the starting model for an inversion using a finer inversion grid. In this way, the final inverted velocity model contains both large-scale and small-scale features. For this study, the inversion grid nodes in longitude and depth are fixed but in latitude we first use coarser grid and then finer grid. In depth, grid nodes are set at 0, 80, 150, 210, 270, 320, 370, 410, 440, 480, 520, 560, 600, 630, 660, 690, 740, 810, and 900 km (Supplementary Figure 3). In longitude, the inversion grid intervals are 0.3° where the 2015 Bonin earthquake happened but are larger near the edge of the regional model.  Figure 3). In latitude, the grid nodes for the coarser grid are set with an interval of 3°-4° and are located at 8.0º, 12.0º, 16.0º, 20.0º, 23.5º, 27.0º, 30.0º, 33.0º, and 36.0º N (Supplementary Figure 3). The finer grid nodes in the second inversion have an interval of 1º between 25º and 30º N and have a varying interval of 1.5º to 2º outside this region. We use the spherically symmetric global model ak135 2 as the initial regional velocity model. The regional model is embedded in a coarser (5° × 5°) global model that extends to 2889 km depth with an average depth interval of ~180 km (Supplementary Figure 3). We use the 3-D model LLNL_G3Dv3 3 for the initial global velocity model.
For the teletomoDD algorithm, large residuals are down-weighted by a bi-weight function [4][5][6] . A first-order smoothing constraint is applied to slowness perturbations to stabilize the tomographic system. The damped least squares system is solved by the LSQR algorithm 7 .
The pseudo-bending ray tracing algorithm 8 that is extended to spherical coordinates 9 is used to calculate travel times. The damping and smoothing parameters are chosen via a trade-off analysis of model perturbations and travel time residuals. For the coarser grid inversion, the root mean square (RMS) value of differential time residuals decreases from 1.089 s to 0.058 s after four simultaneous and two location-only iterations. In the first two simultaneous iterations, weights on absolute arrival time are set larger than differential arrival times to derive the large-scale structure. In the last two simultaneous ones, weights on differential arrival times are set much larger than the absolute arrival times to refine the near source structure. Supplementary Figure 4 shows cross sections of the inverted Vp model along six profiles in Figure 1.
In the second stage of the inversion, the initial velocity structure on the finer grid corresponds to the inverted Vp model using the coarser inversion grid. As for the coarser grid inversion, four simultaneous and two location-only iterations were performed and final RMS value of differential travel time residuals is 0.049 s, slightly lower than the coarser grid inversion. Figure 2 shows cross sections of the final Vp model using finer inversion grid along the 6 profiles in Figure 1. In comparison, the two models look very similar but the model from the 2 nd stage of inversion with finer grid has larger anomaly amplitudes ( Figure 2 and Supplementary Figure 4).
To evaluate the resolution of the Vp model around the subducting slab, we conducted checkerboard resolution tests 10 . Positive and negative 5% velocity anomalies are added to the 1D regional velocity model used in the inversion at alternating grid nodes to create the checkerboard velocity model. Synthetic absolute and differential travel times were calculated from the checkerboard velocity model using the same source-receiver geometry as the real data. Synthetic data were subsequently fed into the inversion system, with similar inversion parameters as the real data inversion. Supplementary Figures 5 and 6 show recovered checkerboard patterns for the coarser inversion grid at different depths and latitudes.

Supplementary Figures 7 and 8 show recovered checkerboard patterns for the finer inversion
grid. The checkerboard tests suggest that the inverted model is well-resolved horizontally and vertically and that the high-velocity structure of slab in the Izu-Bonin area is well-defined below 100 km for the selected inversion grid. We also tried using finer inversion grid nodes in latitude with an interval of 1° from 8° N to 36° N. As shown by the checkerboard test results, the resolution becomes degraded compared to the inversion grid used for real data inversion (Supplementary Figure 9). For this reason, we think the currently selected inversion grid is a good compromise between resolution and robust retrieval of velocity anomalies. To quantitatively assess the resolvability of checkerboard models, we calculate the semblance value at each grid node between true and recovered checkerboard models using nearest nodes in three directions following the way of Guo et al. 11 which is modified from the method of Zelt 12 . We calculated the model resolvability for both coarser and finer grids. Supplementary  Figure 4). For the model inverted using the finer inversion grid, we combine resolvability values obtained for coarser and finer grids to determine the well-resolved model region. This is because the model using finer inversion grid is inverted using coarser inversion grid model as the starting model so that any artifact from the coarser grid model would be brought onto the finer grid model. The model at finer grid node is treated as well resolved if the resolvability with the finer inversion grid is larger than 0.75 and the resolvability with the coarser inversion grid is larger than 0.8 ( Figure 2).
In addition to the checkerboard resolution tests, we also conducted additional resolution tests by constructing various synthetic models with different characteristic features based on the final inverted velocity model. Supplementary Figure 12 shows the synthetic model reflecting the main features of slab morphology resolved by the real data inversion, such as slab overturn to the south of 27° N (Figure 2 and Supplementary Figure 4). Here the synthetic model is constructed following high velocity anomalies of the inverted model. As a result, the "slab" in the synthetic model does not have smooth edges (Supplementary Figure 12). As for the checkerboard resolution test, we calculated absolute and differential travel times with the same distribution of sources and receivers as the real data. Then we applied the same inversion scheme as the real data for the synthetic data. It can be seen that the main features can be well recovered, but with a slightly lower amplitude (Supplementary Figure 13). This is expected because of the damping and smoothing regularizations used to stabilize the inversion system. This also indicates that the velocity anomalies in the inverted Vp model (Figure 2) are likely underestimated.
To further test the robustness of the low velocity anomaly around 400 km in depth with respect to the shallower and deeper parts of the slab appearing in the inverted model ( Figure   2), we constructed a synthetic slab model similar to that in Supplementary Figure 12 but containing a low velocity zone around depth 400 km (Supplementary Figure 14). The inverted model shows that not only slab overturn can be well resolved but also the low velocity zone around 400 km can also be resolved (Supplementary Figure 15). Two other endmember slab models further show the robustness of main slab features resolved in the real data inversion (Supplementary Figures 16-19). The model in Supplementary Figure 17 shows that slab overturn cannot be retrieved as an inversion artifact. The model in Supplementary   Figure 19 shows that the inversion procedure can also recover the case where the slab overturns throughout the area. In the slab model shown in Supplementary Figure 16, we include the (relatively low) velocity anomaly in the slab around depth 400 km. This feature is also well resolved (Supplementary Figure 17). These synthetic tests show that the velocity anomalies that we discuss in the main text are robust and can be resolved by the real data ( Figure 2).
The last resolution test that we conducted is the restoration test 13 , which is a traditional and widely used method for model resolution estimation. It can directly evaluate the reliability of the shape and amplitude of velocity anomalies from the real data inversion as well as some features seen from earthquake relocations. Firstly, we set earthquake relocations and velocity models obtained from the real data inversion as "true" earthquake locations and "true" velocity models (Figure 2). The synthetic travel times are then calculated based on these "true" earthquake locations and velocity models. We then add Gaussian distributed noise of mean of zero and standard deviation of 1 second to the synthetic travel times. The same inversion strategy is applied to these synthetic travel times as the real data. It can be seen that the recovered model (Supplementary   Stixrude and Lithgow--Bertelloni 16 ). Therefore, they are expected to result in an increase in slab velocities at a given pressure relative to the ambient mantle. The presence of metastable olivine within a cold slab could lower velocities relative to surrounding mantle at depths greater than the "410 km" transition 17 (Supplementary Figure 23).
With respect to the proposed shear zone in the Izu--Bonin slab, there is little change in predicted slab temperatures along the slab based on plate age and subduction velocity alone 19 . Therefore, major changes in velocity structure along the slab are unexpected. In particular, the ability to metastably preserve olivine is unlikely to be very different in the south relative to the north. We therefore prefer to explain the localised reduction of velocities about the tight fold in the slab as a consequence of deformation, rather than differences in phase assemblage resulting from different thermal conditions. Table 1: Flow law parameters for the shear heating model. Figure 22: Estimates of temperature increase and viscosity within a 25 km thick shear zone at 13 GPa (∼400 km depth). The flow laws used in the calculations correspond to experimental values for dry olivine 14,15 .

Supplementary Figure 23:
The effect of slab temperature and metastability on seismic anomalies relative to ambient mantle. a) Simplified slab temperature structure (in Kelvin) calculated according to the model in Frohlich 18 . The age at the trench is 130 Myr, the velocity is 7 cm/yr, the diffusivity 5e--7 m 2 /s, and the slab dip is 60 degrees. Black contours correspond to depths in km. Temperatures are obtained from analytically-calculated potential temperatures by adding a component due to adiabatic heating, using a PerpleX lookup table for pyrolite (the choice of depleted or undepleted slab material has little influence on the results). b) Equilibrium P--wavespeeds for a pyrolitic slab in km/s. c) P--wavespeeds for a model where olivine remains metastable below 973 K. For ease of calculations, a simplified harzburgite composition is used wherever the slab temperature is <973 K. This is obviously non--physical (it implies that slab material "magically" transforms from harzburgite to pyrolite at 973 K), but simplifying the bulk composition simplifies the definition of metastability (allowing Fe--Mg exchange in olivine and pyroxene, but inhibiting transformation to wadsleyite, ringwoodite and majoritic garnet) while having only a small effect on slab velocities (see the 1% difference in velocities between (b) and (c) at <300 km). Metastable olivine is visible as a wedge of anomalously low velocities between 400--600 km depth. For both (b) and (c), velocities have been smoothed with a Gaussian filter with sigma= 5 km.