## Abstract

How mantle materials flow and how intraslab fabrics align in subduction zones are two essential issues for clarifying material recycling between Earth’s interior and surface. Investigating seismic anisotropy is one of a few viable technologies that can directly answer these questions. However, the detailed anisotropic structure of subduction zones is still unclear. Under a general hexagonal symmetry anisotropy assumption, we develop a tomographic method to determine a high-resolution three-dimensional (3D) *P* wave anisotropic model of the Japan subduction zone by inverting 1,184,018 travel time data of local and teleseismic events. As a result, the 3D anisotropic structure in and around the dipping Pacific slab is firstly revealed. Our results show that slab deformation plays an important role in both mantle flow and intraslab fabric, and the widely observed trench-parallel anisotropy in the forearc is related to the intraslab deformation during the outer-rise yielding of the subducting plate.

## INTRODUCTION

The subduction of an oceanic plate beneath a continental plate causes notable structural heterogeneities and complicated dynamic processes. Although the heterogeneities can be well constrained by three-dimensional (3D) isotropic seismic velocity and attenuation models from high-resolution seismic tomography (*1*, *2*), the dynamic processes are generally indicated by seismic anisotropy (*3*–*13*), which reflects direction-dependent variations of seismic velocity, originating from lattice-preferred orientation (LPO) or shape-preferred orientation (SPO) of anisotropic materials. Seismic anisotropy can be detected by measuring shear-wave splitting (*3*, *4*, *6*–*8*), which is an integrated effect of seismic anisotropy along seismic ray paths. However, the poor depth resolution of shear-wave splitting measurements often leads to conflicting explanations even for a same feature of the observation, because it is difficult to determine where the anisotropy originates from (*3*). For example, trench-parallel fast velocity directions are widely observed at forearc seismic stations in subduction zones, but it is debated whether they come from trench-parallel mantle flow in the forearc mantle wedge (*14*) or beneath the slab (*3*), or hydrated outer-rise faults and strain in the subducted slab (*15*). The inconsistency between the observed anisotropy and theoretical fossil fabric indicates that subduction has influenced fabric alignment in the slab (*4*), but whether the anisotropy reflects the overprinted fossil fabric or simply the serpentinized outer-rise faults in the top portion of the slab (*15*) is unclear. Therefore, to fully understand subduction dynamics, it is necessary to clarify how the materials flow in the mantle and how the fabrics align in the slab. The key to answering these questions is to investigate the 3D distribution of seismic anisotropy in a subduction zone, which can directly constrain the detailed dynamic processes.

With the increasing number of seismic stations deployed and high-quality seismic data accumulated, body-wave anisotropic tomography has shown its advantages in revealing the 3D anisotropic structure of subduction zones (*11*–*13*). To simplify the calculation, the previous studies often assume an anisotropy with a hexagonal symmetry for materials in the modeling space, i.e., there exist a fast velocity plane (FVP) parallel to which a seismic wave propagates the fastest, and an FVP-normal symmetry axis along which the seismic wave propagates the slowest. This assumption is a good approximation for typical mantle minerals such as olivine (*9*, *11*, *12*). For a further simplification, the symmetry axis is assumed to be either horizontal (*11*, *13*) or vertical (*12*) for studying azimuthal anisotropy or radial anisotropy, respectively. However, the assumption of horizontal or vertical symmetry axis is not reasonable at least in two cases that are common in subduction zones: (i) complex mantle flow with a flow direction neither parallel nor perpendicular to any of the modeling axes, and (ii) intraslab anisotropy affected by the geometric inclining of a subducting slab. Therefore, the horizontal or vertical symmetry axis assumption should be released for better revealing the intraslab fabric alignment and the 3D mantle flow field in the subduction zones.

## RESULTS

