Seismic Anisotropy Tomography and Mantle Dynamics

Abstract

Seismic anisotropy tomography is the updated geophysical imaging technology that can reveal 3-D variations of both structural heterogeneity and seismic anisotropy, providing unique constraints on geodynamic processes in the Earth’s crust and mantle. Here we introduce recent advances in the theory and application of seismic anisotropy tomography, thanks to abundant and high-quality data sets recorded by dense seismic networks deployed in many regions in the past decades. Applications of the novel techniques led to new discoveries in the 3-D structure and dynamics of subduction zones and continental regions. The most significant findings are constraints on seismic anisotropy in the subducting slabs. Fast-velocity directions (FVDs) of azimuthal anisotropy in the slabs are generally trench-parallel, reflecting fossil lattice-preferred orientation of aligned anisotropic minerals and/or shape-preferred orientation due to transform faults produced at the mid-ocean ridge and intraslab hydrated faults formed at the outer-rise area near the oceanic trench. The slab deformation may play an important role in both mantle flow and intraslab fabric. Trench-parallel anisotropy in the forearc has been widely observed by shear-wave splitting measurements, which may result, at least partly, from the intraslab deformation due to outer-rise yielding of the incoming oceanic plate. In the mantle wedge beneath the volcanic front and back-arc areas, FVDs are trench-normal, reflecting subduction-driven corner flows. Trench-normal FVDs are also revealed in the subslab mantle, which may reflect asthenospheric shear deformation caused by the overlying slab subduction. Toroidal mantle flow is observed in and around a slab edge or slab window. Significant azimuthal and radial anisotropies occur in the big mantle wedge beneath East Asia, reflecting hot and wet upwelling flows as well as horizontal flows associated with deep subduction of the western Pacific plate and its stagnation in the mantle transition zone. The geodynamic processes in the big mantle wedge have caused craton destruction, back-arc spreading, and intraplate seismic and volcanic activities. Ductile flow in the middle-lower crust is clearly revealed as prominent seismic anisotropy beneath the Tibetan Plateau, which affects the generation of large crustal earthquakes and mountain buildings.

Article Highlights

• Advances in the theory and application of seismic anisotropy tomography are reviewed. The results provide unique constraints on the 3-D seismic structure and geodynamic processes

• Slab deformation plays an important role in both mantle flow and intraslab fabric. Toroidal mantle flow is widely observed in and around a slab edge or slab hole

• Trench-normal fast anisotropy generally occurs in the mantle wedge and subslab mantle, reflecting subduction-driven corner flow and asthenospheric shear deformation, respectively

• Significant anisotropy occurs in the big mantle wedge (BMW) under East Asia, reflecting convection associated with deep subduction of the Pacific plate and its stagnation in the mantle transition zone. Geodynamic processes in the BMW caused craton destruction, back-arc spreading, and intraplate volcanic and seismic activities

1 Introduction

Our knowledge of the Earth’s interior structure mainly comes from observations and analyses of seismic waves from natural earthquakes and artificial explosions. Seismic imaging is the most powerful tool to investigate the three-dimensional (3-D) structure of the Earth’s interior. A great number of studies have been made using this approach, which have provided us with new insights into the depth extension of surface geological features and the control of mantle dynamics on lithospheric processes such as earthquakes, volcanic eruptions and mountain building (e.g., Anderson 1989; Zhao 2015). Seismologists mainly adopt three kinds of physical parameters to express the structure of the Earth’s interior, i.e., seismic velocity, attenuation, and anisotropy. Seismic velocity is P or S wave propagating speed (Vp or Vs) and can be obtained by using body-wave travel times or surface-wave dispersions measured from seismograms. The loss of energy experienced by seismic waves as they travel can be described in terms of seismic attenuation, which is affected by temperature, lithology, melt and fluid content of rocks through which seismic waves propagate. Seismic anisotropy is a phenomenon that some properties of seismic waves change with the wave propagating direction. Both seismic observations and laboratory experiments show that the solid Earth contains strongly anisotropic materials and its structure and dynamics can be better understood by studying seismic anisotropy, because anisotropy contains important information on past and present deformations in the lithosphere, asthenosphere, and subducting slabs. Adopting a tomographic approach to analyze seismological data (e.g., travel times, amplitudes and waveforms), we are able to determine 3-D distributions of the three physical parameters in the Earth’s interior, which are called seismic velocity tomography, attenuation tomography, and anisotropic tomography (e.g., Zhao 2015).

To date, a majority of tomographic studies are concerned about isotropic velocity tomography, i.e., the Earth structure is assumed to be isotropic for the propagation of seismic waves. This is because in many cases seismic data sets used in tomography do not have enough data and/or ray-path coverage to resolve the anisotropic part of the signal. In addition, this isotropic assumption is made for mathematical convenience. Although most of seismic data can be modeled well with this isotropic assumption, it does not mean that the Earth has an isotropic structure. Seismic anisotropy can cause the biggest variations in seismic velocity, which can be even larger than those caused by changes in temperature, composition or mineralogy. Therefore, seismic anisotropy is a first-order effect to the seismological data (Anderson 1989).

Now it is widely accepted that seismic anisotropy exists almost everywhere in the solid parts of the Earth, i.e., the crust, mantle, and the inner core (see Silver 1996; Fouch and Rondenay 2006; Plomerova et al. 2012; Long 2013 for detailed reviews). Seismic anisotropy is mainly caused by lattice-preferred orientation (LPO) and shape-preferred orientation (SPO) of the rocks in the Earth (e.g., Babuska and Cara 1991; Karato et al. 2008). Both body-wave and surface-wave data are exploited to investigate seismic anisotropy. The body-wave methods include shear-wave splitting (SWS) measurements, receiver functions, and travel-time tomography, which have higher spatial resolution than that of the surface-wave methods (e.g., Fouch and Rondenay 2006; Plomerova et al. 2012; Zhao 2015). SWS, also called seismic birefringence, is a feature that a shear wave splits into two polarized shear waves exhibiting different speeds when propagating through an anisotropic medium. Two splitting parameters (φ, δt) can be measured from horizontal-component seismograms, which are the polarization direction of the fast quasi-shear phase (φ) and the time difference (δt, also called delay time) between the fast and slow shear waves (e.g., Long and Silver 2009; Huang et al. 2011a, 2011b). To date, measuring SWS has been the most popular method for detecting seismic anisotropy, partly because it is quite simple and straightforward and can be made for a single station. By analyzing short-period local shear waves or long-period waves passing through the core (SKS, SKKS, etc.), the observed SWS is usually regarded as reflecting seismic anisotropy at a depth or in a depth range of the crust or mantle beneath a seismic station. Thus, the two SWS parameters (φ, δt) can be adopted to reveal azimuthal anisotropy in the crust and/or the mantle. So far, many researchers have made SWS measurements to investigate seismic anisotropy in the crust and/or mantle beneath different parts of the world (see detailed reviews by Long 2013, 2016).

Although measuring SWS is popular and effective to detect seismic anisotropy, this method has poor depth resolution, that is, the depth or depth range of an anisotropic source under a seismic station cannot be determined. Hence, the SWS results are usually hard to be interpreted uniquely. This shortcoming of the SWS measurements can be overcome by conducting seismic anisotropic tomography. To date, many researchers have used Pn wave travel-time data to study P-wave azimuthal anisotropy (e.g., Hearn 1996; Mutlu and Karabulut 2011; Du et al. 2022). Because Pn waves refract at the Moho discontinuity and travel in the uppermost mantle, Pn tomography determines two-dimensional (2-D) Vp variations and azimuthal anisotropy in the uppermost mantle directly beneath the Moho. In contrast, P-wave anisotropic tomography (PWAT) can map the 3-D distribution of Vp azimuthal and radial anisotropies in the crust and upper mantle beneath a seismic network (e.g., Wang and Zhao 2008, 2013; Liu and Zhao 2017; Ishise et al. 2018). PWAT is a relatively new technique and has been mainly developed and applied in the past one and a half decades (see Zhao 2015 for a detailed review). PWAT has a potential problem, i.e., there could be trade-off between Vp heterogeneity and anisotropy. However, recent studies have shown that the trade-off issue can be resolved when the study volume is densely and uniformly covered by ray paths used in a tomographic inversion (e.g., Huang et al. 2015; Liu and Zhao 2017; Gou et al. 2019a; Zhao 2021) and if the theoretical travel time is calculated more precisely (Wang and Zhao 2021). To improve the depth resolution of SWS measurements, Mondal and Long (2019, 2020) conducted probabilistic finite-frequency SKS splitting intensity tomography, which provides both lateral and depth constraints on anisotropic structure of the upper mantle. Their technique is based on a Markov chain Monte Carlo approach to searching parameter space, and they used finite-frequency sensitivity kernels to relate model perturbations to splitting intensity observations.

Twenty-one independent elastic moduli are required to fully express a general anisotropic medium, which are hard to handle, particularly in the case of tomographic imaging, because the available seismological data are far from complete and ideal. Hence, we have to make some approximations to reduce the number of elastic moduli for the anisotropic medium. Hexagonal symmetry is such an approximation to the Earth material, which can greatly reduce the number of physical parameters describing seismic anisotropy (e.g., Christensen 1984; Park and Yu 1993; Maupin and Park 2007). In practical studies, researchers usually assume horizontal hexagonal symmetry when azimuthal anisotropy is concerned (Fig. 1a) in SWS analysis (e.g., Crampin 1984; Silver 1996; Savage 1999; Long 2013), Vp studies (e.g., Hess 1964; Backus 1965; Raitt et al. 1969; Hearn 1996; Eberhart-Phillips and Henderson 2004; Wang and Zhao 2008, 2013), and surface-wave studies (e.g., Forsyth and Li 2005; Liu and Zhao 2016a). On the other hand, vertical hexagonal symmetry is assumed for radial anisotropy (Fig. 1b) in surface-wave studies (e.g., Fichtner et al. 2010; Nettles and Dziewonski 2008; Yuan et al. 2011) and body-wave Vp studies (e.g., Wang and Zhao 2013; Ishise et al. 2018; Zhao and Hua 2021). From a computational point of view, it is necessary to assume hexagonal anisotropy with a horizontal or vertical symmetry axis, because seismic data usually have no resolving power to constrain more complicated anisotropy.

