Introduction

Double seismic zones, where intraslab earthquakes occur in two distinct dipping planes at depths of ~50–~250 km, are observed widely in subduction zones regardless of the oceanic plate age, convergence rate, or stress orientation1,2,3. Proposed hypotheses for the physical mechanism of intermediate-depth (50–250 km) earthquakes include transformational faulting4, thermal shear instability5, 6, and dehydration embrittlement7, 8. Currently, the hypothesis of dehydration embrittlement is considered to be the leading mechanism, which suggests that earthquakes in double seismic zones are linked to brittle failure associated with dehydration reactions of hydrous minerals in the slab crust and uppermost mantle9,10,11,12. This popular hypothesis has been tested mainly by thermal-petrologic models and laboratory deformation experiments for hydrous minerals in subducting slabs11, 13. Although the existence of hydrous minerals in subducting slabs has been demonstrated by studies of seismic tomography14 and receiver functions15, there still remain questions where and how water migrates in a subduction zone.

Seismic anisotropy in the upper mantle results from crystallographic preferred orientations (CPOs) of minerals (prominently olivine) induced by dislocation creep16,17,18. Laboratory experiments have revealed that the varying types of olivine fabrics (e.g., A-, B-, C-, D-, and E-type), which relate to different olivine slip systems, cause different seismic anisotropy structures due to water content, temperature and stress/pressure17, 18. For example, as decribed by Karato et al.17, A- and D-type olivine fabrics generally exist in the oceanic and continental lithosphere with poor water content; B-type olivine fabrics are induced mainly by increased water content at high stress and cause azimuthal anisotropy normal to mantle flow, which is proposed to explain trench-parallel azimuthal anisotropy in the fore-arc mantle wedge of subduction zones; C-type olivine fabrics are induced by increased water content at lower stress than that of B-type olivine fabrics, and possibly cause negative radial anisotropy (i.e., vertical velocity being faster than horizontal one) in horizontal mantle flow and weak positive radial anisotropy (i.e., horizontal velocity being faster than vertical one) in vertical cylindrical flow; E-type olivine fabrics are induced by moderate water content and observed in island arc environments19. Hence, seismic anisotropy may provide important insights into intermediate-depth seismicity because both of them are closely associated with water migration and distribution in subduction zones.

The subduction processes of the old Pacific plate and the young Philippine Sea (PHS) plate (Fig. 1) control seismotectonics of the Japan Islands. The Japan subduction zone is an ideal place to study the physical mechanism of intermediate-depth intraslab seismicity, because the earthquake distribution exhibits significant changes during the subduction processes (Supplementary Fig. S1). In the old Pacific slab, intraslab seismicity occurs down to the mantle transition zone20. In the young PHS slab, however, the intraslab seismicity occurs down to ~60 km depth and the aseismic PHS slab exists down to ~150 km depth under Chubu, Kinki and Shikoku20,21,22. Beneath Kyushu, the intraslab seismicity occurs down to ~200 km depth, whereas the aseismic PHS slab is revealed down to ~500 km depth20, 21, 23.

Figure 1
figure 1