In this work, we develop a method of *P* wave anisotropic tomography to reveal the detailed 3D heterogeneous and anisotropic structure of the Japan subduction zone, under the hexagonal anisotropic assumption with a freely oriented symmetry axis. Our method is modified from the effective grid-discontinuity tomographic method for isotropic medium (*1*). The surface geometry of the subducting Pacific slab (Fig. 1) is considered in modeling the subduction zone. We adopt three anisotropic parameters, in addition to an isotropic velocity parameter for heterogeneity, to describe the 3D seismic anisotropic property of rock medium. To accurately calculate theoretical travel times and reduce the trade-off between the isotropic velocity and anisotropy (*16*), we approximate the theoretical travel times using the second-order Taylor series with respect to the model parameters and use the nonlinear approximate model to fit the observed travel time data. This makes the tomographic problem a nonlinear inversion, which is apparently different from the previous studies of anisotropic tomography (*11*–*13*, *17*). Because of a large amount of computational cost caused by the second-order terms, we adopt a highly efficient parallel computation scheme, which involves 192 central processing units (CPUs) simultaneously, to accelerate the calculation. More details about the tomographic method can be found in Materials and Methods.

Our *P* wave travel time data were recorded at 1952 Kiban seismic stations deployed on the Japan islands (Fig. 1). The dataset contains a total of 692,248 arrival times of 7098 local earthquakes (Fig. 2A) and 491,770 relative travel time residuals measured from vertical component seismograms of 747 teleseismic events (Fig. 2B). The local earthquakes are precisely located by the dense seismic network (Fig. 1). The teleseismic events (*M* ≥ 6.0) used are located in a distance range of ~30° to 90° from central Japan (Fig. 2B). Many of the local events occurred in the subducting Pacific slab close to the middle portion of the slab (>25-km depth below the slab surface; Fig. 2A), which can notably improve the tomographic resolution in the Pacific slab. As a result, under a slow symmetry axis assumption (see Materials and Methods), different patterns of anisotropy and isotropic *P* wave velocity (Vp) are revealed in the crust, mantle wedge, subducting slab, and subslab mantle (Figs. 3 to 5). Our results reveal trench-normal and upright FVPs with an anisotropic amplitude of <5% in the mantle wedge (Figs. 3B and 4, A and B), trench-parallel FVPs normal to the slab surface (amplitude 3 to 8%) in the slab top portion (<10 km below the slab surface; Figs. 3, D to F; 4, D to F; and 5A), complex FVPs (amplitude 1 to 5%) in the slab middle portion (>20 km below the slab surface; Figs. 3, D to F; 4, D to F; and 5B), and subhorizontal FVPs (amplitude <8%) in the subslab mantle where the slab bends laterally (at depths of 150 to 300 km; Figs. 3C and 4C). Figure 6 shows the estimated mantle flows from the obtained Vp anisotropic tomography (Figs. 3 and 4).

We also conducted a tomographic inversion under a fast symmetry axis assumption to further determine the distribution of three olivine axes (i.e., [100], [010], and [001]) in the mantle wedge, subducting slab, and subslab mantle (see Materials and Methods for details). In the top portion of the slab, the olivine [100] axis is generally trench parallel, whereas in the mantle wedge, the [100] axis is generally vertical (Fig. 7A). In the subslab mantle, the olivine [100] axis is generally parallel to the slab in the obtained subhorizontal FVPs (Fig. 7B).

## DISCUSSION

Trench-parallel mantle flows are suggested by previous studies to interpret the measured trench-parallel seismic anisotropy from *S* wave splitting observations (*3*, *14*, *18*), which can hardly be explained by 2D mantle flows in a trench-normal plane driven by the oceanic plate subduction. The inferred models of those studies are based on an assumption that the trench-parallel anisotropy originates from the mantle above or below the slab (*3*, *14*, *18*). Our tomographic results show that the FVPs revealed in the mantle wedge align in a consistent way: They are generally upright (i.e., normal to the horizontal plane) and perpendicular to the trench (at depths of 30 to 150 km; Figs. 3B; 4, A and B; and 6C), indicating vertical and trench-normal mantle flows in the planes perpendicular to the trench (Fig. 6B). The upright and trench-normal FVPs are highly consistent with the existence of low-velocity anomalies in the mantle wedge extending from depth beneath the back-arc to the roots of active arc volcanoes (30- to 150-km depths; Figs. 3, D to F, and 4, D to F), which unquestionably reflect upwelling hot materials driven by corner flows in the mantle wedge (*1*, *12*, *13*). Considering that a vertical olivine [100] axis is generally revealed by the fast symmetry axis inversion, it can be inferred that upwelling flows dominate the mantle wedge dynamics (Fig. 7A) and the slip planes are parallel with the obtained FVPs if A-type olivine is assumed in the wedge. Similar upright FVPs are also revealed in the subslab mantle (at ~300-km depth and in the latitude range of 38° to 41°; Fig. 4C), probably reflecting entrained flow beneath the subducting Pacific slab. These results suggest that classical 2D slab-driven mantle flows dominate mantle dynamics in the wedge and beneath the slab (Fig. 6B). Numerical (*19*) and laboratory experimental (*20*) studies suggest that complex return flows occur around a slab edge, which interchange the materials above and below the subducting slab. However, our tomographic results do not show such a feature, probably because the Pacific slab is continuous without edge under Japan and the slab is old, thick, relatively intact, and has subducted down to the mantle transition zone (*21*), providing no channels or windows for material flows between the mantle wedge and the subslab mantle.