Fig. 1

The coordinate system specifying a ray path (dashed line) and the hexagonal symmetry axis (red line) for azimuthal anisotropy (a) and radial anisotropy (b). V is propagation vector of a ray with incident angle i and azimuth $\phi$. θ is angle between the propagation vector and the symmetry axis. For azimuthal anisotropy tomography (a), the hexagonal symmetry axis is normal to the fast-velocity direction (FVD) in the horizontal plane, and ψ’ is azimuthal angle of the hexagonal symmetry axis

Laboratory experiments and studies of LPO in naturally deformed rocks (e.g., Karato 1995; Jung and Karato 2001; Kneller et al. 2005; Karato et al. 2008; Jung 2017) have shown that the geometry of anisotropy is controlled by the geometry of flow (i.e., macroscopic aspect) and by the geometry of LPO that is determined by the nature of plastic and elastic anisotropy of minerals (i.e., microscopic aspect). In a majority of the upper mantle, anisotropy is formed in the asthenosphere where flow is nearly horizontal and strain is large, and the $2\phi$ terms (see Eq. 1 in the next section) dominate. But if the flow geometry is more complex, then seismic anisotropy would become more complicated. These mineral physics studies suggest that at least some of the trench-parallel fast anisotropy can be explained assuming 2-D flow geometry by invoking the B-type fabrics. Browaeys and Chevrot (2004) showed that hexagonal symmetry accounts for ~ 75% and > 50% of the total anisotropy for olivine and enstatite, respectively. In addition, the olivine alignment along all three axes is not perfect in the mantle, and the inherent anisotropy of several minerals combines in such a way as to make the typical mantle and crustal petrofabric even closer to hexagonal. The hexagonal approximation for seismic anisotropy may capture ~ 80% of the total anisotropy in the upper mantle (Becker et al. 2006).

In the following sections, we first describe the methods of azimuthal and radial anisotropy tomography and their practical applications in the past decade to understand the 3-D mantle structure and dynamics, then introduce more recent advances of regional-scale anisotropic tomography, and finally we discuss future perspectives of this research field.

2 Azimuthal Anisotropy Tomography

Seismic travel-time tomography for 3-D isotropic velocity structure relates observed travel-time residuals (r) to perturbations of hypocentral parameters (i.e., latitude $\Delta {\varphi }_{e}$, longitude $\Delta {\lambda }_{e}$, focal depth $\Delta {h}_{e}$, and origin time $\Delta {T}_{e}$) of local earthquakes and 3-D velocity perturbations (ΔV) relative to a 1-D reference velocity model. In isotropic tomography, at a point in the study volume, P-wave is assumed to propagate with the same velocity (${Vp}_{0}$) in all directions (the corresponding slowness is defined as ${Sp}_{0}=1/{Vp}_{0}$). For weak azimuthal anisotropy (AAN), Vp and Vs variations with the azimuth of wave propagation ($\phi$) can be expressed as a linear combination of terms proportional to $\mathrm{sin}2\phi$, $\mathrm{cos}2\phi$, $\mathrm{sin}4\phi$ and $\mathrm{cos}4\phi$ (Backus 1965; Raitt et al. 1969; Crampin 1981). For SH waves, the $4\phi$ terms could be important. But, for P waves, the $4\phi$ terms are always ignored due to their relatively small effects on the observed seismic data as compared with the $2\phi$ terms. Hence, weak hexagonal anisotropy (Fig. 1a) is usually assumed to exist in the modeling space to facilitate the calculation due to the limited spatial coverage and resolution of the seismic data used. Thus, a simple expression for the velocity Vp (or slowness Sp) of a horizontally propagating P-wave (e.g., Pn wave traveling in the uppermost mantle) is as follows (Hearn 1996):

$$S_{P} = S_{P0} + M\cos 2\theta = S_{P0} + A\cos 2\phi + B\sin 2\phi$$

(1)

where θ is the angle between the propagation vector and the symmetry axis, A and B are two parameters describing AAN, and $M=\sqrt{{A}^{2}+{B}^{2}}$. The anisotropic amplitude ($\alpha$) and fast-velocity direction (FVD, $\psi$) are given by:

$$\alpha = \frac{{V_{{{\text{Pf}}}} - V_{{{\text{Ps}}}} }}{{2V_{P0} }} = \frac{M}{{S_{P0} - M^{2} /S_{P0} }}$$

(2)