Map of the Japan Islands and surrounding regions. (a) Tectonic settings. The red triangles denote active arc volcanoes. (b) Distribution of 1656 seismic stations (black triangles) used in this study, which belong to the Japanese national universities, Japan Meteorological Agency (JMA), and the High-sensitivity seismic network (Hi-net). The red lines show locations of the four profiles shown in Fig. 3. (c) Epicentral distribution of 7350 local earthquakes (color circles) used in this study. The colors denote the focal depths. (d) Epicentral distribution of the 376 teleseismic events (red crosses) used in this study. This Figure is generated using GMT 4.5.3 (http://www.soest.hawaii.edu/gmt/) developed by Wessel and Smith59.

Full size image

Previous studies have revealed significant seismic anisotropy in the upper mantle at depths <150 km beneath Japan24,25,26,27, however, little is known about seismic anisotropy in the deep upper mantle at depths >~200 km. The existence of seismic anisotropy in the deep upper mantle is controversial16. The deformed materials in the deep upper mantle are suggested to be isotropic because diffusion creep rather than dislocation creep is dominant there, which may not develop preferred mineral orientation28. Although the mechanism causing seismic anisotropy in the deep upper mantle is still in debate, anisotropic structures may exist at least down to the mantle transition zone29 and the lower mantle30.

We have developed a tomographic method for determining 3-D distribution of P-wave velocity (Vp) anisotropy27, in which Vp radial anisotropy is estimated with assumed vertical hexagonal symmetry by adding one anisotropic parameter. In this study, we apply this method to obtain the first model of 3-D Vp radial anisotropy in the upper mantle beneath the entire Japan Islands using a large number of high-quality arrival-time data of local earthquakes and teleseismic events. Our present results shed new light on the relation between seismic anisotropy and water distribution in the upper mantle.

Results

We used 901,726 P-wave arrival times of 7350 local earthquakes and 39,175 P-wave relative travel-time residuals of 376 teleseismic events recorded by 1656 stations (Fig. 1b–d) during a period of 1998–2012. The relative travel-time residuals of teleseismic events are calculated from waveforms using the multichannel cross-correlation technique31. A slightly modified J-B model32 with varying depths of the Conrad and Moho discontinuities33,34,35 is adopted as the starting Vp model (Supplementary Fig. S2). A 3-D grid with a grid interval of 0.5° in latitude, 0.4° in longitude and at depths of 8, 20, 40, 60, 80, 120, 160, 200, 240, 280, 320 and 380 km is set up in the modeling space to express the 3-D isotropic Vp variation and Vp radial anisotropy. The maximum ray-azimuthal gap angle36 (MRAGA), which denotes the azimuthal coverage of ray paths around a grid node (Supplementary Fig. S3), is set to be 45° during inversions. The optimal values of damping and smoothing parameters for both the isotropic Vp variation and Vp radial anisotropy are obtained by considering the balance between reduction of the root-mean-square travel-time residual and smoothness of the inverted 3-D Vp model (see Supplementary Fig. S4). The local earthquakes are relocated during the inversion process.

Figures 2 and 3 show the obtained images of Vp isotropic variation and Vp radial anisotropy at eight selected depths and along four vertical cross-sections, respectively (more images are shown in the Supplementary Figs S5 and S6). The subducting Pacific slab is revealed as a high-velocity (high-V) zone down to the mantle transition zone, and the subducting PHS slab is imaged as a high-V zone with a complex geometry beneath Southwest Japan. Low-velocity (low-V) anomalies appear in the mantle wedge beneath the volcanic front. The estimated upper boundary of the Pacific slab is generally consistent with the previous results denoted as the red dashed lines in Figs 2 and 3, which was estimated from hypocentral locations of intermediate-depth earthquakes37 and converted waves at the slab upper boundary38.

Figure 2
figure 2

Tomographic images at selected depths beneath Japan. (a) Isotropic P-wave velocity images. The red and blue colors denote low and high velocities, respectively. (b) Images of P-wave radial anisotropy. The red and blue colors denote negative and positive radial anisotropies, respectively. The scales for the isotropic Vp anomaly and radial anisotropy amplitude are shown on the right. A 3-D grid with a grid interval of 0.5° in latitude, 0.4° in longitude and at depths of 8, 20, 40, 60, 80, 120, 160, 200, 240, 280, 320, and 380 km is set up in the study region for the tomographic inversion. The red dashed line in each map denotes the location of the upper boundary of the subducting Pacific slab at each depth37, 38. The red triangles denote active arc volcanoes. The black dots denote the seismicity during a period of 2002–2007 that occurred within a 10 km depth of each layer. This Figure is generated using GMT 4.5.3 (http://www.soest.hawaii.edu/gmt/) developed by Wessel and Smith59.

Full size image
Figure 3
figure 3

Vertical cross-sections of tomography along the four profiles shown in Fig. 1b. (a) Isotropic P-wave velocity images. (b) Images of P-wave radial anisotropy. The definitions of the color scales for the isotropic velocity anomaly and radial anisotropy amplitude are the same as those in Fig. 2. The red dashed lines show the Moho discontinuity33, 35 and the upper boundary of the subducting Pacific slab37, 38, respectively. The red triangles and white dots denote active arc volcanoes and seismicity during a period of 2002–2007 within a 10 km width of each profile. This Figure is generated using GMT 4.5.3 (http://www.soest.hawaii.edu/gmt/) developed by Wessel and Smith59.

Full size image

We conducted a checkerboard resolution test (CRT) and two synthetic tests for the radial anisotropy tomography. In the input model of the CRT, isotropic Vp anomaly and Vp radial anisotropy with an amplitude of ±4% were alternatively assigned to the grid nodes. Random errors in a normal distribution with a standard deviation of 0.15 s were added to the theoretical arrival times. The test results (Supplementary Figs S7 and S8) indicate that the resolution is very good at depths <380 km. In the first synthetic test, the inversion results (Supplementary Figs S5 and S6) were adopted as the input model, and random errors in a normal distribution with a standard deviation of 0.15 s were added to the synthetic arrival times. The second synthetic test is similar to the first one, in the input model the pattern of Vp variations is the same but the anisotropy is opposite. The output results of the synthetic tests (Supplementary Figs S9–S14) indicate that both the isotropic Vp variation and Vp radial anisotropy are well resolved in the upper mantle, suggesting that main features of our tomographic results are reliable.

Compared with our previous study27, which shows a good resolution at depths <150 km using only P-wave arrival times of local earthquakes, our present results reveal complex seismic anisotropy structures in the deeper upper mantle, and a close relation between the intraslab seismicity and radial anisotropy. The seismic PHS slab exhibits a positive radial anisotropy (i.e., horizontal Vp being faster than vertical Vp), but the aseismic PHS slab shows a negative radial anisotropy (i.e., vertical Vp being faster than horizontal Vp) beneath Kinki and Kyushu (see Fig. 2 and vertical cross-sections c and d in Fig. 3). The intermediate-depth earthquakes in the Pacific slab occurred mainly in areas of positive radial anisotropy (Fig. 2 and cross-sections a and b in Fig. 3).

Discussion

Although some small amounts of hydrous minerals (serpentine, antigorite, talc, chlorite, clinoclore, etc.) in the subducting slab are very anisotropic for both P and S waves39, we think that seismic anisotropy is induced mainly by olivine CPOs in the upper mantle considering the resolution scale (~40 km × 50 km × 50 km) of our tomographic model. Recent experiments have suggested that the same olivine fabric shows the same seismic anisotropy for S-wave17 and P-wave40. Since the olivine fabric changes due to different water content, temperature and stress/pressure, we interpret our Vp radial anisotropy results by following and extrapolating the experiment measurements on olivine crystal-fabrics17, 40. Similar to the straightforward method by Michibayashi et al.40, we use the structural framework (x-y-z axes) to denote P-wave velocities in three directions of horizontal flow: Vx denotes P-wave velocity in the mantle flow direction (trench-normal direction), Vy denotes P-wave velocity perpendicular to the flow direction in the horizontal plane (trench-parallel direction), and Vz denotes P-wave velocity in the vertical direction.

When the oceanic plate is produced (marked by 1 in Fig. 4), its lithospheric mantle is anhydrous because water prefers to dissolve in melting. Hence, the olivine fabric of anhydrous lithospheric mantle is A-type, in which seismic anisotropy is Vx > Vy > Vz, that is, the fast Vp direction is trench-normal in azimuthal anisotropy, and positive radial anisotropy occurs. This explanation is consistent with the observation that Pn waves travel faster normal to the mid-ocean ridge41.

Figure 4
figure 4

Schematic vertical cross-sections depicting relations between seismic anisotropy, water migration, and double seismic zone in cold (a) and warm (b) subduction zones. They show the anisotropic structures of the subducting slabs (marked by 1, 2 and 4), the overlying B-type mantle wedge, and the hot mantle upwelling associated with arc volcanism (marked by 3). The dashed arrow denotes the slab-driven corner flow direction. The solid arrow denotes water migration induced by dehydration embrittlement of hydrous materials in the slabs. Possible types of olivine fabrics in different parts are also shown. In each of the big circles below the vertical cross-sections, the solid square, open triangle and open circle denote the directions of the maximum, intermediate and minimum P-wave velocities, respectively. The two open triangles in the circle ③ denote that the axis directions and amounts of the intermediate and minimum P-wave velocities cannot be distinguished by this study. The plus (+) and minus (−) symbols beside the big circles denote positive and negative radial anisotropies, respectively. The velocity features for the types of olivine fabric are derived from the previous studies17, 40. This figure is generated by free and open source INKSCAPE 0.91 (https://inkscape.org/en/).

Full size image

Large amounts of water entered the oceanic lithospheric mantle from the outer-rise normal faults42 and was stored in hydrous minerals when the oceanic plate began to subduct. The dehydration of hydrous minerals in the slab crust and uppermost mantle at intermediate depths may cause brittle failures, which result in double seismic zones within the subducting slab. As shown in Fig. 3, the areas of positive radial anisotropy exist roughly above the high-V zones (the subducting slabs), probably suggesting that the water has been expelled from the subducted slabs into the overlying mantle wedge. The expelled water hydrates the surrounding slab mantle and the overlying mantle wedge, resulting in B-type olivine fabric there (marked by 2 in Fig. 4). Thus the seismic anisotropy is Vy > Vx > Vz, that is, trench-parallel azimuthal anisotropy and positive radial anisotropy in the upper part of the subducting slabs and the overlying mantle wedge. This result agrees well with the observed trench–parallel azimuthal anisotropy in the fore-arc mantle wedge and the uppermost parts of the slab mantle beneath Japan by shear-wave splitting measurements24, 43 and P-wave anisotropic tomography27. However, the conditions inducing the B-type olivine fabrics or serpentinization are still unclear when the water migrates into the overlying mantle wedge, for example, the CPOs of serpentine rather than B-type olivine in a hydrated mantle wedge are proposed to explain a strong trench-parallel anisotropy (1–2 s shear-wave delay time) beneath the Ryukyu arc26. Since the anisotropic amplitudes beneath Japan24, 27, 43 are smaller than those beneath the Ryukyu arc, the anisotropic features of the hydrated slab mantle (marked by 2 in Fig. 4) and the overlying mantle wedge are induced mainly by B-type olivine fabrics.

The water, migrated from the hydrated overlying mantle wedge44, mainly entered melts rather than olivine and induced the mantle upwelling beneath the arc volcanoes. Thus, E-type olivine may be dominant in the mantle upwelling, resulting in negative radial anisotropy, because Vx is the maximum P-wave velocity in the vertical flow direction (marked by 3 in Fig. 4). The result is consistent with the observed E-type olivine fabrics of some samples in the island arc environments19.

Hydrated materials of the oceanic plate subsequently release fluids induced by the rising temperature of slab during its subduction45. The hydrous phases may be completely dehydrated at some critical depth (~200 km), and hydrogen could not be transferred to a greater depth when the slab temperature is relatively high, however, a cold slab may transfer hydrogen even into the lower mantle with staged dehydrations during the subduction process39. Because the PHS slab is relatively young (~10–40 m.y.)46, it may be anhydrous and could not produce seismicity due to the absence of dehydration embrittlement at its intermediate depths. The olivine fabrics in the anhydrous aseismic PHS slab may change to D-type because of high stress induced by slab bending, for which seismic anisotropy is Vx > Vy > Vz. In a number of subduction zones on Earth, the slab geometry generally shows an increasing dip angle at a greater depth in the mantle, possibly affected by trench retreat and slab rollback47. Thus, the above-mentioned x and z axes for the olivine fabric in the slab rotated counterclockwise and replaced each other with unchanged y axis, and Vx is the velocity in the vertical direction while Vz and Vy are the velocities in the horizontal plane. As a result, the aseismic PHS slab shows negative radial anisotropy (marked by 4 in Fig. 4b). The same explanation could be applied to the lower part of the cold Pacific slab (marked by 4 in Fig. 4a) at depths of ~240−~320 km where negative radial anisotropy is dominant and induced by D-type olivine fabrics. Our results also show that the intraslab seismicity mainly occurs in areas of positive radial anisotropy (Figs 2 and S5) even down to 380 km depth in the cold Pacific slab. Although it is suggested that the cold Pacific slab can transfer hydrogen into the lower mantle, our results cannot resolve the mechanism of deep earthquakes (>250 km depth) because little is known about the water content and dehydration reaction in the subducting slab in the deeper upper mantle and at the mantle transition zone depth48. It is still in question whether the deeper intraslab earthquakes could be generated by dehydration embrittlement, because the water storage capacity for hydrous phases in the cold slab would increase at greater depths9, 49, 50. Little correlation is found between seismicity rate and calculated slab dehydration flux at depths >240 km51. The hypothesis of transformational faulting on the upper and lower edges of a metastable olivine wedge is applied to explain the seismicity in the deeper double seismic zones observed in Tonga at depths of 350–460 km52 and in Izu-Bonin at depths of 300–450 km53.

Conclusions

Seismic anisotropy provides a wealth of new information and direct constraints for understanding geodynamics of the Earth’s interior. We present the first images of 3-D high-resolution P-wave radial anisotropy in the upper mantle beneath the entire Japan Islands. Our results reveal a close relationship between anisotropy and deformation of hydrous minerals and provide new information on the water migration process in the upper mantle of the Japan subduction zone. The present results provide the first direct evidence from seismic anisotropy for dehydration embrittlement of hydrous minerals as the main contributor to generate intermediate-depth earthquakes in subducting slabs.

Methods

We used P-wave radial anisotropy tomography27, 54 to invert for the 3-D anisotropic structure down to a depth of 380 km beneath the Japan Islands. Rocks in the Earth are generally assumed to have the anisotropy with hexagonal symmetry and the P-wave slowness can be written as55

$$S={S}_{0}(1+{M}_{1}\,\cos (2\theta )),$$
(1)

where S is the total slowness, S 0 is the average slowness (i.e., isotropic component), M 1 is the parameter for anisotropy, and θ is the angle between the propagation vector and the symmetry axis. Thus, for weak radial anisotropy with a vertical hexagonal symmetry axis, equation (1) could be rewritten as27, 54

$$S={S}_{0}(1+{M}_{1}\,\cos (2i)),$$
(2)

where i is the incident angle of a ray path. Following equation (2), the total travel time T of the ray can be expressed as

$$\begin{array}{rcl}T & = & \sum _{k}{T}_{k},\\ {T}_{k} & = & dS=d(1+{M}_{1k}\,\cos (2{i}_{k}))/{V}_{k},\end{array}$$
(3)

where T k is the travel time of the kth ray segment with a length d, V k is the isotropic velocity at the middle point of the kth ray segment, M 1k is the parameter for the radial anisotropy at the middle point of the kth ray segment, and i k is the incident angle of the kth ray segment. As a result, the partial derivatives to the velocity and radial anisotropy are expressed as

$$\begin{array}{rcl}\frac{\partial T}{\partial {V}_{k}} & = & -d(1+{M}_{1k}\,\cos (2{i}_{k}))/{V}_{k}^{2},\\ \frac{\partial T}{\partial {M}_{1k}} & = & \frac{d}{{V}_{k}}\,\cos (2{i}_{k}).\end{array}$$
(4)

Similar to the isotropic tomography method56, the parameters V k and M 1k are defined at (φ, λ, h), and f(φ, λ, h) is calculated by using a linear interpolation function

$$f(\varphi ,\lambda ,h)=\sum _{i=1}^{2}\sum _{j=1}^{2}\sum _{k=1}^{2}\,f({\varphi }_{i},{\lambda }_{j},{h}_{k})[(1-|\frac{\varphi -{\varphi }_{i}}{{\varphi }_{2}-{\varphi }_{1}}|)(1-|\frac{\lambda -{\lambda }_{j}}{{\lambda }_{2}-{\lambda }_{1}}|)(1-|\frac{h-{h}_{k}}{{h}_{2}-{h}_{1}}|),$$
(5)

where φ is latitude, λ is longitude, and h is the depth from the Earth’s surface, φ i , λ j and h k represent the coordinates of the eight grid nodes surrounding the point (φ, λ, h). Thus, a travel-time residual dT for the nth local earthquake (or a relative travel-time residual dT for the nth teleseismic event) to the mth station can be expressed as

$$\begin{array}{rcl}dT & = & {(\frac{\partial T}{\partial \varphi })}_{mn}{\rm{\Delta }}{\varphi }_{n}+{(\frac{\partial T}{\partial \lambda })}_{mn}{\rm{\Delta }}{\lambda }_{n}+{(\frac{\partial T}{\partial h})}_{mn}{\rm{\Delta }}{h}_{n}+{\rm{\Delta }}{T}_{0n}\\ & & +\sum _{p}(\frac{\partial T}{\partial {V}_{p}}{\rm{\Delta }}{V}_{p})+\sum _{q}(\frac{\partial T}{\partial {M}_{1q}}{\rm{\Delta }}{M}_{1q})+{E}_{mn},\end{array}$$
(6)

where φ n , λ n , h n and T 0n are the latitude, longitude, focal depth and origin time of the nth event, respectively; the symbol Δ denotes perturbation of a parameter; V p is the isotropic velocity at the pth node of the 3-D grid net for the isotropic velocity structure; and M 1q is the radial anisotropy parameter at the qth node of the 3-D grid net for anisotropy. The grid intervals of the two grid nets can be the same or different. E mn represents higher-order terms of perturbations and errors in the data. The first four terms on the right side of equation (6) are contributions of the hypocentral parameters, which can be calculated using the method57. The equation (6) is solved using the LSQR algorithm58 with damping and smoothing regularizations. The amplitude β of the radial anisotropy is expressed as

$$\beta =\frac{{V}_{ph}-{V}_{pv}}{2{V}_{0}}=\frac{{M}_{1}}{1-{M}_{1}^{2}},$$
(7)

where V 0 denotes the average isotropic velocity, V ph and V pv are P-wave velocities in the horizontal and vertical directions, respectively. Hence, β > 0 denotes that the horizontally propagating P-wave travels faster than the vertical one, while β < 0 denotes that the vertically propagating P-wave travels faster than the horizontal one.