In addition to the upright FVPs, oblique and flat FVPs are also observed in the mantle wedge (at 40- to 70-km depths; Figs. 3B and 4A) and subslab mantle (at 200- to 300-km depths; Figs. 3C and 4C), which are located close to the outer and inner corners where the slab bends in the lateral direction. It has been suggested that lateral deformation of the slab can cause trench-parallel anisotropy in the inner side of the lateral curve close to the slab due to strong trench-parallel stretching (*22*). The flat FVPs in the subslab mantle (Fig. 3C), indicating lateral mantle flow or stretching (Fig. 6B), are possibly caused by this mechanism. The trench-parallel olivine [100] axes revealed by the fast symmetry axis inversion also support the mechanism of trench-parallel stretching (Fig. 7B). The oblique FVPs in the mantle wedge may be caused by very local perturbations of mantle flow associated with variations of the slab geometry (Fig. 6B).

Slab deformation can influence not only the mantle flow but also the fabric alignment in the slab itself. To demonstrate this, we rotate FVPs in the middle portion of the slab as “restored FVPs,” based on the morphology of the slab surface (the rotation method is explained in Fig. 5C, and the obtained restored FVPs are shown in Fig. 5B). Seismic anisotropy in the slab is considered to reflect a fossil fabric forming at the mid-ocean ridge and originating at the time of plate formation (*12*, *23*). Such a fabric within the slab may be also overprinted because of the slab deformation during subduction, as inferred by a recent study of the Nazca slab (*4*). For the Pacific slab beneath Japan, if the fossil fabric is the only factor causing the anisotropy in the slab, then the FVPs should be simply distributed with a consistent orientation parallel to the slab surface (i.e., flat restored FVPs) (*23*). However, our results do not show such a feature (Fig. 5B). To further clarify whether the slab deformation has influenced the fabric alignment in the slab, we compare the distribution of the slab curvature (Fig. 5F) with that of the included angle between the symmetry axis of the FVPs and the outward normal direction of the slab surface (Fig. 5E). The areas with high included angles (Fig. 5E), where the FVPs are approximately perpendicular to the slab surface, exhibit a pattern very similar to that of high-curvature areas (Fig. 5F). We think that slab surface normal shear strain associated with the slab bending is likely responsible for the slab surface normal LPO of the fabric in the slab. For example, numerical modeling studies [e.g., (*24*)] show that outer-rise bending can cause faulting and strain in the very deep part of the slab. Because the outer-rise bending–related anisotropic amplitudes in the top portion of the slab can reach 3 to 8% (as discussed below), which are obviously greater than those caused by corner flows in the mantle wedge (<5%), we deem that the anisotropy caused by the slab bending can overprint the fossil anisotropy formed at the mid-ocean ridge.