$$\psi \left( {A,B} \right) = \left\{ {\begin{array}{*{20}l} {\frac{1}{2}\tan^{ - 1} \left( \frac{B}{A} \right) + \left\{ {\begin{array}{*{20}l} {\frac{\pi }{2}, A > 0} \hfill \\ {0,A < 0} \hfill \\ \end{array} } \right.} \hfill \\ { - \frac{\pi }{4},B > 0,A = 0} \hfill \\ {\frac{\pi }{4},B < 0,A = 0} \hfill \\ \end{array} } \right.$$

(3)

where $V_{{{\text{Pf}}}}$ and $V_{{{\text{Ps}}}}$ are P-wave velocities in the fast and slow velocity directions, respectively.

For a general P-wave ray with an incidence angle i, Eq. (1) becomes (e.g., Eberhart-Phillips and Henderson 2004):

$$S_{P} = S_{P0} - \sqrt {A^{2} + B^{2} } \cdot \cos^{2} i + \left( {A\cos 2\phi + B\sin 2\phi } \right) \cdot \sin^{2} i$$

(4)

for media with hexagonal symmetry having a slow symmetry axis oriented horizontally (Fig. 1a). Thus, the observation equation for weak P-wave AAN tomography is:

$$r = \frac{\partial T}{{\partial \varphi_{e} }} \cdot {\Delta }\varphi_{e} + \frac{\partial T}{{\partial \lambda_{e} }} \cdot {\Delta }\lambda_{e} + \frac{\partial T}{{\partial h_{e} }} \cdot {\Delta }h_{e} + {\Delta }T_{e} + \mathop \sum \limits_{n = 1}^{N} \left( {\frac{\partial T}{{\partial S_{n} }} \cdot {\Delta S}_{n} } \right) + \mathop \sum \limits_{m = 1}^{M} \left( {\frac{\partial T}{{\partial A_{m} }} \cdot {\Delta }A_{m} + \frac{\partial T}{{\partial B_{m} }} \cdot {\Delta }B_{m} } \right)$$

(5)

where $\partial T/\partial {\varphi }_{e}$, $\partial T/\partial {\lambda }_{e}$ and $\partial T/\partial {h}_{e}$ are partial derivatives of travel time (T) with respect to the hypocentral parameters of a local earthquake, and $\partial T/\partial {S}_{n}$ is partial derivative of travel time with respect to slowness at each node (1 … N) of a 3-D grid arranged in the study volume for 3-D isotropic Vp structure, $\partial T/\partial {A}_{m}$ and $\partial T/\partial {B}_{m}$ are partial derivatives of travel time with respect to the anisotropic parameters (A and B) at each node (1 … M) of another 3-D grid arranged in the study volume for Vp anisotropy (Wang and Zhao 2013). Note that the grid interval for the anisotropic Vp structure should be greater than that for the isotropic Vp structure, because Vp anisotropy is usually resolved with lower resolution than that of isotropic Vp structure. In a practical study, usually tens of thousands of or even over one million travel-time data are used to conduct 3-D anisotropic tomography beneath a seismic network. Thus, the observation Eq. (5) forms a huge but sparse system of linear equations, which can be resolved using an efficient matrix solver such as the LSQR algorithm (Paige and Saunders 1982).

Zhao et al. (1992, 1994) developed a tomographic inversion method for studying 3-D isotropic Vp and Vs structures, which can deal with complex-shaped seismic discontinuities in a study volume (e.g., the Conrad and Moho discontinuities, and a subducting slab boundary) and can use a great number of P and S wave arrival times and later phase data (i.e., reflected or converted waves at velocity discontinuities) from local earthquakes and teleseismic events. In the past three decades, this tomographic method has been widely applied to many regions for studying 3-D velocity structures from local, regional to global scales (see Zhao 2015 for a detailed review). Wang and Zhao (2008, 2013) modified the tomographic codes of Zhao et al. (1992, 1994) to study Vp anisotropic tomography.

In the past one and a half decades, a number of researchers have used the AAN tomographic methods to study the 3-D Vp structure and dynamics of many subduction zones, including New Zealand (Eberhart-Phillips and Henderson 2004; Eberhart-Phillips and Reyners 2009), NE Japan (Wang and Zhao 2008; Cheng et al. 2011; Ishise et al. 2018, 2021), SW Japan (Ishise and Oda 2008; Wang and Zhao 2012, 2013), South Kuril (Koulakov et al. 2015a; Niu et al. 2016; Gou et al. 2019b; Liu et al. 2013), NW Pacific (Wei et al. 2015; Ma et al. 2019; Liang et al. 2022), Taiwan–Philippine (Koulakov et al. 2015b; Fan and Zhao 2019, 2021a), SE Asia (Koulakov et al. 2009; Hua et al. 2022; Wei et al. 2022), Alaska (Tian and Zhao 2012; Gou et al. 2019a), Cascadia (Zhao and Hua 2021), Central America (Rabbel et al. 2011), Chile (Huang et al. 2019), the Eastern Mediterranean and Middle East (Wei et al. 2019), and Alps (Hua et al. 2017; Rappisi et al. 2022).

Figure 2 shows a recent example of subduction-zone AAN tomography, which was determined by a joint inversion of local and teleseismic travel-time data recorded by dense seismic networks deployed on the Cascadia margin (Zhao and Hua 2021). In the crust of the overriding North American plate, the FVDs are generally trench-parallel, which reflects N-S compression along the Cascadia margin. In the subducting slab under the forearc, the FVDs are generally trench-parallel, suggesting that the slab anisotropy may be modified at the outer rise by serpentinization of outer-rise faults (e.g., Faccenda et al. 2008). In northern Cascadia, a prominent low-velocity (low-V) zone below the slab exhibits NE-SW FVDs, reflecting thermally buoyant mantle flows from nearby oceanic hotspots, which move toward the northeast and gradually accumulate under northern Cascadia, leading to decompression melting.

Fig. 2

a–j Map views of P-wave azimuthal anisotropy (AAN) tomography of the Cascadia subduction zone (Zhao and Hua 2021). The layer depth is shown atop each map. Red and blue colors denote low and high isotropic Vp perturbations, respectively. The orientation and length of the white bars represent the AAN fast-velocity direction (FVD) and amplitude, respectively. The scales for the isotropic Vp perturbation and the AAN amplitude are shown at the bottom. Red triangles in a, b denote Cascadia volcanoes. The black solid and dashed lines denote locations of the upper and lower boundaries, respectively, of the subducting Juan de Fuca slab at each depth. k Tectonic setting of the Cascadia subduction zone. The arrows and associated annotations show the motion direction and rate of the Juan de Fuca plate relative to the North American plate. The red triangles denote the Cascadia volcanoes. The green crosses denote seismicity. The red dots denote three hotspots (Bowie, Cobb and Anahim). The red box shows the study region. The black dashed lines denote approximate boundaries of Northern Cascadia (NC), Central Cascadia (CC) and Southern Cascadia (SC). VI, Vancouver Island; OP, Olympic Peninsula; PS, Puget Sound; CR, Cascade Range; OCR, Oregon Coast Range; WV, Willamette Valley; HLP, High Lava Plain; KM, Klamath Mountains

Two simplified endmember models of mantle flow in subduction zones have been proposed based on a great number of studies of seismic anisotropy and numerical modeling (e.g., Long 2016). In one model (Fig. 3a), 2-D corner flow in the mantle wedge and 2-D entrained flow below the slab are induced by viscous coupling between the subducting slab and the ambient mantle. In another model (Fig. 3b), a 3-D flow field is induced by trench migration and slab rollback, and the flow is deflected horizontally around a slab edge or slab tear, producing trench-parallel flow beneath the slab and focused trench-normal flow in the mantle wedge beneath the arc and back-arc areas. The two models are very simplified, and the actual flow field in any subduction zone may reflect a combination of the two (Long 2016). The AAN tomographic images of the Cascadia margin (Fig. 2) and other subduction zones as discussed below may reflect joint effects of the subduction-driven flow (Fig. 3a) and toroidal mantle flow around the slab edge (Fig. 3b).

Fig. 3

Cartoon sketch of two simplified endmember models for subduction-zone mantle flow: a 2-D entrained flow model and b 3-D toroidal flow model. In a, viscous coupling between the downgoing slab and surrounding mantle leads to 2-D corner flow in the mantle wedge and 2-D entrained flow beneath the slab (black arrows). Black triangles on the overriding plate indicate arc volcanoes. In b, the trench migration (white arrow) and slab rollback induce a 3-D flow field (black arrow), in which flow is deflected horizontally around the slab edge, leading to a significant component of trench-parallel flow beneath the slab and focused flow parallel to the trench migration (i.e., generally perpendicular to the trench) in the central portion of the wedge back-arc. These sketches indicate highly simplified endmember models; the actual flow field in any given subduction zone likely reflects a combination of the two. After Long (2016)

An excellent example of toroidal mantle flow around the slab edge is that occurring in North Kamchatka (Fig. 4). The junction of the Aleutian Island and the Kamchatka peninsula defines a sharp turn in the boundary of the Pacific and North American plates, terminating the NW Pacific subduction zones. The regional SWS pattern near the junction shows that trench-parallel strain follows the Wadati–Benioff deep seismic zone but rotates to trench-normal beyond the slab edge, suggesting that asthenospheric mantle flows around and below the disrupted slab edge and may influence the shallowing dip of the Wadati–Benioff zone at the Aleutian junction (Peyton et al. 2001). Zhao et al. (2021) determined 3-D Vs AAN tomography down to ~300 km depth beneath Kamchatka using teleseismic Rayleigh-wave data (Fig. 4a–e). The subducting Pacific slab is clearly imaged as a dipping high-Vs zone beneath South Kamchatka, whereas this high-Vs zone is absent in the north, suggesting that the slab edge exists under North Kamchatka where low-Vs anomalies are revealed (Fig. 4a–e). The slab edge was also revealed by isotropic Vp tomography (Fig. 4f; Jiang et al. 2009). The Pacific slab exhibits an NE–SW FVD in the area south of ~54° N latitude, whereas the dominant FVD rotates to ~N–S or ~NNW–SSE in the area of ~54°–57° N latitude (Fig. 4e). The FVDs in the slab are generally parallel to the Kamchatka trench, which may reflect SPO of the stratified oceanic lithosphere. At shallow depths (< ~150 km), the low-Vs mantle wedge exhibits ~NW–SE FVDs, whereas ~E–W FVDs appear around the slab edge (Fig. 4a–e) which may reflect a combination of corner flow in the mantle wedge (Fig. 3a) and toroidal mantle flow around the slab edge (Fig. 3b). The toroidal flow may deform, heat or even melt the slab near its edge (Zhao et al. 2021).

Fig. 4

a–e Map views of S-wave velocity (Vs) azimuthal anisotropy tomography beneath Kamchatka (Zhao et al. 2021). The layer depth is shown at the lower-right corner of each map. The red and blue colors denote low and high isotropic Vs perturbations, respectively. The orientation and length of white bars represent the fast-velocity direction (FVD) and the anisotropic amplitude, respectively. The scales for the isotropic Vs perturbations and the anisotropic amplitude are shown beside (c). The white circles denote local seismicity within a ± 10 km depth range of each layer. The thick purple line on each map denotes the upper boundary of the subducting Pacific slab at each depth. The purple dashed line shows the 110 Ma isochron of the incoming Pacific plate. The gray arrow denotes the moving direction of the incoming Pacific plate. The red triangles denote active volcanoes. The black sawtooth lines denote oceanic trenches. f A cartoon of the Pacific slab edge beneath North Kamchatka (Jiang et al. 2009)

P-wave AAN tomography has been also applied to study the 3-D structure and dynamics beneath continental regions (Figs. 5, 6 and 7), including Mainland China (Tian and Zhao 2013; Huang et al. 2014; Wei et al. 2016; Jiang et al. 2021; Jia et al. 2022), the Tibetan Plateau (Wei et al. 2013; Zhang et al. 2017; Yang et al. 2022a), western and central USA (Huang and Zhao 2013; Yu and Zhao 2018; Wang and Zhao 2019a; Wang et al. 2022a; Yang et al. 2022b), Africa (Yu et al. 2020a), Turkey (Wang et al. 2020), and Greenland (Toyokuni and Zhao 2021). Figure 5 shows Vp AAN tomography beneath China (Wei et al. 2016), revealing that the northern limit of the subducting Indian plate has reached the Jinsha River suture in eastern Tibet. The FVD is NE–SW beneath the Indian continent, while the FVD is arc-parallel beneath the Himalaya and Tibetan Plateau, which may reflect re-orientation of minerals due to lithospheric extension, in response to the India-Eurasia collision. Multiple layers with variable anisotropies may exist in some parts of the Tibetan Plateau, resulting in null SWS in those areas (e.g., Chen et al. 2010). Around the Philippine Sea (PHS) slab beneath SE China, FVDs exhibit a circular pattern that may reflect asthenospheric strain due to toroidal mantle flow around the PHS slab edge (Fig. 5e, f).

Fig. 5

Map views of Vp azimuthal anisotropy tomography of China and surrounding regions (Wei et al. 2016). The red triangles denote active intraplate volcanoes. The purple lines denote sutures and major tectonic boundaries. Other labels are the same as those in Fig. 2

Fig. 6

a–d Map views of Vp azimuthal anisotropy tomography of the NE Tibetan Plateau (the blue box on the inset map in e). Yellow lines denote active faults. Thick gray lines denote major block boundaries. Other labels are the same as those in Fig. 2. e Tectonic background of NE Tibet. The yellow star denotes the 1920 Haiyuan earthquake (M8.5) epicenter. Yellow and white lines denote active faults and boundaries of tectonic blocks, respectively: the Alxa block (ALXB), Ordos basin (ODB), eastern boundary zone of Tibetan Plateau (AEB), Qilian block (QLB), eastern Kunlun‐Qaidam block (EKQB), Bayan Har block (BYHB), Qinling fold zone (QLFZ), and Sichuan basin (SCB). Sepia arrows denote GPS vectors. A gray scale of the surface topography is shown at the bottom. Black arrows on the inset map show the crustal channel flow. After Sun and Zhao (2020)

Fig. 7

Map views of Vp azimuthal anisotropy tomography of Southern California (Yu and Zhao 2018). The red circles and the red thick line in d–f depict a circular pattern of azimuthal anisotropy and deflected mantle flow, respectively. Black lines denote major faults. Red triangles are Cenozoic volcanoes. Red stars denote large earthquakes (M > 6.0) that occurred from 1933 to 2017. GV, Great Valley; PR, Peninsular Ranges; SN, Sierra Nevada; ST, Salton Trough; ECSZ, Eastern California Shear Zone; IB, Inner Borderland. Other labels are the same as those in Fig. 2

Prominent low-V anomalies are revealed in the upper mantle under NE Asia where the dominant FVD is NW–SE (Fig. 5). This feature reflects mantle flows in the big mantle wedge (BMW) above the northwestward subducting Pacific slab and the stagnant slab in the mantle transition zone (MTZ) under East Asia (Zhao et al. 2004, 2009; Lei and Zhao 2005; Chen et al. 2017; Liang et al. 2022; Cloetingh et al. 2022). The hot and wet upwelling flows and related geodynamic processes in the BMW have caused the destruction of the North China craton, back-arc spreading and the Japan Sea opening, as well as intraplate seismic and volcanic activities in East Asia (see a recent review by Zhao 2021).

Recently, Sun and Zhao (2020) revealed more detailed AAN features beneath the NE Tibetan Plateau (Fig. 6). They found extensive low‐V anomalies in the middle-lower crust, which may account for most of crustal anisotropy. The dominant FVDs are correlated with the strain field revealed by GPS observations and focal mechanism solutions in the transition zones among the Ordos basin, the Alxa block and the Tibetan Plateau. This feature is attributed to regional crustal flow intruding northeastward into NE Tibet and possibly affecting vertical ground motions. The flow has been blocked by the surrounding rigid areas and so failed to further extrude eastward between the Sichuan basin and the Ordos basin. The crustal flow may cause intra-crustal and crust‐mantle decoupling in the transition zones of NE Tibet. In the Helan–Liupan area in west-central China, ductile flow in the middle-lower crust was also revealed, and fluids in the flow affected the generation of large crustal earthquakes there (Cheng et al. 2016).

Significant depth-dependent anisotropy is revealed beneath Southern California (Fig. 7). FVDs in the lithospheric mantle closely follow the strike of the San Andreas fault, whereas FVDs in the asthenosphere exhibit a circular pattern centered in the high-V Isabella anomaly under the Great Valley (Yu and Zhao 2018). The Isabella anomaly may be a remnant of the fossil Farallon slab and is experiencing lithospheric downwelling now, contributing to the development of a circular asthenospheric flow. High-V anomalies exist below 300 km depth under areas surrounding the Great Valley, probably reflecting delaminated lithospheric segments. Different rifting processes may be occurring in the Salton Trough and the Inner Borderland, which may be related to lithospheric stretching and regional mantle upwelling, respectively.

3 Radial Anisotropy Tomography

Assuming that the hexagonal symmetry axis is vertical in the medium under study (Fig. 1b), Eq. (1) can be rewritten as a function of the ray incidence angle i (Wang and Zhao 2013):

$$S_{P} = S_{P0} + M\cos 2i$$

(6)

where M is a radial anisotropy (RAN) parameter. The anisotropic amplitude $\beta$ is then defined as:

$$\beta = \frac{{V_{ph} - V_{pv} }}{{2V_{P0} }} = \frac{M}{{S_{P0} - M^{2} /S_{P0} }}$$

(7)

where Vph and Vpv denote horizontally propagating Vp and vertically propagating Vp, respectively. Thus, positive RAN (i.e., $\beta$ > 0) means Vph > Vpv, and vice versa.

Similar to the case of AAN tomography, the observation equation for weak P-wave RAN tomography is:

$$r = \frac{\partial T}{{\partial \varphi_{e} }} \cdot {\Delta }\varphi_{e} + \frac{\partial T}{{\partial \lambda_{e} }} \cdot {\Delta }\lambda_{e} + \frac{\partial T}{{\partial h_{e} }} \cdot {\Delta }h_{e} + {\Delta }T_{e} + \mathop \sum \limits_{n = 1}^{N} \left( {\frac{\partial T}{{\partial S_{n} }} \cdot {\Delta S}_{n} } \right) + \mathop \sum \limits_{k = 1}^{K} \left( {\frac{\partial T}{{\partial M_{k} }} \cdot {\Delta }M_{k} } \right)$$

(8)

The right-hand side of Eq. (8) includes terms for local-earthquake hypocentral parameters, isotropic Vp perturbations, and anisotropic Vp parameters (M_k) at 3-D grid nodes arranged in a study volume. Then the system of observation Eq. (8) can be resolved (Wang and Zhao 2013), in a way similar to that of AAN tomography as described in the last section.

During the past decade, the RAN tomographic method has been applied extensively to study the 3-D Vp structure and dynamics of the crust and upper mantle in many regions, including NE Japan (Wang and Zhao 2013; Ishise et al. 2018), SW Japan (Wang and Zhao 2013), South Kuril (Niu et al. 2016), Eastern China (Wang et al. 2014; Jiang et al. 2021; Liang et al. 2022), Tibet (Zhang et al. 2021), the Philippines (Fan and Zhao 2019), Alaska (Gou et al. 2019a), Cascadia (Zhao and Hua 2021), Southern California (Yu and Zhao 2018; Yu et al. 2020b), Alps (Hua et al. 2017), Africa (Yu et al. 2020a), Greenland (Toyokuni and Zhao 2021), and Antarctica (Zhang et al. 2020).

Figure 8 shows a recent example of subduction-zone RAN tomography beneath the Cascadia margin (Zhao and Hua 2021). Negative RAN (i.e., V_pv > V_ph) is revealed in the crust and mantle wedge under the Cascadia volcanoes and back-arc area, which may reflect hot and wet upwelling flow associated with fluids from dehydration reactions of the young and warm Juan de Fuca-Gorda plate that is descending northeastward under the North American plate. A similar feature of negative RAN was revealed under arc volcanoes in Japan (Wang and Zhao 2013). Note that the relationship between RAN and upwelling of hot materials depends on the type of LPO. Negative RAN would correspond to vertical flow for LPO developed under relatively dry conditions (A, D and E type LPO) but not for B and C type LPO (Karato et al. 2008). The RAN patterns in the subducted Juan de Fuca slab and the subslab mantle are nearly the same (Fig. 8), i.e., negative RAN occurs under the Cascadia volcanoes and the back-arc, whereas positive RAN appears beneath the Cascadia forearc. This feature may indicate strong coupling of the subducting slab with the subslab asthenosphere, and that the entrainment of asthenosphere with the overriding slab results in the subslab mantle flow. In northern Cascadia, a low-V zone appears in the subslab mantle, and it exhibits negative RAN and NE–SW FVDs (Figs. 2 and 8). The subslab low-V zone may reflect thermally buoyant mantle flows from nearby oceanic hotspots (Cobb, Bowie and Anahim), which may affect the nucleation of megathrust earthquakes on the slab interface (Bodmer et al. 2018; Fan and Zhao 2021b). The hot and buoyant mantle materials are transported toward the northeast by the overriding slab and gradually accumulate beneath northern Cascadia, leading to decompression melting. A subslab low-V zone exhibiting positive RAN and NE–SW FVDs is also revealed under southern Cascadia (Figs. 2 and 8), reflecting that the asthenosphere there is entrained by the overriding slab.

Fig. 8

East–west vertical cross sections of Vp radial anisotropy (RAN) tomography of the Cascadia subduction zone along the profiles shown on the inset map (Zhao and Hua 2021). Red and blue colors denote low and high isotropic Vp perturbations, respectively. Horizontal and vertical bars represent positive and negative RAN results, respectively. Scales for the isotropic Vp perturbation and the RAN amplitude are shown at the bottom-right. The surface topography along each profile is shown above each panel. Open circles denote local seismicity within a 10-km width of each profile. Red triangles denote the Cascadia volcanoes. White lines denote the Moho discontinuity and the slab upper boundary. The white dashed line denotes the estimated lower boundary of the slab

Gou et al. (2019a) determined high-resolution 3-D Vp AAN and RAN images down to 900 km depth beneath Alaska by inverting a large number of arrival-time data from local earthquakes and teleseismic events (Fig. 9). Their results show that the high-Vp Pacific slab has subducted down to 450–500 km depths. A prominent slab gap is revealed at depths of 65–120 km near the Wrangell volcanic field, which is likely a slab tear acting as a channel that provides ascending mantle materials to generate magmas feeding the surface volcanoes. In the back-arc mantle wedge near the eastern slab edge, the AAN exhibits trench-parallel FVDs, which may reflect along-strike mantle flow. The FVDs in the subducting Pacific slab are nearly east–west, which may indicate fossil anisotropy formed at the mid-ocean ridge. Negative RAN is revealed within the subducting slab, which may be caused by the fast plate subduction with a steep dip angle. Trench-normal FVDs of the AAN are revealed in the mantle below the Pacific slab, which may reflect mantle flow entrained by the subducting slab. Positive RAN is revealed in the mantle beneath the Yakutat slab, indicating that its shallow subduction flattens the mantle flow below the slab to be subhorizontal. Along-strike FVDs of the AAN around the eastern slab edge may indicate the edge-induced toroidal mantle flow.

Fig. 9

Vertical cross sections of Vp radial anisotropy (RAN) tomography of the Alaska subduction zone along the profiles shown on the inset map (Gou et al. 2019a). Red triangles denote active volcanoes. Two red dashed lines denote the 410 and 660 km discontinuities. Other labels are the same as those in Fig. 8

Hua et al. (2017) determined the first tomographic images of 3-D Vp AAN and RAN in the crust and upper mantle beneath the Alps by conducting a joint inversion of arrival-time data of local earthquakes and teleseismic events (Fig. 10). Their results show the southward dipping European plate with a high-V anomaly beneath the western-central Alps and the northward dipping Adriatic plate with a high-V anomaly beneath the Eastern Alps, indicating that the subduction polarity changes along the strike of the Alps. The Vp AAN is characterized by mountain chain-parallel FVDs in the western-central Alps and NE-SW FVDs in the eastern Alps, which may be caused by mantle flow induced by the slab subductions. Their results show negative RAN (i.e., V_ph < V_pv) within the subducting slabs and positive RAN (i.e., V_ph > V_pv) in the low-V mantle wedge, which may reflect the subvertical plate subduction and its induced mantle flow. The results of anisotropic tomography provide important new information on the complex mantle structure and dynamics of the Alps and adjacent regions.

Fig. 10

Vertical cross sections of Vp radial anisotropy tomography of the Alps region along the profiles shown on the inset map (Hua et al. 2017). Other labels are the same as those in Fig. 9

4 Joint Inversion for P and S Wave Anisotropy

Both P and S wave arrival times of earthquakes are abundant and precise data in seismology. To date, most body-wave anisotropic studies have used only P-wave travel-time data, whereas S-wave arrival times are seldom used. Liu and Zhao (2016b) presented 3-D images of AAN tomography of the Japan subduction zone (Fig. 11), which are determined by using both P and S wave arrival-time data of local earthquakes and teleseismic events recorded by the dense seismic networks deployed on the Japan Islands. They modified the inversion method for Vp AAN tomography (Wang and Zhao 2013) and extended it to invert S-wave travel times for 3-D Vs AAN tomography. They also performed a joint inversion of the P and S wave data to constrain the 3-D AAN structure beneath Japan.

Fig. 11

Map views of Vp azimuthal anisotropy tomography along three slices in a the middle mantle wedge, b the subducting Pacific slab, and c the subslab mantle beneath NE Japan as shown in the sketch beside (b). This 3-D Vp model is obtained by a joint inversion of P and S wave travel-time data of both local and teleseismic events. The slice (a) is located at the middle between the Moho and the slab upper boundary. Slices (b) and (c) are located at 50 km and 200 km depth below the slab upper boundary, respectively. Red and blue colors denote low and high isotropic Vp perturbations, respectively, whose scale (in %) is shown below (c). The orientation and length of the black bars represent the horizontal FVD and anisotropic amplitude, respectively, whose scale is shown below (c). The white lines to the east of the Japan Trench denote isochrones of the Pacific Ocean floor. The white arrow in each map shows the moving direction of the Pacific plate relative to NE Japan. Red triangles denote active arc volcanoes. After Liu and Zhao (2016b)

Different from the previous AAN tomographic studies (e.g., Wang and Zhao 2008, 2013) assuming hexagonal anisotropy with a horizontal symmetry axis, Liu and Zhao (2016b) assumed that the medium under study exhibits weak orthorhombic anisotropy that has three orthogonal symmetry axes, two of which are located in the horizontal symmetry plane. When the two horizontal symmetric axes represent the fast and slow velocity directions, for any P-wave ray with an incident angle i, Eq. (1) can be rewritten as:

$$S_{P} = S_{P0} + C \cdot \cos^{2} i + \left( {A\cos 2\phi + B\sin 2\phi } \right) \cdot \sin^{2} i$$

(9)

where C is a constant and represents the anisotropic strength in the vertical direction. When C = 0, the slowness in the vertical direction is equal to the isotropic slowness, then for any S-wave ray with an incident angle i, its slowness Ss can be approximately expressed as:

$$S_{S} = S_{S0} + \frac{1}{2}\left( {\sin^{2} i + \cos^{2} j \cdot \cos^{2} i - \sin^{2} j} \right) \cdot \left( {A\cos 2\phi + B\sin 2\phi } \right)$$

(10)

where S_S0 is isotropic S-wave slowness, and j is the angle between the S-wave and SV-wave polarization directions, both of which have the same propagating direction. In the practical analysis, j is hard to determine precisely, so it is estimated with a grid-search method finding its optimal value that leads to the minimum root-mean-square travel-time residual (Liu and Zhao 2016b). Because the same AAN parameters A and B appear in Eqs. (9) and (10), a joint inversion of P and S wave travel-time data was conducted to constrain the 3-D AAN structure. The system of observation equations is similar to Eq. (5), and the FVD and anisotropic amplitude of AAN are estimated using equations similar to those in Eqs. (2) and (3).

The results of Liu and Zhao (2016b) show that the subducting Pacific and PHS slabs beneath the Japan Islands mainly exhibit trench-parallel FVDs (Fig. 11b), which may reflect frozen-in LPO of aligned anisotropic minerals formed at the mid-ocean ridge as well as SPO due to hydrated faults in the oceanic plate produced at the outer-rise area near the trench axis (e.g., Faccenda et al. 2008; Wang et al. 2022b). Trench-normal FVDs are generally revealed in the mantle wedge (Fig. 11a), which reflect corner flows in the mantle wedge due to the oceanic plate subduction and dehydration. These results are generally consistent with those of the earlier studies (Wang and Zhao 2008, 2013). Roughly trench-normal FVDs exist in the subslab mantle (Fig. 11c), possibly reflecting asthenospheric shear deformation related to the overlying slab subduction, similar to that in the subslab mantle beneath the Cascadia margin (Fig. 2).

Liu and Zhao (2016b) also revealed toroidal mantle flows in and around a hole (window) in the subducting PHS slab under SW Japan that was firstly found by isotropic Vp and Vs tomography (Zhao et al. 2012; Asamori and Zhao 2015). This result suggests that the PHS slab window may have led to a complex flow pattern in the mantle wedge above the Pacific slab. As mentioned above, a similar circular FVD pattern appears around the PHS slab in SE China (Fig. 5e, f), reflecting asthenospheric strain caused by toroidal mantle flow around the PHS slab edge (Wei et al. 2016).

5 Joint Inversion for Azimuthal and Radial Anisotropy

Due to the limitations of data sets used, most P-wave anisotropic studies perform separate inversions for either Vp AAN tomography or Vp RAN tomography, because ray paths in a data set usually do not have complete coverages in both azimuth and incident angle. Taking advantage of abundant P-wave travel-time data of both local shallow and deep earthquakes as well as teleseismic events recorded by the dense seismic networks on the Japan Islands, Liu and Zhao (2017) attempted to perform a joint inversion for 3-D Vp AAN and RAN variations, as well as isotropic Vp tomography, by assuming orthorhombic anisotropy with a vertical symmetry axis and two horizontal symmetry axes in their study volume of the Japan subduction zone.

As mentioned in the last section, for any P-wave ray with an incident angle i in a medium exhibiting weak orthorhombic anisotropy, its slowness can be expressed in Eq. (9). Liu and Zhao (2016b) assumed C = 0 in Eq. (9), and inverted for only the AAN parameters A and B, as well as 3-D isotropic Vp variations. In contrast, Liu and Zhao (2017) took C as an unknown parameter representing RAN, and made a joint inversion for the three anisotropic parameters (A, B and C) in addition to isotropic Vp variation. After the joint inversion, the AAN FVD and the anisotropic amplitudes of both AAN and RAN were determined using Eqs. (2), (3) and (7).

The result of Liu and Zhao (2017) shows that the AAN tomography determined by the joint inversion (Fig. 12d–f) is quite consistent with that obtained by inverting for only the two AAN parameters (A and B) (Fig. 12a–c). This result suggests that the P-wave ray paths used have a quite complete azimuthal coverage, hence the AAN tomography is obtained stably. However, the RAN tomography determined by the joint inversion (Liu and Zhao 2017) is somehow different from that obtained by inverting for only the RAN parameter C (e.g., Wang and Zhao 2013). This difference in the RAN tomography is mainly due to the still incomplete ray coverage in the ray incident angle, so the inversion becomes less stable because of the increase of the unknown parameters (see a detailed discussion by Liu and Zhao 2017). These results indicate that the data quality (including both the data picking accuracy and the ray-path coverage and crisscrossing) determines the quality of a final tomographic result (e.g., Zhao 2015). Only when a good enough data set is used, a reliable tomographic result can be obtained.

Fig. 12

Map views of isotropic Vp image and Vp azimuthal anisotropy obtained by conducting a–c Vp azimuthal anisotropy tomography (the VA inversion) and d–f a joint inversion for 3-D isotropic Vp and 3-D Vp azimuthal and radial anisotropies (the VAR inversion). The layer depth is shown at the upper-left corner of each map. The blue line shows location of the upper boundary of the subducting Pacific slab at each depth. Red and blue colors denote low and high isotropic Vp perturbations, respectively, whose scale (in %) is shown in (a). Orientation and length of the black bars represent the horizontal FVD and azimuthal anisotropy amplitude (α_P), respectively. The α_P scale is shown in (d). Red triangles denote active volcanoes. Solid sawtooth lines and the black dashed line denote plate boundaries at the surface. After Liu and Zhao (2017)

Because olivine, the most abundant mineral in the mantle, is a kind of orthorhombic mineral, the Vp tomographic results obtained with the orthorhombic anisotropy assumption (Liu and Zhao 2017) may well reflect the mineral fabric patterns in the upper mantle. According to a simple flow field and the types of olivine fabric classified by mineral physics studies (e.g., Karato et al. 2008), Liu and Zhao (2017) used their new models of Vp AAN and RAN to estimate the distribution of olivine fabrics in the mantle wedge above the subducting Pacific and PHS slabs beneath Japan. Their results show that the forearc mantle wedge above the Pacific slab under NE Japan exhibits AAN with trench-parallel FVDs and Vhf > Vv > Vhs (here Vv is Vp in the vertical direction, Vhf and Vhs are Vp values in the fast and slow directions in the horizontal plane, respectively), where B-type olivine fabric with vertical trench-parallel flow may dominate. However, this anisotropic feature is not clear in the forearc mantle wedge above the PHS slab beneath SW Japan, possibly because of higher temperatures and more fluids there associated with the young and warm PHS slab subduction (e.g., Salah and Zhao 2003; Asamori and Zhao 2015; Hua et al. 2018). Trench-normal FVDs and Vhf > Vv > Vhs generally occur in the mantle wedge under the arc and back-arc in Japan, where E-type olivine fabric with FVD-parallel horizontal flow may dominate. Under western Honshu where four potentially active volcanoes exist (e.g., Hua et al. 2018; Zhao et al. 2018), the mantle wedge exhibits an anisotropic feature of Vv > Vhf > Vhs, so C-type olivine fabric may dominate there. This result suggests that the water content is high in the mantle wedge beneath western Honshu, because the Pacific and PHS slabs coexist there and their dehydration reactions release abundant fluids to the overlying mantle wedge (e.g., Zhao et al. 2012; Asamori and Zhao 2015).

The general AAN features of the Japan subduction zone revealed by the P-wave tomographic studies (Figs. 11 and 12) are confirmed by a recent S-wave AAN study. Liu and Zhao (2016a) determined 3-D Vs AAN tomography down to ~ 300 km depth beneath the Japan Islands and the Japan Sea (Fig. 13) using teleseismic Rayleigh-wave data recorded by the dense Japanese seismic networks. The subducting Pacific slab is revealed clearly as a dipping high-Vs zone with trench-parallel FVDs, which may indicate the slab anisotropy arising from hydrated faults produced at the outer-rise area near the Japan trench axis, overprinting the slab fossil fabric (e.g., Wang and Zhao 2021; Wang et al. 2022b), whereas the mantle wedge generally exhibits lower Vs with trench-normal FVDs which reflect subduction-driven corner flow and anisotropy. A similar 3-D Vs AAN result was obtained recently beneath the Mariana arc and back-arc using teleseismic Rayleigh-wave data (Qiao et al. 2021).

Fig. 13

Map views of Vs azimuthal anisotropy tomography beneath the Japan Islands and the Japan Sea determined by inverting Rayleigh-wave dispersion data (Liu and Zhao 2016a). The layer depth and average Vs at the depth are shown at the upper-left corner of each map. Red and blue colors denote low and high isotropic Vs perturbations, respectively. The orientation and length of black bars represent the fast-velocity direction (FVD) and the anisotropic amplitude, respectively. Their scales are shown at the bottom. The two white arrows in each map denote moving directions of the Pacific and Philippine Sea plates. Red solid and dashed lines denote plate boundaries. Red triangles denote active volcanoes

Depth variations of AAN are also revealed in the big mantle wedge (BMW) under the Japan Sea (Fig. 13), which may reflect past deformations in the Eurasian lithosphere related to the back-arc spreading during 21–15 Ma and complex current convection in the BMW induced by active subductions of both the Pacific and PHS plates (Zhao et al. 2004; Liu et al. 2017; Liang et al. 2022).

6 Inversion for Tilting-Axis Anisotropy

In the heterogeneous and dynamic Earth’s interior, inclined hexagonal symmetry is not unusual (e.g., Plomerova et al. 2008, 2012; Eken et al. 2012). Levin and Park (1998) investigated the influence of 1-D anisotropy, its type, layering and orientation, on the coda of teleseismic P waves, particularly on P-SH converted waves. They employed a reflectivity technique to compute the transmission response of a flat-layered medium with arbitrarily oriented hexagonally symmetric anisotropy, and showed main features of synthetic seismograms for a variety of simple models and discussed how to interpret converted phases in observations. Park and Levin (2016) applied harmonic stacking to migrated receiver functions from 471 teleseismic events recorded at a station in Saudi Arabia. They found a strong two-lobed back-azimuth variation in Ps converted wave from a Moho near 40 km depth and from an interface 25–30 km deeper, both consistent with a north-striking tilted axis of symmetry. The identical polarity of these two signals indicates that they are not Ps conversions from the top and bottom of a single subcrustal anisotropic layer. These results confirm the general features reported by Levin and Park (1998).

Seismic anisotropy with inclined hexagonal symmetry may occur particularly in subduction zones, due to the existence of a dipping slab and 3-D mantle flow around the slab (Fig. 3). As mentioned above, the PWAT studies made so far generally assumed the symmetry axis in the anisotropic medium to be either horizontal or vertical to estimate Vp azimuthal or radial anisotropy (Fig. 1). Thus, the real anisotropy in the subducting slab or the mantle around the slab is hard to clarify, because the anisotropic symmetry axis in those areas may be neither horizontal nor vertical (Fig. 3).

To resolve this issue, Munzarova et al. (2018a) developed a PWAT code that inverts relative travel-time residuals of teleseismic P waves for 3-D distribution of isotropic Vp perturbations and Vp anisotropy in the upper mantle. They considered weak hexagonal anisotropy with the symmetry axis oriented generally in 3-D. Munzarova et al. (2018b) used their code to determine an anisotropic Vp model of the upper mantle beneath northern Fennoscandia that is a tectonically stable Precambrian region with a thick anisotropic mantle lithosphere without significant thermal heterogeneities (Plomerova et al. 2012). Their model shows the strongest anisotropy and the largest Vp perturbations at depths corresponding to the mantle lithosphere, whereas the lateral variations are insignificant in the underlying asthenosphere. The retrieved domain-like anisotropic structure of the mantle lithosphere in northern Fennoscandia was attributed to preserved fossil fabrics of the Archean microplates accreted during the Precambrian orogenic processes (Munzarova et al. 2018b). Due to the lack of local seismicity, Munzarova et al. (2018a, 2018b) used only teleseismic P-wave data, so the ray-path coverage is less ideal in their study volume, particularly at shallow depths (< 50 km).

Recently, Wang and Zhao (2021) developed a new PWAT code to invert for 3-D Vp anisotropic structure under an assumption of general hexagonal-symmetry anisotropy. Because their code conducts a joint inversion of local and teleseismic travel-time data, their study volume is densely covered by crisscrossing ray paths, which can well constrain the complex 3-D anisotropic structure of the Japan subduction zone (Figs. 14 and 15). According to the assumption, a seismic wave propagates slowest along the hexagonal-symmetry axis (HSA) direction and propagates fastest along the directions normal to the axis. Thus, the FVDs constitute a plane (fast-velocity plane, FVP) that is parallel to fault planes in fractured rocks or the two faster directions of orthotropic olivine, depending on whether the anisotropy is caused by SPO or LPO (Fig. 15d).

Fig. 14

Updated Vp anisotropic tomography of East Japan with tilting symmetry axes (Wang and Zhao 2021). Map views at A 10 km, B 70 km and C 200 km depths. The blue line denotes the slab upper boundary at each depth. Red triangles denote active volcanoes. Background colors denote isotropic Vp perturbations whose scale is shown in A. T-shaped symbols denote anisotropy whose amplitude scale is shown in C. The cartoon in B shows how a dipping FVP (fast-velocity plane) is displayed on the 2-D paper plane (PP) as a T-shaped symbol. FVPs in and outside the subducting Pacific slab are shown in yellow and red, respectively

Fig. 15

a–c The same as Fig. 14 but vertical cross sections along the three profiles shown in Fig. 14b. The black curved line denotes the Moho discontinuity. Blue solid and dashed lines denote the slab upper and lower boundaries, respectively. Black crosses and green dots show background seismicity (M 2.0–4.0) and low-frequency microearthquakes (M 0.0–2.0), respectively, within a 20-km width of each profile. d A schematic diagram showing the relationship between the hexagonal symmetry axis (HSA) and fast-velocity plane (FVP) and causes of anisotropy. LPO, lattice-preferred orientation; SPO, shape-preferred orientation. After Wang and Zhao (2021)

In an anisotropic medium with a HSA, P-wave anisotropic slowness (Sp) can be expressed as

$$S_{P} = S_{P0} + M\cos 2\theta$$

(11)

where Sp₀ is the average P-wave slowness (i.e., the isotropic component), M is a positive parameter quantifying the anisotropic amplitude, and $\theta$ is the angle between the HSA and the ray propagation direction. Equation (11) is rewritten by Wang and Zhao (2021) as follows:

$$S_{P} = S_{P0} \left[ {1 + 2\left( {dU + eW + fR} \right)^{2} - \left( {U^{2} + W^{2} + R^{2} } \right)} \right]$$

(12)

where (d, e, f) is P-wave ray vector, and (U, W, R) is the HSA direction vector. The azimuth $\gamma$ and the angle q from the positive depth direction of HSA, and the anisotropic amplitude $\alpha$ can be expressed as follows:

$$\gamma = \left\{ {\begin{array}{*{20}l} {\frac{\pi }{2}} \hfill & {W = 0} \hfill \\ {\arctan \frac{U}{W}} \hfill & {W \ne 0} \hfill \\ \end{array} } \right.$$

(13)

$$q = \left\{ {\begin{array}{*{20}l} {\frac{\pi }{2}} \hfill & {R = 0} \hfill \\ {\arctan \frac{{\sqrt {U^{2} + W^{2} } }}{R}} \hfill & {R \ne 0} \hfill \\ \end{array} } \right.$$

(14)

$$\alpha = \frac{{V_{{{\text{Pf}}}} - V_{{{\text{Ps}}}} }}{{2V_{P0} }} = \frac{{U^{2} + W^{2} + R^{2} }}{{1 - \left( {U^{2} + W^{2} + R^{2} } \right)^{2} }}$$

(15)

To obtain a more accurate result, Wang and Zhao (2021) rewrote the observation equation by adopting the second-order Taylor series to express theoretical travel time as follows:

$$r = \frac{\partial T}{{\partial \varphi_{e} }} \cdot \Delta \varphi_{e} + \frac{\partial T}{{\partial \lambda_{e} }} \cdot \Delta \lambda_{e} + \frac{\partial T}{{\partial h_{e} }} \cdot \Delta h_{e} + \Delta T_{e} + \mathop \sum \limits_{n = 1}^{N} \left( {\frac{\partial T}{{\partial Q_{n} }} \cdot \Delta Q_{n} } \right) + \frac{1}{2}\mathop \sum \limits_{{n_{1} ,n_{2} }} \left( {\frac{{\partial^{2} T}}{{\partial Q_{{n_{1} }} \partial Q_{{n_{2} }} }} \cdot \Delta Q_{{n_{1} }} \Delta Q_{{n_{2} }} } \right)$$

(16)

where r is travel-time residual, N is the total number of model parameters defined at the 3-D grid nodes, Qn is the nth model parameter (i.e., Sp₀, U, W or R at the 3-D grid nodes). This approach can resolve the problem of inaccurate travel-time calculation due to the very nonlinearity of the anisotropic model and can greatly reduce the trade-off between the isotropic and anisotropic model parameters (Wang and Zhao 2021).

However, expressing the theoretical travel time with the second-order Taylor series leads to two problems. One is that the linear system of observation equations becomes nonlinear. As a result, the well-used tools to solve the linear system of equations, such as the LSQR (Paige and Saunders 1982) or the singular value decomposition (SVD) techniques, cannot be used to solve Eq. (16). The other problem is that the memory storage and computational cost are increased by one to two orders, making a single computer hard to handle the problem of 3-D anisotropic tomography. To resolve the first problem, Wang and Zhao (2021) adopted a nonlinear solver, the L-BFGS-B algorithm (Zhu et al. 1997), to solve Eq. (16). To resolve the second problem, they adopted parallel processing to improve their codes to realize high-performance computations.

Wang and Zhao (2021) applied their new PWAT code to determine an updated 3-D Vp anisotropic model of the Japan subduction zone by inverting more than one million travel-time data of local and teleseismic events recorded by the dense seismic networks on the Japan Islands, a high-quality data set similar to those used by Liu and Zhao (2016b, 2017). As a result, the 3-D anisotropic structure in and around the dipping Pacific and PHS slabs is revealed precisely (Figs. 14 and 15). Their results show that the deformation of the subducting slabs plays an important role in both mantle flow and intraslab fabric. An important new finding of their study is that the trench-parallel anisotropy in the eastern Japan forearc region observed by the SWS measurements (e.g., Huang et al. 2011a; Nakajima et al. 2006; Okada et al. 1995) and the above-mentioned Vp studies results, at least partly, from the intraslab deformation during outer-rise yielding of the subducting Pacific plate (Wang and Zhao 2021).

More recently, the new technique by Wang and Zhao (2021) was applied to arrival-time data of local earthquakes recorded at both onshore Hi-net and offshore S-net seismic stations in NE Japan to determine detailed Vp anisotropic tomography of the Tohoku forearc area (Wang et al. 2022b). The new results show that trench-parallel intraslab FVPs of anisotropy intersect the Pacific slab upper surface at high angles (~ 45°–90°), reflecting aligned hydrated faults in the slab. Ruptures of the hydrated faults may cause large intraslab earthquakes, such as the 2022 M7.4 Fukushima-Oki earthquake (Wang et al. 2023). The hydrated-fault associated water entering a large near-trench asperity in the megathrust zone could have triggered the great 2011 Tohoku-oki earthquake (Mw 9.0). These results suggest that outer-rise hydrated faults subducting with the oceanic lithosphere are very important for the water cycle in subduction zones, which play a critical role in seismogenesis and arc magmatism (e.g., Zhao et al. 2010; Yu and Zhao 2020).

7 Anisotropic 3-D Ray Tracing

In the above-mentioned PWAT studies, theoretical ray paths and travel times are generally calculated with a 3-D ray-tracing technique for the isotropic medium (e.g., Zhao et al. 1992). In other words, the influence of Vp anisotropy on the travel time and ray path is not fully considered. This effect could become serious, if strong anisotropy exists in the study region. Hence, it is necessary to develop a method of 3-D anisotropic ray tracing and use it to conduct anisotropic tomography. However, differences between phase velocity and group velocity in the anisotropic media make it difficult to perform 3-D anisotropic ray tracing. In general, phase velocity describes the advancing speed of wave front, whereas group velocity describes the propagating speed of energy along a ray trajectory. In isotropic media, when there is no dispersion, the phase and group velocity vectors have the same direction and amplitude, but they are usually different in anisotropic media. This difference leads to the complexity of 3-D anisotropic ray tracing. Gou et al. (2018) investigated this issue by modifying the isotropic 3-D ray-tracing scheme of Zhao et al. (1992) to trace P-wave rays in a 3-D weak RAN medium that contains complex-shaped velocity interfaces such as the Conrad and Moho discontinuities and the subducting slab boundary (Zhao et al. 2012).

Among the two velocities, group velocity (v) is directly related to the ray trajectory. Hence, group velocity (or slowness) is more concerned to perform the anisotropic ray tracing (Gou et al. 2018). In anisotropic media with hexagonal symmetry, the P-wave group velocity vector, phase velocity vector and the symmetry axis are located in the same plane. Taking advantage of this co-plane property, Gou et al. (2018) expressed the phase slowness vector $\vec{p}$ by combining the group slowness u (i.e., 1/v) and the unit vector of anisotropic symmetry axis ($\widehat{m}$) as follows:

$$\vec{p} = u\left( {\eta \hat{m} + \zeta \hat{u}} \right)$$

(17)

$$\eta = \left\{ {\begin{array}{*{20}l} 0 \hfill & {\Psi = 0,\frac{\pi }{2},\pi } \hfill \\ {\frac{\tan \Phi }{{\sin \Psi }}} \hfill & {\Psi \ne 0,\frac{\pi }{2},\pi } \hfill \\ \end{array} } \right.$$

(18)

$$\zeta = \left\{ {\begin{array}{*{20}l} 1 \hfill & {\Psi = 0,\frac{\pi }{2},\pi } \hfill \\ {1 - \frac{\tan \Phi }{{\tan \Psi }}} \hfill & {\Psi \ne 0,\frac{\pi }{2},\pi } \hfill \\ \end{array} } \right.$$

(19)

where $\widehat{u}$ is the unit vector of the group slowness, $\Phi$ is the angle between the phase and group slowness vectors, and $\Psi$ is the angle between the group slowness vector and the anisotropic symmetry axis. On the basis of these relations and considering $\nabla u=d\overrightarrow{p}/ds$, Gou et al. (2018) rewrote P-wave ray equation as:

$$- \frac{{d^{2} \vec{X}}}{{ds^{2} }} = \frac{1}{\zeta v}\left[ {\nabla v + \left( {v\frac{d\eta }{{ds}} - \eta \frac{dv}{{ds}}} \right)\hat{m} + \left( {v\frac{d\zeta }{{ds}} - \zeta \frac{dv}{{ds}}} \right)\frac{{d\vec{X}}}{ds}} \right]$$

(20)

where $\overrightarrow{X}$ is position vector at a point on the ray trajectory, and s is the ray length. Then they derived anisotropic schemes of the pseudo-bending technique (Um and Thurber 1987) and Snell’s law (Zhao et al. 1992) and used them iteratively to find the fastest ray path.

Applying the anisotropic ray-tracing code to calculate P-wave travel times and ray paths in the 3-D Vp RAN model of the Japan subduction zone (Liu and Zhao 2017), Gou et al. (2018) found that the travel-time difference caused by the Vp RAN is generally < 0.1 s and the ray-path difference is generally < 15 km (Fig. 16). Gou et al. (2018) also conducted tomographic inversions of the local and teleseismic travel-time data (Liu and Zhao 2017) using both the isotropic ray-tracing code (Zhao et al. 1992) and their anisotropic ray-tracing code. The two models of 3-D Vp RAN tomography they obtained are quite similar to each other (see Gou et al. 2018 for details).

Fig. 16

An isotropic ray (blue dashed line) and an anisotropic ray (red line) in a 3-D Vp radial anisotropy model in two vertical cross sections (a, b) and in a map view at 70 km depth (c). Locations of the vertical cross sections are shown in c. The red star denotes an earthquake hypocenter, and the green triangle denotes a seismic station. Red triangles denote active volcanoes in NE Japan (c) and those within a 50-km width of the profile AA’ (a). Black solid curves in a, b denote four velocity discontinuities including the Moho and the slab upper boundary. Note that in b the horizontal and vertical scales are quite different. After Gou et al. (2018)

The study of Gou et al. (2018) indicates that, at this stage, the use of isotropic ray-tracing codes is acceptable to conduct anisotropic tomography, taking into account the limited spatial resolution of the current models of anisotropic tomography (e.g., those shown in Figs. 2, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 and 15) and the influences of regularizations (e.g., smoothing and damping). However, for anisotropic tomography models with higher resolution in the near future, the 3-D anisotropic ray tracing will be necessary.

8 Discussion

The above-mentioned studies indicate that seismic anisotropy is a very useful and important physical parameter, because it can provide a wealth of information regarding dynamic processes in the Earth’s interior and a unique constraint on past and present deformations in the lithosphere and sub-lithospheric mantle. These studies have shown that seismic anisotropy and heterogeneity coexist in all parts of subduction zones and continental margins, which are related to complex mantle flows and dynamic processes associated with plate subductions and collisions. Lithospheric plates descend into the mantle and drive hot mantle upwelling, causing large earthquakes and vigorous magmatism along island arcs and continental margins. These processes induce different stress and strain regimes at different depths and generate varying seismic anisotropies. Hence, measurements of seismic anisotropy provide straightforward constraints on the subduction processes, in particular, mantle convection in subduction zones. So far, SWS measurements have been the most important tool to studying seismic anisotropy in all tectonic settings including subduction zones. However, a major drawback of the SWS method is that it has no depth resolution, so the obtained SWS results usually cannot provide a unique answer to key geodynamic issues.

Body-wave anisotropic tomography is a new but powerful tool for mapping 3-D variations of seismic anisotropy in the crust and mantle. It has excellent depth resolution, so it can resolve the geodynamic issues unanswered by the SWS studies. As introduced in the last sections, high-resolution images of seismic anisotropy tomography obtained in the past decade have provided essential and direct information on the 3-D anisotropic structure of subduction zones and continental regions. These new findings well constrain deformations of the crust and mantle lithosphere associated with the surface geology and tectonics, seismogenic processes in the subduction channel associated with mechanical coupling between the subducting slab and the overlying plate, interplate and intraplate magmatism associated with slab dehydration and subduction-driven corner flow in the mantle wedge, anisotropy in the subslab mantle related to entrained flow caused by subduction and 3-D flow caused by trench retreat, as well as olivine fabric transitions due to changes in water content, stress and temperature. The most significant scientific achievement of the anisotropic tomography approach is its constraints on anisotropy within the slab in subduction systems. Although the slab anisotropy was discussed briefly by some SWS studies, the new scientific insights into intraslab deformation and modification of preexisting slab fabric, among other processes, have been gained only from this aspect of tomographic studies. Therefore, the systematic studies of seismic anisotropy tomography in the past decade have revolutionized our understanding of the 3-D structure and dynamics of subduction zones and continental regions.

While the technical advances as mentioned in the last sections are significant, there are still several avenues to pursue in the improvement of seismic anisotropy tomography, which are summarized as follows.

• (1) S-wave arrival times are a very important source of data in seismology, because their quality and quantity have been comparable to those of P-wave arrival times in the past decade and from now, thanks to the deployment of three-component seismometers in modern seismic networks in various regions of the world. Hence, the S-wave data should be fully exploited in the study of anisotropic tomography. Liu and Zhao (2016b) used S-wave travel-time data to determine AAN tomography. It is necessary to study whether the S-wave data can also be used to determine RAN tomography. If yes, then the S-wave data can also be used to conduct a joint inversion for both AAN and RAN tomography, as was done by using P-wave travel-time data (Liu and Zhao 2017). Wang and Zhao (2021) used only P-wave arrivals to invert for 3-D tilting-axis anisotropy, it is necessary to investigate whether S-wave arrivals can be also used to do that.

• (2) As mentioned in Sect. 2, in almost all the practical studies of seismic anisotropy, only the $2\phi$ terms are concerned (see Eq. 1) because of the limited spatial coverage and resolution of the seismic data set used. When better data sets are available in future studies, the $4\phi$ terms should be also investigated.

• (3) In recent years, many researchers have attempted to conduct simultaneous inversion of different kinds of seismological data sets for tomographic images. For example, Liu and Zhao (2016c) determined Vs tomography of the Japan subduction zone by performing a joint inversion of local and teleseismic travel times and surface-wave data. Recently, 3-D Vs tomography down to 400 km depth beneath the Cape Verde hotspot was determined by jointly inverting teleseismic S-wave travel-time residuals, Rayleigh-wave phase velocities and receiver functions (Liu and Zhao 2021). However, the two studies determined only isotropic Vs tomography. It is necessary to extend their approaches to include anisotropy in the tomographic inversions. As mentioned above, Liu and Zhao (2016a) determined 3-D Vs AAN tomography of the upper mantle beneath the Japan Islands and the Japan Sea using surface-wave data (Fig. 13). It would be easy to combine the approaches of Liu and Zhao (2016a, 2016b) to determine Vs AAN tomography using both body-wave travel times and surface-wave data.

• (4) The 3-D anisotropic ray-tracing method of Gou et al. (2018) is applicable to P-wave rays in RAN media. It is necessary to extend their approach to both P and S wave rays in 3-D AAN media, as well as in the 3-D anisotropic model with tilting symmetry axes proposed by Wang and Zhao (2021). From now, with the increasing quality and quantity of seismological data recorded by dense seismic networks deployed in many regions, the spatial resolution of anisotropic tomography will become higher and higher, then the 3-D anisotropic ray-tracing codes will be very necessary (Gou et al. 2018).

• (5) In recent years, many 3-D P and S wave attenuation (Qp, Qs) models of the crust and upper mantle have been determined, particularly beneath the Japan Islands (e.g., Liu and Zhao 2015; Wang and Zhao 2019b; Hua et al. 2019). The spatial resolution of the Q tomographic models is comparable to that of velocity tomographic models. It would be very interesting to include anisotropy in the Q tomographic inversion to investigate whether Q anisotropy occurs, and if yes, whether it is consistent with the velocity anisotropy or not.

Seismic data quality and coverage control the first-order features of a tomographic model (Zhao 2015). The existing seismograph networks are still very sparse with an inhomogeneous distribution of seismic stations. For example, in East Asia and the western Pacific region, permanent seismic stations are mainly deployed on the Japan Islands, East China, and South Korea, whereas the other areas are still poorly instrumented, in particular, the oceanic areas. The sparse distribution of seismic stations is the major handicap to seismologists for studying the detailed 3-D structure and dynamics of this broad region. Gradual installation of permanent or portable seismic stations on the seafloors and in those less instrumented land areas such as the Tibetan Plateau will be the most important task for seismologists from now. For this purpose, international cooperation of seismologists in the related countries will be very necessary and important.

In future studies of seismic imaging, high-quality P and S wave arrivals and later phases, surface-wave and waveform data recorded by dense seismograph networks on land and seafloors should be collected and exploited. The different kinds of seismic data could be used simultaneously to conduct joint inversions. Thus, we will be able to obtain more detailed 3-D tomographic images of P and S wave velocity (Vp, Vs), Poisson’s ratio, attenuation (Qp, Qs), and anisotropy. The new results will further improve our understanding of the geofluid migration, seismogenic and magmatic processes in the crust and upper mantle, as well as the fate of subducting slabs and birth of mantle plumes. At the same time, we should also pay attention to updated findings in other branches of the Earth sciences so as to better interpret the seismological results.

9 Conclusions

Advances in the theory and applications of seismic anisotropy tomography in the past decade are reviewed. The obtained results have greatly advanced our understanding of the 3-D seismic structure and geodynamic processes in the crust and mantle of subduction zones and continental regions. The main findings of these anisotropic studies are summarized as follows.

• (1) The greatest scientific achievement of the anisotropic tomography approach is its unique constraints on the slab anisotropy. The subducting oceanic slabs generally exhibit trench-parallel fast-velocity directions (FVDs) of azimuthal anisotropy, which may reflect frozen-in lattice-preferred orientation of aligned anisotropic minerals and/or shape-preferred orientation due to transform faults formed at the mid-ocean ridge and normal faults produced at the outer-rise area near the trench axis.

• (2) Deformation of the subducting slabs plays an important role in both mantle flow and intraslab fabric. The widely observed trench-parallel anisotropy in the forearc by shear-wave splitting measurements and P-wave analyses results, at least partly, from the intraslab deformation during outer-rise yielding of the incoming oceanic plate. Toroidal mantle flow is widely observed in and around a slab edge or slab window (hole).

• (3) Trench-normal FVDs generally occur in the mantle wedge beneath the volcanic front and back-arc area, reflecting subduction-driven corner flows. FVDs are generally trench-normal in the mantle below the slab, reflecting asthenospheric shear deformation due to active subduction of the overlying slab.

• (4) Significant low-velocity (low-V) anomalies with negative radial anisotropy are revealed in the subslab mantle, which may reflect thermally buoyant materials ascending from the deep mantle. The subslab low-V anomalies may affect the nucleation and rupture processes of great megathrust earthquakes on the slab interface.

• (5) Ductile flow exhibiting significant seismic anisotropy is clearly revealed in the middle-lower crust beneath the Tibetan Plateau. Fluids in the ductile flow and the flow itself can affect seismic and volcanic activities, mountain buildings and other tectonic events.

• (6) Significant azimuthal and radial anisotropies occur in the big mantle wedge (BMW) under East Asia, which reflect convective flows associated with the deep subduction of the western Pacific plate and its stagnation in the mantle transition zone. The geodynamic processes in the BMW have caused the craton destruction, back-arc spreading, and intraplate volcanic and seismic activities in East Asia.