Our tomographic results reveal no systematic trench-parallel anisotropy in the mantle beneath the forearc (Fig. 4A) and the middle portion of the slab (Fig. 5B). Hence, the possible source of trench-parallel anisotropy from shear-wave splitting measurements at the forearc stations, which has attracted wide attention and has been debated (*3*, *14*, *15*, *18*, *22*, *25*, *26*), can only exist in the crust of the overriding Okhotsk plate or the top portion of the subducting Pacific slab. The anisotropy in the crust can be caused by multiple geological features, such as lithological variations, active faults, and aligned fractures formed by long-term shear strain (*11*). Because these geological features can be very different from one area to another, the consequent anisotropy in the crust may exhibit an irregular distribution, which is well confirmed by the disordered FVPs shown in our result (Fig. 3A). Although the anisotropic amplitudes of the crustal FVPs can reach 10%, they are not likely to form the systematic trench-parallel anisotropy because of their irregular distribution. Another potential source of the trench-parallel anisotropy is in the top portion of the slab, which exhibits a very clear feature of trench-parallel and slab surface normal FVPs (see the images at 3-km depth below the slab surface in Fig. 5, A and D). In the top portion of the slab, hydrated faults filled with highly anisotropic hydrous minerals are predicted to form due to the strain in the slab when the plate passes through the outer-rise yielding zone (*15*, *27*, *28*). These faults are thought to be trench parallel and steeply dipping (*15*). The trench-parallel FVPs revealed by our result may reflect effects of the hydrated faults in the slab. Furthermore, the olivine [100] axes revealed by our fast symmetry axis inversion are also trench parallel (Fig. 7A), indicating that the trench-parallel direction is the fastest direction of anisotropy in the top portion of the slab. Although the anisotropic layer containing the hydrated faults may be thin (<10 km), the integrated anisotropy sampling the layer can be large because alignment of the filling minerals can result in a very strong anisotropy [e.g., ~32% for serpentine (*25*) and >10% for amphibole (*15*)]. Our results show that the anisotropic amplitudes (3 to 8%) in the top portion of the slab are generally greater than those in the mantle wedge (<5%) and in the middle portion of the slab (1 to 5%). Considering that the anisotropic amplitudes could have been somewhat suppressed by damping and smoothing regularizations adopted in the tomographic inversion (see Materials and Methods), the actual anisotropic amplitudes are possibly even greater. If we assume that a seismic wave propagates the slowest in the direction perpendicular to the fault planes, then the slab surface normal FVPs indicate that the hydrated faults are also perpendicular to the slab surface.

The anisotropic tomography method developed in this study requires a greater amount of arrival-time data than those for the isotropic tomographic inversions, because the orientation of an FVP is determined by a combination of three parameters (see Materials and Methods). Not only the amount of arrival times but also the crisscrossing of seismic rays used are of crucial importance for obtaining a robust 3D model of anisotropic tomography. Fortunately, with the increasing number and density of seismic stations deployed around the world in the past decade, this requirement is met in more and more regions.

## MATERIALS AND METHODS

### Model parameterization

We determine 3D isotropic Vp and anisotropy simultaneously by jointly inverting travel time data of local and teleseismic events. Our method is modified from the tomographic inversion code Tomog3d for studying the 3D isotropic velocity structure (*1*, *29*). To take into account the lateral depth variations of the surface of the subducting Pacific plate (Fig. 1), we partition the modeling space into two parts: one part is the subducting Pacific slab, and the other part is the region surrounding the Pacific slab (*30*). The geometry of the slab surface (Fig. 1) has been well constrained by many previous studies (for details, see (*30*) and (*31*)], but the lower boundary of the Pacific slab beneath Japan has not been well determined yet. In this study, we assume a consistent slab thickness of 85 km as estimated by teleseismic tomography (*30*). Two 3D grids with the same lateral grid interval of 0.33° are set up for expressing the 3D isotropic Vp structures of the two parts of the study volume. Thirteen meshes of the first 3D grid are set at depths of 10, 25, 40, 70, 100, 150, 200, 300, 400, 500, 600, 700, and 900 km in the crust and mantle outside the Pacific slab, whereas three meshes of the second 3D grid are set within the Pacific slab at depths of 5, 25, and 60 km beneath the slab surface. Similar to the grid settings for the isotropic Vp structure, additional two 3D grids with a lateral grid interval of 0.5° are also set up to parameterize the seismic anisotropy structure outside and inside the Pacific slab. Meshes of the anisotropic grids are set at depths of 10, 40, 70, 100, 150, 200, 300, 400, and 700 km in the crust and mantle and within the Pacific slab at 5- and 40-km depths below the slab surface. The grid settings for anisotropy in the vertical direction are carefully determined by a series of resolution tests. Because very few intermediate-depth and deep earthquakes occur in the lower part of the slab, the anisotropy there is hard to be constrained if we set a grid mesh there (e.g., at 60-km depth below the slab surface). Therefore, we only arrange two grid meshes at 5- and 40-km depths below the surface of the Pacific slab to represent anisotropy in the slab upper portion and the anisotropy in the slab middle-lower portion, respectively. Hence, the obtained anisotropy at 40-km depth below the slab surface actually represents a composite effect of the anisotropy in the middle and lower parts of the slab. We adopt the CRUST1.0 model (http://igppweb.ucsd.edu/~gabi/rem.html) as the starting model for the crust. For the mantle part, we start with the IASP91 Earth model (*32*) for our tomographic inversion. Following previous studies (*1*, *13*, *29*, *30*), in the starting model, we add +4% isotropic Vp anomaly to the grid nodes in the Pacific slab so that the slab structure can be better constrained.

### Anisotropy with a tilting hexagonal symmetry axis

In a medium with anisotropy having a tilting hexagonal symmetry axis (HSA), *P* wave anisotropic slowness can be expressed as (*11*–*13*)*S*_{0} is the average *P* wave slowness (i.e., the isotropic component), *M* is a positive parameter quantifying the anisotropic amplitude, and θ is the angle between the HSA and the ray propagation direction. A unit vector along the ray propagation direction can be expressed as ** l**(sin

*p*sin λ, sin

*p*cos λ, cos

*p*) (fig. S1A), whereas a unit vector along the HSA can be expressed as

**( sin**

*n**q*sin γ, sin

*q*cos γ, cos

*q*)(fig. S1B). λ and γ are azimuths of the two vectors, whereas

*p*and

*q*are incident angles of the two vectors. To ensure the HSA orientation to be unique, we assume 0 ≤

*q*< π and 0 ≤ γ < π. Thus, cosine of the two vectors (fig. S1C) can be written as

*P*wave slowness can be rewritten as

*d*= sin λ sin

*p*,

*e*= cos λ sin

*p*,

*f*= cos

*p*,

*A*,

*B*, and

*C*are three anisotropic parameters to be determined by tomographic inversion. According to the definitions of

*q*and γ,

*B*and

*C*can be either negative or positive, whereas

*A*should be zero or positive. Once the three anisotropic parameters are determined, the orientations of FVP and HSA can be easily obtained (fig. S1, B and D). In the following, we adopt a T-shaped symbol to represent the FVP, whose long axis denotes the intersecting line between the FVP and the paper plane (i.e., the fast velocity direction in the paper plane) and whose short bar normal to the long axis is determined by the angle between the FVP and the paper profile, as shown in Figs. 3B and 4B. To further clarify how the 3D FVPs look like in a horizontal plane and a trench-normal vertical profile, we show three typical FVPs in fig. S2 as an illustration.

### Nonlinear observation equations and tomographic inversion

Slowness (*S*_{0}) at the isotropic grid nodes and anisotropic parameters (i.e., *A*, *B*, and *C*) at the anisotropic grid nodes are determined by tomographic inversion. Previous studies of anisotropic tomography linearize theoretical travel times by the first-order Taylor expansion referencing to a starting model and then iterate the calculation and generate a model sequence to approach the optimal model (*11*–*13*, *17*, *33*). These linearization methods are effective when (i) the starting model is close enough to the optimal model or (ii) the nonlinearity of the tomographic problem is weak, so that a linearization using the first-order Taylor series is accurate enough to approximate the theoretical model. Unfortunately, for the anisotropy with a freely oriented HSA (Eq. 3), neither of the two preconditions is satisfied, because we do not have a good starting model for the 3D anisotropic structure, and the anisotropic model is highly nonlinear due to the complex relations of the anisotropic parameters (i.e., *A*, *B*, and *C* in Eq. 3). In this case, the conventional linear methods may easily fail to converge to the optimal solution, unless initial models and parameters are chosen very carefully for the inversion (*17*).

To address this problem, in the present study, we adopt the second-order Taylor series to better estimate theoretical travel times that are used to determine the slowness and anisotropy parameters simultaneously by fitting the observed data in the least-squares sense. The theoretical travel time can be written as*j*th local earthquake to the *i*th station for the starting model, and *m _{n}* represents the

*n*th unknown in the variable list of the parameterized anisotropic model, i.e.,

*S*

_{0}at the isotropic grid nodes or

*A*,

*B*, and

*C*at the anisotropic grid nodes. To take advantage of the information of the second-order derivatives in Eq. 4, we adopt a quasi-Newton method, which is called the limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm for bound constrained optimization (L-BFGS-B) (

*34*,

*35*), to find an optimal solution for best explaining the observed data in the least-squares sense. The L-BFGS-B is a limited-memory algorithm for solving a large nonlinear optimization problem subject to simple bounds on the variables and is adapted to problems in which information on the Hessian matrix is difficult to extract or the matrix is dense (

*34*). These characteristics of the algorithm make it particularly suitable for solving the nonlinear fitting problem of this study, because (i) our tomographic problem is very large (~1.2 million nonlinear equations); (ii) the large number of second-order terms make the Hessian matrix very dense; and (iii) nonnegative bounds have to be added to the anisotropic parameter

*A*.

Because we use the second-order derivatives of travel times, a great number of second-order terms have to be calculated and stored. As a result, the ~1.2 million nonlinear equations lead to a total of ~0.38 billion first-order terms and ~10 billion second-order terms, which require a great computational capacity and storage memory for the computer. To accelerate the calculation, we conduct the anisotropic tomography by adopting a high-efficiency parallel computation scheme, in which the computation is almost equally partitioned for each CPU. Our final inversion involves 192 CPUs simultaneously and costs only 3.2 hours for a single tomographic iteration. The complete inversion procedure of this study is shown in fig. S3.

Damping and smoothing regularizations are used to obtain an optimal solution by regularizing the tomographic problem. We first determine an optimal combination of smoothing parameters for the isotropic Vp and anisotropy (fig. S4A) without considering damping regularization, and then determine an optimal damping parameter using the trade-off curve that is obtained by performing a number of trial inversions (fig. S4B). To take account of ray path changes in the obtained 3D anisotropic model, we conduct several iterations, and in each iteration, we recalculate the ray paths with the new 3D velocity model obtained by the previous iteration. It takes four iterations for the inversion procedure to reach an optimal tomographic model. After four iterations, the root mean square (RMS) travel time residual decreases to its minimum and keeps stable, and the distribution of travel time residuals is notably improved. The final RMS travel time residual is reduced to 0.317 from 0.534 s before the inversion, and the variance reduction is 64.8%.

### Resolution analysis

We performed extensive checkerboard resolution tests (CRTs) with synthetic input models to verify our anisotropic tomography method and to evaluate the data capacity for resolving the 3D underground structure beneath the seismic network. To construct an input checkerboard model, isotropic *P* wave slowness perturbations of ±6% (i.e., +6.4 and −5.7% isotropic Vp perturbations) are assigned to adjacent nodes of the 3D isotropic grid. In addition, anisotropic anomalies with the same amplitude of 5% but different HSA orientations are assigned to adjacent nodes of the 3D anisotropic grid. Synthetic travel time data are generated for the input checkerboard model. Gaussian noise (−0.2 to +0.2 s) with an SD of 0.1 s is added to the synthetic data. A series of checkerboard models with lateral grid intervals of 0.33° to 1.0° are tested. The test results show that the resolution for the isotropic Vp part is ~0.33° in the crust, mantle wedge, slab, and subslab mantle, whereas the resolution of the anisotropic part is ~0.5° in the crust, mantle wedge, and slab and ~1.0° in the subslab mantle. Figures S5 and S6 show the CRT results for the crust and upper mantle with a lateral grid interval of 0.33° for the isotropic Vp and 0.5° for the anisotropy. Figure S7 shows the same CRT results for the slab part, which indicates that the ray coverage is good enough to resolve the isotropic Vp and anisotropy structures of the Pacific slab.

To confirm the main features of the obtained isotropic Vp model and anisotropy, we conducted a restoring resolution test (RRT) whose procedure is similar to that of the CRT, but the obtained 3D tomographic model (Figs. 3 and 4) is adopted as the RRT input model. Gaussian noise (−0.2 to +0.2 s) with an SD of 0.1 s is also added to the synthetic data calculated for the input model. The RRT results show that the main features of our tomographic images (both isotropic Vp and anisotropy) can be well recovered (fig. S8).

Trade-off between the obtained isotropic Vp and anisotropy should be carefully considered for a model of anisotropic tomography (*16*, *36*). We conducted two synthetic tests to evaluate the trade-off. In the first test, only the isotropic Vp part of the obtained model (Figs. 3 and 4) is retained to construct a synthetic input model, but the synthetic data are inverted for both isotropic Vp and anisotropy. The second test is similar to the first one, but only the anisotropy part of the obtained model (Figs. 3 and 4) is retained to make an input model. The output models of the two tests (figs. S9 and S10) show that the trade-off between the isotropic Vp and anisotropy is negligible and that the main features as mentioned above are not caused by the trade-off.

To evaluate the influence of the initial model on our final result, we adopted a very different initial model to conduct a tomographic inversion of our dataset. In this case, we set *A* = 0.1731, *B* = 0, and *C* = 0 at each node of the anisotropic grid. As a result, the FVPs in the initial model are all upright with a north-south fast direction in the horizontal plane and with an anisotropic amplitude of 3%. The obtained result of this inversion is displayed in figs. S11 and S12, showing main features very similar to those obtained with the isotropic 1D initial model (Figs. 3 and 4). The most remarkable difference of the two results is in the anisotropy of the subslab mantle, where slab-parallel FVPs with an amplitude of ~3% appear in the latitude range of 37° to 41° at 200-km depth (fig. S11C). While this feature may reflect flows entrained by the subducting slab (*3*), we tend to not explain it because of its dependence on the initial model.

### Comparison with previous models

If we retain the azimuthal part of our anisotropic model (Figs. 3 and 4), i.e., restraining the FVPs upright with respect to the horizontal plane by setting parameter *C* = 0 as previous studies did, then the obtained model is generally consistent with the previous results [e.g., (*13*, *33*, *37*, *38*)], for example, the trench-normal fast velocity directions in the mantle wedge beneath the volcanic front and the back-arc and a general trench-parallel fast velocity direction in the slab. However, our model reveals more information because the symmetric axis is allowed to tilt, especially in the dipping Pacific slab. Because the previous studies (*13*, *33*, *37*, *38*) strictly limited the symmetric axis to be either horizontal or vertical, they could not reveal some important anisotropic features such as the slab surface normal FVPs in the upper portion of the slab as revealed by this study (Figs. 3, 4, and 6E).

### Determination of three olivine axes

Note that seismic waves travel at different speeds along the three orthogonal axes of olivine (*9*, *11*, *33*, *39*). Equations 1 to 3 are derived on the basis of a slow symmetry axis assumption, which means that seismic waves travel slowest along the direction of the symmetry axis, i.e., the olivine [010] axis along the symmetry axis and the olivine [100] and [001] axes in the FVP. Another alternative is the fast symmetry axis assumption, which identifies the fastest direction ([100]) as the symmetry axis. Under this assumption, the olivine [010] and [001] axes are laid out in the slow velocity plane that is orthogonal to the symmetry axis. We adopt the slow symmetry axis assumption because this parameterization is more suitable for materials with aligned cracks or faults (*9*, *11*), which play a crucial role in generating SPO-related anisotropy in the subducting slab. To better determine the three axes of olivine, we also conducted an inversion under the fast symmetry axis assumption. Thus, the olivine [100] and [010] axes can be determined from the fast and slow symmetry axes obtained by the two inversions under the two different assumptions (fig. S13). The olivine [010] axis is generally orthogonal to the [100] axis (fig. S13), although they are determined by separate inversions.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/4/eabc9620/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank the data centers of the Japanese seismic networks and the JMA Unified Earthquake Catalog (www.hinet.bosai.go.jp) for providing the high-quality waveform and arrival-time data used in this study. Some of the arrival-time data were measured from the original seismograms by the staffs of the Research Center for Prediction of Earthquakes and Volcanic Eruptions, Tohoku University. We are very grateful to W. Wei and G. Toyokuni for the helpful discussions and to X. Liu for collecting part of the data. M. Ritzwoller (the editor) and two anonymous referees provided thoughtful review comments and suggestions, which have improved this paper.

**Funding:**This work was supported by a research grant (no. 19H01996) from the Japan Society for the Promotion of Science to D.Z. and a grant (no. 41904045) from the National Natural Science Foundation of China to Z.W.

**Author contributions:**Z.W. and D.Z. wrote the manuscript and prepared the figures, and both authors reviewed the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2021 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).