The 3D Seismic Azimuthal Anisotropies and Velocities in the Eastern Tibetan Plateau Extracted by an Azimuth‐Dependent Dispersion Curve Inversion Method

An azimuth‐dependent dispersion curve inversion (ADDCI) method is applied to Rayleigh waves to extract 3D velocity and azimuthal anisotropy. The synthetic tests show that the ADDCI method is able to extract azimuthal anisotropy at different depths. The errors of the fast propagation direction (FPD) and the magnitude of the anisotropy (MOA) are less than 10° and 1–2%, respectively. The 3D anisotropic model shows large variations in the FPDs and MOAs with depth and blocks; strong contrasts are observed across major faults, and the average MOA in the crust is approximately 3%. The FPDs are positively correlated with the GPS velocity direction and the strikes of regional faults in most of the blocks. The low‐velocity zones (LVZs) in middle to lower crust are widely observed in the Songpan‐Ganze terrane, the north Chuan‐Dian block, and surprisingly in the Huayingshan thrust and fold belt. The LVZs in middle crust are also positively correlated with a region of low velocities in the uppermost mantle. These observations may suggest that large‐scale deformation is coupled vertically from surface to uppermost mantle. Crust shortening by the pure shearing process, which involves the thrusting and folding of the upper crust and the lateral extrusion of blocks, may be the major mechanism causing the growth of the eastern Tibetan Plateau.

Tectonic history and structural complexities in the eastern Tibetan Plateau have been further complicated by a number of unresolved scientific issues. At the core of these issues is the geodynamic mechanism behind the structural deformation responsible for the growth of the eastern Tibetan Plateau. Two competing geodynamic models have been debated. The lower crust flow model (Clark & Royden, 2000;Klemperer, 2006;Royden, 1996;Royden et al., 1997;Royden et al., 2008) was proposed to account for the uplift of the eastern Tibetan Plateau largely by the vertical expansion of the weak lower crust without much crustal shortening. Tapponnier et al. (1982), Shaw (2009), andHubbard et al. (2010) proposed a pure shearing model to explain the growth. In this model, compression between the eastern Tibetan Plateau and the Yangtze Craton resulted in folding and thrusting of the detached upper crust and lateral expansion in the north-south direction. The differential expansions controlled by block geometry (i.e., more expansions to the north than the south) may have resulted in large-scale strike-slip faults, such as the Xianshuihe, Anninghe, and Xiaojiang faults (Figure 1).
Seismic velocity is one of the key parameters used to resolve these problems. For example, surface wave tomography based on noise data has become routine for regional structural studies since it was introduced into seismological studies in the early 2000s (e.g., Campillo & Paul, 2003;Shapiro & Campillo, 2004;Wapenarr, 2004). Particularly, for the Tibetan Plateau, there are a number of publications dedicated to ©2020 The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. seismic noise tomography (e.g., Bao et al., 2015;Li et al., 2009;Legendre et al., 2015;Yao et al., 2006Yao et al., , 2008. These works have made unique contributions to our understanding of the geodynamics of the eastern Tibetan Plateau. In addition to seismic velocity, seismic anisotropy is of particular value for evaluating the mechanism of tectonic deformation and the status of regional stresses because anisotropies are often attributed to the preferential alignment of planar or elongated structures, such as rock layers, cracks, or anisotropic minerals, such as olivine and mica, under regional stresses (Chang et al., 2015;Kong et al., 2018;Lev et al., 2006;León Soto et al., 2012;McMamara et al., 1994;Silver & Chan, 1991;Sol et al., 2007;Zhao et al., 2010;Wu et al., 2015). Anisotropy due to the preferential alignment of planar or elongated structures is called shape-preferred orientation (SPO), while anisotropy due to the preferential alignment of mineral lattice is called lattice-preferred orientation (LPO) (Stein & Wysession, 2003).
There are two types of anisotropy, that is, radial anisotropy and azimuthal anisotropy, that can be constrained by seismic data. Radial or polarization seismic anisotropy accounts for the different wave speeds of differently polarized shear waves. These waves can be estimated by inverting SH and SV wave velocities independently (Agius & Lebedev, 2014;Huang et al., 2010;Xie et al., 2013). Polarization anisotropy may also be measured from shear wave splitting of the SKS/SKKS phases (Chang et al., 2015;Hirn et al., 1995;Huang et al., 2015;Kong et al., 2018;León Soto et al., 2012;Lev et al., 2006;McMamara et al., 1994;Silver & Chan, 1991;Sol et al., 2007;Wang et al., 2008;Wu et al., 2015;Zhao et al., 2010), Pms phases (converted S waves from P waves at the Moho) from receiver functions (Chen et al., 2013;Levin & Park, 1998;Sun et al., 2015), and direct shear waves (Gao et al., 2013;Singh et al., 2015). Anisotropy from  Deng et al. (2003), and the white thick lines are the major block boundaries. The two thick blue lines mark the surface projection of the two profiles in Figure  SKS/SKKS splitting is often attributed to the mantle even though both the mantle and crust are sampled; Pms phases of receiver functions sample the crust only, and therefore, the resulting anisotropies reflect the accumulated effects of the whole crust beneath the station.
Azimuthal seismic anisotropies can be extracted from Pn travel times (e.g., Liang et al., 2004;McNamara et al., 1997). Pn waves travel mostly in the uppermost mantle, and therefore, the resulting azimuthal anisotropies are from the uppermost mantle. Pms phases from different directions can also be fitted with cosine functions to obtain the azimuthal anisotropies beneath stations (Sun et al., 2012;Sun et al., 2015;Zheng et al., 2018). Surface wave dispersions are the most widely used data to extract azimuthal anisotropies for different periods (Ceylan et al., 2012;Legendre et al., 2015;Su et al., 2008).
The aforementioned methods measure the accumulated azimuthal anisotropies in the crust or crust-mantle column with limited depth resolution. Several inversion algorithms have been used in the literature to extract 3D radial and azimuthal anisotropies. After obtaining phase or group velocity maps for each geographic point, Montagner and Nataf (1986) and Yao et al. (2010) inverted the dispersion curve to obtain the depth-dependent elastic moduli for anisotropic materials. By combining all 1D depth-varying velocities and anisotropies at different geographical grids, 3D velocity and 3D anisotropy models were constructed. Yao et al. (2010) applied the method to the southeastern region of the Tibetan Plateau, while Chen et al. (2016) applied the method to Yunnan Province of China (part of the southeast Tibetan Plateau). Xie et al. (2017) applied a similar technique but included tilting parameters to find the depth-dependent elastic moduli. Priestley et al. (2006) inverted the surface waveforms to obtain a path-specific 1D model at the first stage, and then the 1D models from all paths were combined to develop models including both 3D velocity and 3D azimuthal anisotropy. The study applied the method to long-period multimode surface waves to study east Asia. Three-dimensional anisotropies can also be obtained by inverting shear wave splitting times from multiple paths Zhang et al., 2007). Ozacar and Zandt (2010) used a waveform modeling method to obtain depth-dependent anisotropic structures beneath each station, and 3D anisotropies were obtained by assembling models for all stations.
Given these achievements, the debate on the deformation mechanism in this region still exists. For example, joint analysis of GPS and SKS/SKKS splitting results shows that deformations in the crust and mantle are largely coupled in the eastern Tibetan Plateau, which is inconsistent with the predictions of the lower crust flow model (Chang et al., 2015;León Soto et al., 2012;Sol et al., 2007;Wang et al., 2008;Yang et al., 2018). On the other hand, based on anisotropies from Pms phases, some studies (Kong et al., 2016;Sun et al., 2012;Sun et al., 2015;Zheng et al., 2018) have suggested the existence of a lower crust flow in the region. Chen et al. (2013) argued that a midlower crust flow may exist beneath the interior of the eastern Tibetan Plateau, but block extrusion resulting from pure shearing may occur in the Central Yunnan subblock. Agius and Lebedev (2017) argued that both mechanisms may occur to form the eastern Tibetan Plateau.
In this study, we propose an azimuth-dependent dispersion curve inversion method (ADDCI) to invert the group dispersion curves of Rayleigh waves to obtain a model with shear wave velocities varying with 3D space and azimuth. By analyzing the azimuthal variations, the anisotropic parameters can be obtained for each 3D grid. Compared with the existing methods, the ADDCI method is more straightforward to implement, and this is the first study to apply it to the eastern Tibetan Plateau. The obtained 3D velocity and anisotropy models are then used to constrain the geodynamic models in question in this region.

Geological Background
The complex tectonic history of the eastern Tibetan Plateau has resulted in a mosaic tectonic environment composed of multiple fragmented microblocks (Burchfiel et al., 1995;Molnar & Tapponier, 1975;Roger et al., 2010;Xu et al., 2008;Yin & Nie, 1993;Zhang et al., 2006). Figure 1 shows the major blocks that are essential to discuss the tectonic activities in this region. The eastern Tibetan Plateau consists of several small blocks from north to south, including the eastern Kunlun block (EKLB), the Songpan-Ganze terrane (SGT), the northern Chuan-Dian block (NCDB), the Qiantang block (QTB), the southern Chuan-Dian block (SCDB), and the southwestern Dian block (SWDB).
To the east of the Tibetan Plateau, the Yangtze Craton serves as a barrier for the eastward expansion of the Tibetan Plateau. The Yangtze Craton also consists of several small units from west to east, including the Sichuan basin (SB), the Huanyingshan thrust and fold belt (HTFB), the Yungui plateau, and the eastern Yangtze Craton.
Between the eastern Tibetan plateau and the Yangtze Craton is the eastern Tibet boundary belt (ETBB). It is a seismically active region with earthquakes (the Ms8.0 Wecnhuan, 2008;Ms7.0 Lushan, 2013;and the Ms7.0 Jiuzhaigou, 2017) that have occurred successively inside or near the belt in the last 11 years.
The Songpan-Ganze terrane is a Triassic orogenic belt with a widespread Triassic flysch complex, Mesozoic volcanism, and diachronous thrusting events Xu et al., 2008). The Longmenshan fault (LMSF) and Mingjiang fault (MJF) mark the eastern boundary of the Songpan-Ganze terrane. Inside the Songpan-Ganze terrane, approximately 200 km west of the LMSF, lies the SW-NE striking dextral Longriba fault (LRBF). This fault is believed to have been aided in the differential kinematic movements of the late Quaternary owing to the intense resistance of the Yangtze Craton. This fault may have absorbed part of the deformation and reduced the deformation on the LMSF Xu et al., 2008). As a result, a relatively small horizontal convergence rate of less than 3 mm/year across the LMSF is observed (Shen et al., 2005).
The northern and southern Chuan-Dian blocks are transition zones between the NW-SE directional motion in the north and southward motion to the east of the eastern Himalayan syntaxis (EHS) with respect to the Sichuan basin (Shen et al., 2005). The rotational movement inside these blocks is aided by large-scale strike-slip faults, such as the NW-SE striking Xianshuihe fault (XSHF) to the north, the N-S striking Anninghe fault (ANHF), the Zemuhe fault (ZMHF), and Xiaojiang fault (XJF) to the east and the Jinshajiang fault (JSJF) and Red River fault (RRF) zone to the west-southwest of the Chuan-Dian block. The Chuan-Dian block is divided into northern and southern parts by the Xiaojinhe fault (XJHF) (i.e., the northern and southern Chuan-Dian blocks). A number of studies have found that these two parts are structurally different (Bao et al., 2015;Chen et al., 2016;Yao et al., 2010;Zheng et al., 2018).
Located in the northwest region of the Yangtze Craton, the Sichuan basin is a foreland basin that has accumulated Paleozoic to Cenozoic sedimentary sequences with a bulk thickness of 6-12 km (Xu et al., 2008). Seismic studies (e.g., Liu et al., 2014;Wang et al., 2015;Zhu et al., 2013) have found that the thickness of the sediments may reach up to 14 km. To the southeast of the Sichuan basin, the Yungui plateau lineated by the Dalou fault (DLF) to the northwest and the Wuling fault (WLF) to the southeast is a large mountainous region with rugged terrain, including steep karst peaks and deep gorges. The Yungui Plateau has undergone two phases of uplifts, that is, the Yanshan orogeny in the Late Mesozoic and the Himalayan orogeny in the Cenozoic (Dong et al., 2018;Wan, 2013). Lying between the Yungui Plateau and the Sichuan basin is a transitional thrust and fold belt: the Huayingshan thrust and fold belt. This belt is characterized by lines of low-altitude mountains and valleys.

Data
The study region spans 21°N to37°N and 96°E to111°E (Figure 1), and it covers the major tectonic blocks in play in the eastern Tibetan Plateau and western Yangtze Craton. We collected ambient noise data from networks maintained by the Sichuan Seismological Bureau and the Yunnan Seismological Bureau, and data from other stations compiled by the China Earthquake Network Center and one subsidiarity of the China Earthquake Administration. In addition to the data from the national network, data from the Lms Gap experiment (the experiment monitored the seismic gap between the 2008 Wenchuan earthquake and the 2013 Lushan earthquake in southern Longmenshan; Liang et al., 2018) are also used. Overall, continuous seismic data recorded at 146 broadband stations in 2015 were used to extract Rayleigh waves.
We generally follow the steps outlined by Langston (2008, 2009) to extract surface waves from ambient noise and further to extract dispersion curves. The processes include decimating continuous data to 10 samples per second, removing trend and mean values, filtering in the period band centered at 5 to 50 s, converting waveforms to a 1-bit signal and cross-correlating. We then used phase-equalized frequency-time analysis (FTAN) algorithms to extract the group velocity curve for every ray.
After visually inspecting the cross-correlation functions and the velocity period spectrum by FTAN analysis, we obtained 5,500-6,700 rays for T = 5-10 s and approximately 6,000 rays for T = 10-30 s, and the number of rays decreased lineally from 4,000 to 6,000 for T = 30-40 s and from 2,000 to 4,000 for T = 40-50 s. The ray path maps of T = 20 s and T = 40 s are shown in the supporting information S1.

The Azimuth-Dependent Dispersion Curve Inversion Method
Compared to surface wave tomography without considering anisotropy, azimuthal anisotropy terms can be added in the same way as anisotropic Pn tomography terms (Liang et al., 2004). Without anisotropy terms, the travel-time residual (time difference between the synthetic and measured times) for each period can be written as where Δt ij is the travel-time residual for the great circle path that links the (i,j) th station pair, d ijk is the length of the (i,j) th ray in the kth cell (d ijk = 0 if the (i,j) th ray does not pass through the kth cell), and s k is the slowness (1/V group velocity ) of the kth cell.
Solving equation 1 yields the group velocity maps for different periods. By combining the group velocities for different periods, a new group velocity dispersion curve is obtained for each cell. A 1D velocity model (Vs velocity varying with depth) is then inverted using linear (Herrmann, 2013) or nonlinear (Liang & Langston, 2009) techniques. In the end, a 3D velocity model is obtained by assembling all 1D models.
The aforementioned technique has been successfully used to extract 3D velocity in many areas (e.g., Bao et al., 2015;Bensen et al., 2007;Liang & Langston, 2009;Li et al., 2009;Yang et al., 2012). However, in cases where strong azimuthal anisotropy exists, a modification to equation 1 is often applied, as shown in equation 2: In this equation, the travel-time perturbation includes the contribution from the background slowness perturbation (the first term on the right-hand side) and the contribution from azimuthal anisotropy (the second and third terms on the right-hand side). a k and b k are the anisotropy parameters in the kth cell, ∅ ijk is the azimuth of the segment of the (i,j) th ray inside cell k, and N is the total number of cells. Here, we assume the form of transverse anisotropy described by a sinusoidal 2∅ ijk azimuthal variation (Backus, 1965). Note that we use the azimuth of the ray inside of each cell (∅ ijk ) instead of the same azimuth (or back azimuth) for all cells along the ray because the azimuth (back azimuth) may change substantially along the ray.
When the anisotropy terms are relatively small compared to the isotropic term, then the velocity-azimuth curve is equivalent to a quasi-sine or cosine function: The isotropic velocity (v k ), the magnitude of anisotropy, that is, MOA (γ k ) and the fastest propagation direction, that is, FPD (∅ k ) are given by Other routine techniques, such as data weighting, 2D Laplacian regularization, and resolution and stability analysis, may also be applied. Implementations of these techniques are similar to the implementations in 10.1029/2019TC005747 Tectonics Liang et al. (2004) and Liang and Langston (2009), but some new aspects for this study are described below. The final equation system is solved using a classic least-square QR decomposition (LSQR) method (Paige & Saunders, 1982a, 1982b. With the anisotropic terms included, for every period, in addition to a group velocity map, an anisotropic map that presents the MOAs and FPDs can be obtained. As shown in Figure 2a, for the kth cell, we can obtain a group velocity and anisotropic parameters (a k and b k ). Both sets of parameters are combined to give a pseudo-velocity ellipse, as shown in Figure 2. Based on the pseudo-velocity ellipse of all periods, we can construct velocity period curves, that is, dispersion curves, for rays in any direction.  Step 1: Obtain the dispersion curve collections. See section 2.1 for details. 10.1029/2019TC005747

Tectonics
Step 2: Apply dispersion curve tomography to obtain 2D group velocity and anisotropy maps for different periods. In this study, the study region is divided into 0.5°× 0.5°cells. Due to the different ray numbers for different periods and to account for the finite frequency effects, higher smoothing parameters are chosen for longer periods. See section 3.1 for details. The group velocity and anisotropic maps are shown in Figure 6.
Step 3: For every cell, the group velocity and anisotropic parameters a k and b k for different periods are used to construct dispersion curves for azimuths between 0°and 170°at a constant interval of 10°. The idea is further illustrated in Figure 2.
Step 4: For an individual cell, the dispersion curves at different azimuths are inverted to obtain 1D velocity models (the velocities vary with depth) for all azimuths. Assuming that the study region is divided into N b cells and that the azimuthal step is 10°, the number of azimuthal bins is N a = 18 (i.e., 0°,10°, …,170°). The total number of 1D models is N b × N a .
Step 5: For an individual cell, the average of all N a velocities at each depth gives the isotropic velocity v 0 at that 3D point. If the angular interval is sufficiently smaller, the maximum and minimum velocities and the corresponding azimuths can be identified directly. The azimuth corresponding to the maximum velocity is the FPD, while the MOA is given by equation 7: If the angular interval is too large to meet the accuracy requirement, then the velocities for N a azimuths can be fitted to a sine or cosine function in the form of equation 3. Then, the MOAs and FPDs can be found using equations 5 and 6, respectively.
The last step is to combine the MOAs and the FPDs in all 2D cells and at all depths to form a 3D velocity and 3D anisotropic model.

Synthetic Testing
In this section, we use synthetic waveforms to test the usability and stability of the method to extract depth-dependent anisotropies. Waveforms are computed using the AXISEM package (Nissen-Meyer et al., 2007, which is based on the spectral element method. Figure 4a shows the top 60 km of the model with different velocities for different azimuths. From 60 km to the center of the earth, the model is identical to the isotropic PREM model (Dziewonski & Anderson, 1981). Different azimuthal anisotropies are added at different depths. For depth ranges of 0-20 km and 20-40 km, the anisotropic parameter (a, b) values are (0.006, 0.0) s/km and (0.002, 0.01) s/km, respectively. Figure 4b shows the velocity-azimuth curves at different depths. The red plus signs mark the maximum velocities at corresponding depths. Figure 4c shows the MOAs at different depths. Based on equations 5 and 6, the FPDs for the two depth ranges are approximately 90°and 129°, respectively, as illustrated by the red plus signs in Figure 4b. The MOAs for the two depth ranges are approximately 4% and 7%, respectively, as shown in Figure 4c. The waveforms for azimuths of 10°to 170°(at an increment of 20°) and constant epicenter distances of 50°are synthesized. The synthetic waveforms are plotted in the supporting information in Figure S2a. The FTAN of the Rayleigh wave for an azimuth of 10°is shown in Figure S2b.  Figure S2b.
The fitted anisotropic parameters, that is, (a, b) for different periods in this case, can be treated as anisotropy inverted for each cell described in section 3.2. Following Step 3 in Figure 3, azimuth-dependent dispersion curves are then constructed for azimuths from 0°to 170°at an increment of 10°. A total of 18 dispersion curves are then inverted to obtain 18 1D models, as shown in Figure 5. The inverted models are similar to the true models in Figure 4 in several aspects: the overall velocity-depth profiles are similar to each other for all azimuths (Figures 4a and 5a), and the inverted FPDs are approximately 90°and 130°for depths of 0-20 km and 25-40 km, respectively (Figure 5b), which are close to the true FPDs in Figure 4b. The inverted MOAs ( Figure 5c) at 0-20 km and 30-40 km are 3-4% and 5-6%, respectively. These numbers are slightly less than the true MOAs of 4% and 7% in the corresponding depth range (Figure 4c). The MOAs decrease to approximately 1% at a depth of 58 km. Although, as is often observed for most geophysical inversions, depth smearing indeed occurs for both FPDs and MOAs. The errors for the FPDs and MOAs are relatively smaller at 0-20 km than at deeper depths.
To summarize, the application of the ADDCI to synthetic data shows that this technique is able to resolve both FPDs and MOAs at different depths for the top 60 km. The errors of the FPD and MOA are approximately 10°and 1-2%, respectively.

The Stability and Resolution of the Surface Wave Tomography
As discussed in Liang and Langston (2009), the regularization parameters may be carefully chosen to guarantee stability of the inversion, but this works at the cost of the resolution. In this study, we use the jackknife method to test the stability of the inversion for different regularization parameters. The best parameters are the ones that yield relatively small errors while also guaranteeing a relatively good resolution. With the regularization parameters for shear wave velocities and anisotropies given in Table 1, the errors of the velocities, FPDs, and MOAs are less than 0.1 km/s, 20°, and 1%, respectively, for most of the study region ( Figure S3). More constraints are placed on the longer periods because compared to the shorter periods, the rays are less dense; this also works when accounting for the finite frequency effects of longer periods. Compared to the parameters used for velocity, higher constraints are placed on the anisotropic parameters to ensure that the errors of the FPDs and MOAs for most of the study region are small enough for meaningful interpretation. Tectonics Figure 6 shows the checkerboard test results for a resolution of 2°× 2°for T = 10, 20, 30, and 40 s. For clarity, the anisotropic checkerboard tests are also plotted separately in Figure S4. The general patterns of both velocity and anisotropy are resolvable for most of the region for T = 10-30 s. The smearing effects are obvious for period T = 40 s, and the resolution is particularly poor in the southwestern part of the study region. With more extensive tests for different pattern sizes, the best resolution for the central part of the study region can be as optimal as 1°× 1°. Based on these tests, the region with a resolution better than 3°× 3°for all periods is analyzed in the 4D inversion shown in later sections. Figure 7 shows the group velocity and anisotropy maps for periods of 10, 20, 30, and 40 s, respectively. The sensitivity kernels are shown in the supporting information in Figure S5.

Group Velocity and Anisotropy Maps for Different Periods
The group velocity of the Rayleigh waves at T = 10 s is sensitive to the Vs structure in the uppermost 10-15 km. Due to the thick sediments overlaying the crystallized basement beneath the Sichuan basin, the group velocity at T = 10 s is predominantly lower than 2.7 km/s, which is a sharp contrast from its adjacent regions, where the group velocity is primarily larger than 2.9 km/s. The velocity in the vicinity of the Sichuan Table 1 The

10.1029/2019TC005747
Tectonics basin is relatively higher than that in regions farther away. The MOAs at T = 10 s are the largest among the four periods plotted in Figure 6. This makes sense because the upper 20 km is more susceptible to elastic deformation under regional stress, which forms anisotropic structures. The FPDs of the anisotropies change dramatically across the region. The FPDs in the north and south Chuan-Dian block are mostly aligned in the NNW-SSE directions and are consistent with the strikes of the large-scale strike-slip faults defining the two blocks. The large-scale faults appear to separate regions with dramatically different FPDs, for example, the LRBF and the LMSF.
For T = 20 s, the velocities in the Sichuan basin, the Huayingshan thrust and fold belt (between the HYSF and DLF), the southwestern Dian block and the southern part of the southern Chuan-Dian block are less than 2.9 km/s. The MOAs are smaller than T = 10 s in most of the study region. The FPDs show significant block characteristics. In the north Chuan-Dian block, the FPDs are nicely aligned with the strike-slip bounding faults, such as the XSHF and the RRF. This is consistent with the movement directions of materials in this region. Similar to those at T = 10 s, relatively large anisotropies are observed to the southeast of the Sichuan basin. For T = 40 s, the Rayleigh waves are mostly affected by the midlower crust and the uppermost mantle. As a response to the strong basement of the Yangtze Craton, the group velocities of the Sichuan basin, the western Yangtze Craton, and the southwestern Dian block are dominantly higher than their adjacent regions. The FPDs in the north Chuan-Dian block and south Chuan-Dian block are similar to those at T = 10-30 s. An interesting feature is that the Sichuan basin is isolated from the rest of the western Yangtze Craton by a linear low-velocity belt between the HYSF and the DLF, that is, the Huayingshan thrust and fold belt.

Dispersion Curve Inversion for One 2D Cell
Dispersion curves are constructed for different azimuths for every 2D cell (Step 3 in Figure 3). A linear dispersion curve inversion (Herrmann, 2013) is conducted to obtain a 1D velocity model for every azimuth. The nonlinear inversion method may also be used, but it can be very time consuming because the number of dispersion curves is N a times the number of isotropic inversions. Figure 8 shows the dispersion curves before inversion and the seismic models after inversion in the 2D cell centered at (103.25°, 30.25°). This approach is comparable to the synthetic test in section 2.3. Figure 8a shows the velocity variation as a function of azimuth for different periods before 4D dispersion inversion. For each period, the velocity-azimuth curve is computed using equation 3. The MOAs for different periods shown in Figure 8b are computed using equations 5 or 7.
After inversion, Figure 8c shows the 18 1D models (velocity varies as a function of depth) for 18 azimuths from 0°to 170°. Figure 8d shows the results in Figure 8c in a different way: the velocities vary with azimuth for different depths; Figure 8e shows the MOAs for different depths. The MOAs are computed using equation 7. In Figures 8a and  8d, the numbers on each function mark the maximum and minimum velocities, and they are generally separated by 90°, representing the peak and trough of the cos(2θ) function.
The angular interval used in this study is 10°, and we believe this value is accurate enough for the discussion of anisotropy in this study. For every depth in this 2D cell, the FPD can be readily determined from the velocity-azimuth function by searching for the maximum velocity. The isotropic velocity v mean is defined as the average velocity of all azimuths. The MOA is obtained using equation 7 (Stein & Wysession, 2003).

10.1029/2019TC005747
Tectonics Surprisingly, without any smoothing schema applied after inversion, the independent inversions of the 18 dispersion curves in Figure 8a yield a smooth-varying cosine-shaped velocity curve for all depths (Figure 8d). The magnitudes of the velocity-azimuth curves vary and are equivalent to half the MOA, as shown in Figure 8e.
The 18 1D models shown in Figure 8c with different colors have similar curvatures. A low-velocity zone is observed at approximately 35 km for all azimuths. The velocities vary slightly at depths of 0-5 km, vary slightly more at depths of 20-30 km, and vary much more at depths deeper than 50 km.
Overall, the amplitudes of the velocity-azimuth functions are larger in the mantle (>50 km) than in the crust. The strong anisotropies near the surface may result from the preferred alignment of the highly deformed structures in the Songpan-Ganze terrane near the LMSF region. The strong anisotropies at depths of 20-35 km may be attributed to the low-velocity zone, which is a typical feature in the Songpan-Ganze terrane. The increasingly larger MOA at the base of the crust and uppermost mantle may be related to the lattice-preferred orientation of olivine in the uppermost mantle.

3D Velocity and Anisotropy Model
Combining the velocity and anisotropy at all depths for all 2D cells, a 3D velocity model and an anisotropy model can be produced. At the 2 km depth (Figure 9a), the Sichuan basin is characterized by extremely low shear wave velocities (<2.7 km/s) due to the thick sediments in the basin. The velocities in the surrounding regions, such as the Songpan-Ganze terrane, the north Chuan-Dian block, the south Chuan-Dian block, and the Yungui Plateau, are approximately 3.1-3.5 km/s. The FPDs of anisotropies in the Songpan-Ganze terrane align dominantly perpendicular to the Longmenshan fault system but parallel to the strike-slip faults inside the terrane. In the north Chuan-Dian block, the FPDs are mostly in the E-W direction, but they are quite complex in the south Chuan-Dian block.
The most interesting feature is that the FPDs of shear waves around the Sichuan basin are dominantly perpendicular to the edges of the basin. For example, to the northwest of the Sichuan basin between the LMSF and the LRBF and to the southeast of the Sichuan basin, large anisotropies with FPDs in the NW-SE direction are observed and are perpendicular to the Sichuan basin edges. The MOAs are relatively smaller in the south Chuan-Dian block than in the north Chuan-Dian block.
At 10 km depth (Figure 9b), the velocity contrast between the basin and surrounding regions is not as significant as that at 2 km depth. The velocities in the northern Sichuan basin are relatively lower (<3.2 km/s) than those in the southern Sichuan basin (3.3 km/s) due to the thicker sediments to the northern part of the basin. The LRBF appears to separate the two different regimes, with FPDs in the E-W direction to the NW and FPDs in the NW-SE direction to the SE direction. The differences in the FPDs may reflect different deformation mechanisms on both sides of the fault. The Huayingshan thrust and fold belt, the region between the HYSF and the DLF, is associated with a low-velocity belt with large range, parallel anisotropies (MOA~10%), which are almost perpendicular to those at a depth of 2 km. To the southeast, the Yungui Plateau and western Yangtze Craton are dominated by high velocities (>3.4 km/s) with fault-parallel anisotropies. Except for the Huayingshan thrust and fold belt and the Yungui plateau, the MOAs at this depth are generally smaller than those at 2 km, and the FPDs inside the Sichuan basin are quite different between the two depths. This may suggest that the sedimentary layers and crystalline basements have different responses to regional stresses.
At 20-32 km depths (Figures 9c and 9d), the relative velocities between the Sichuan basin and the Songpan-Ganze terrane are reversed. The velocities of the Sichuan basin (3.6-3.7 km/s) are primarily higher than those of the surrounding region (<3.5 km/s), which indicates that the basement of the Sichuan basin is relatively stronger than that of the Songpan-Ganze terrane. The basement of the Yungui Plateau is also extremely strong. Similar to the 10 km depth, the region between the HYSF and the DLF, that is, the Huayingshan thrust and fold belt, is associated with a very low-velocity belt (<3.3 km/s at a 20 km depth) with FPDs parallel to the strike of the belt. The general pattern and velocity range at a depth of 30 km are consistent with the 35-km velocity map by Xie et al. (2013).
The MOAs at 32 km are generally smaller than those at 2 and 10 km depths for most regions. The FPDs for the region northwest of the LRBF seem to change gradually from the E-W direction at 10 km depth to the SW-NE direction (

10.1029/2019TC005747
Tectonics suggest that the direction of maximum stress changes with depth. For the Chuan-Dian block, the FPDs are generally parallel to the faults and the movement direction of the materials.
At 40-50 km depth (Figures 9e and 9f), the Sichuan basin and the whole Yangtze Craton are composed of the uppermost mantle, while the Songpan-Ganze terrane and the north Chuan-Dian block are still within the lower crust. Although the anisotropies are strong beneath both the Sichuan basin and the Songpan-Ganze terrane, they are relatively weak beneath the north Chuan-Dian block and the south Chuan-Dian block. There appears to be a continuous low-velocity zone (LVZ) extending from the Songpan-Ganze terrane to the north Chuan-Dian block. The velocity of the south Chuan-Dian block is relatively higher than that of the north Chuan-Dian block.

The Magnitudes of Anisotropies as a Function of Depth
As shown in Figure 8d, the FPDs annotated by plus signs (the maximum velocity on a cosine-shaped curve) vary significantly with depth. The FPDs are approximately 50°near the surface. At a depth of~6-7 km, as the velocity gradient ( Figure 8c) and the MOAs (Figure 8e) start to make sharp changes, the FPDs increase from 50°to 120°. The rapid changes in the FPD and the MOA may indicate that a physical boundary separating layers with different compositions may exist at this depth. A prominent FPD contrast is also observed at approximately 40 km, where a 40°contrast is observed, and the MOAs also start to increase beyond this depth. We speculate that this depth might be where the Moho begins. In this sense, the anisotropic features may serve as a powerful tool to identify the Moho and other physical boundaries, and it may also help to classify rock types at different depths, although further analyses are needed to establish meaningful correlations.
The 3D anisotropies obtained in this study not only show lateral variations but also show significant depth-dependent variations. To analyze how the MOAs vary with depth in the entire region, Figure 10 shows the statistics of the MOAs at different depth ranges. For a depth range, an effective velocity is calculated for each azimuth to form a new cosine function to extract the FPDs and MOAs using equation 7. The result is equivalent to the weighted average of all curves within the depth range. For example, the weighted average of the top 3 curves in Figure 8d yield the velocity-azimuth curve for the depth range from 0 to 10 km for this cell.
The MOAs are especially large at depths of 20-30 km (Figure 10c), with a median MOA of approximately 6%. This depth is characterized by a widespread distribution of an LVZ. This makes sense because under high temperature and high pressure, the materials with low velocity in the midlower crust (deeper than 20 km for most of the region) are more ductile than the materials in the upper crust. Therefore, this region is more susceptible to tectonic stresses to form either SPO or LPO (Stein & Wysession, 2003). This is also the case for the 2D cell shown in Figures 8c and 8e, although the depth and extent of the LVZ may vary for individual cells.
Relatively large MOAs (~5.5%) are also observed at 0-10 km depths. This may suggest that, without the confining stresses exerted at certain depths, the materials near the surface can be easily deformed to yield the SPO type of anisotropy.
The increase in the MOA across the Moho at~40-50 km depths is also significant, with increases from approximately 4.3% at the lower crust (40-50 km beneath the eastern Tibetan Plateau) to 5.6% at the uppermost mantle (Figure 10f). This sharp contrast may be attributed to the strong anisotropy of olivine in the uppermost mantle. The minimum and maximum shear wave velocities of olivine are approximately 4.42 and 5.53 km/s, respectively. The MOA is approximately 22% (Stein & Wysession, 2003). Yao et al. (2010) found that the azimuthal anisotropies are approximately 3-4% at the lower crust, which is close to the 4.3% result found at 30-50 km (Figures 10d and 10e) in this study. The median MOA of the 0-50 km depths is approximately 3.1%, as shown in Figure 10g. A majority of the MOAs are between 2% and 4%. This number is consistent with Su et al. (2008), who found that the azimuthal anisotropies in the crust of eastern Tibetan are larger than 2%. However, this number is much smaller than the MOAs at any single depth range shown in Figures 8a-8e. The velocity-azimuth curves in Figure 8d may help to explain the collectively lower average MOAs for the depth range of 0-50 km. Because the FPDs vary significantly with depth, the cosine functions at different depths are out of phases to yield small stacked amplitudes.
Therefore, the average MOAs of the 0-50 km depths are collectively smaller. This could be the major reason the measured crustal azimuthal anisotropies based on Pms splitting are small in many studies (Chen et al., 2013;Sun et al., 2012Chang et al., 2015. Ozacar and Zandt (2010) found that the anisotropies along the INDEPT-III profile are approximately 10%, 5%, and 18% at depth ranges of 0-3, 10-16, and 24-32 km, respectively. Although the numbers are much higher in their study than in our study, it does show that the highly deformed layer near the surface and the weak middle crust are more anisotropic than those at other depths, as shown in our result.

Tectonics
In summary, significant anisotropic changes are observed at different depths. How the depth-dependent anisotropies can be integrated with stresses and temperatures to study the compositions and geodynamics at various depths will need further investigation.
However, we also should realize that the absolute MOAs may be affected by the relative regularization between the anisotropic and isotropic terms in equation 2. Extensive tests have shown that the relative regularization may influence the absolute MOAs but have little influence on the FPDs and the relative MOAs in the lateral and vertical directions. Figure 11 shows the weighted average velocities and anisotropies of two depth ranges. The depth ranges of 0-20 and 20-50 km roughly represent the upper crust and the midlower crust, respectively. The average velocities of the two depth ranges are dramatically different. The Sichuan basin shows a strong midlower crust with a low-velocity upper crust due to its thick sediments. Nonetheless, the Yungui Plateau and western Yangtze Craton are the strongest at both depth ranges compared to other regions. The existence of the Emeishan large igneous province in the north Chuan-Dian block makes the south Chuan-Dian block appear to be stronger than the north Chuan-Dian block. The Huayingshan thrust and fold belt has a lower velocity upper crust and midlower crust compared to the Sichuan basin and the Yungui Plateau (on both sides).

The Average Crustal Anisotropies Compared to the GPS Velocity Field
However, the FPDs of the anisotropies at the two depth ranges are largely consistent. For example, the FPDs are mostly parallel to the local faults inside the north Chuan-Dian block, south Chuan-Dian block, and the western Yangtze Craton (including the Sichuan basin, the Huayingshan thrust and fold belt, and the Yungui Plateau) for the two depth ranges. In particular, strong NE-SW FPDs are observed in the Huayingshan thrust and fold belt, and they are parallel to the strikes of the faults inside this narrow low-velocity belt. To the southeast of the Huayingshan thrust and fold belt, the strikes of the faults on the Yungui Plateau are nearly perpendicular to those inside the Sichuan basin and the Huayingshan thrust and fold belt, and the FPDs in this unit also take sharp turns to align with the strikes of the faults at both depth ranges.
The alignment of the FPDs with major faults may imply that the anisotropies mostly result from the preferential alignment of the planar fault fabrics created by the tectonic stresses. The overall consistency of the FPDs at the two depths may imply that the deformations are coupled between the two depths, at least in the north Chuan-Dian block, the south Chuan-Dian block, the Huayingshan thrust and fold belt, and the Yungui Plateau.
The major inconsistency between the two depths occurs at the western Songpan-Ganze terrane (west of the LRBF), where the FPDs in the upper crust are mostly parallel with the strikes of the faults, while the FPDs in the midlower crust are mostly in the SW-NE direction, which is more consistent with the direction of ductile deformation. That is, the anisotropies in the upper crust of the western Songpan-Ganze terrane are more SPO originated, while the anisotropies in the midlower crust are more LPO originated (Stein & Wysession, 2003) due to the preferred orientation of minerals such as mica in the midlower crust.
Anisotropy can be measured from shear wave splitting of the Pms, SKS, and SKKS phases. The first two measure anisotropy in the crust, while the last one measures anisotropy resulting from the whole mantle and the crust. Chen et al. (2013), Sun et al. (2015), and Zheng et al. (2018) measured the azimuthal anisotropies of the crust based on Pms phases from different azimuths, and their results are compared with the weighted average FPDs of the top 50 km obtained in this study ( Figure 12). Figure 12 shows the overlapping of three datasets: the average anisotropy of the top 50 km determined from this study (red bars); the azimuthal anisotropy determined from Pms phases obtained by Chen et al. (2013), Sun et al. (2015), and Zheng et al. (2018) (cyan bars); and the GPS velocity field (Gan et al., 2007) with respect to the Eurasia-fixed reference frame (blue bars). The angular differences between the FPDs determined in this study and that based on Pms phases (Figure 12b), or the GPS directions (Figure 12c) are dominantly less than 30°.
The median delay time based on the Pms phases shown in Figure 12 is 0.38 s ( Figure S6). Assuming an average crust shear wave velocity of Vs = 3.5 km/s and an average crust thickness of H = 50 km, this delay time corresponds to a MOA = 2.7%, which is close to the median MOA = 3.1% of the top 50 km, as shown in Figure 10g. Figure 12. (a) Comparison of the average anisotropy of the top 50 km determined from this study (red bars), the anisotropies determined by shear wave splitting analysis (cyan bars, based on Pms phase) and the GPS velocity with respect to the Eurasia-fixed reference frame (Gan et al., 2007); (b) angular differences between the FPDs determined by this study and by shear wave splitting analysis; (c) angular differences between the FPDs determined by this study and the GPS directions; (d) scales for different datasets.

10.1029/2019TC005747
Tectonics It is very interesting to find that in some cases, the directions of the FPDs vary significantly across major faults. The FPDs in the western Songpan-Ganze terrane are dominantly NWW-SEE and are consistent with the GPS velocity directions in this region, while FPDs are mostly perpendicular to the LMSF between the LRBF and the LMSF and are also consistent with the GPS velocity field; the FPDs change across the LMSF from NW-SE in the Songpan-Ganze terrane to NE-SW beneath the LMSF and in the Sichuan basin. The XJHF divides the Chuan-Dian block into the northern part (north Chuan-Dian block) and southern part (SCDB). Although the FPDs change gradually from NW-SE in the north Chuan-Dian block to the N-S direction in the south Chuan-Dian block across the XJHF fault, they are mostly parallel to the strike-slip faults in both blocks.
Inside the Songpan-Ganze terrane, north Chuan-Dian block, Qiang Tang block, and the Yungui Plateau (western Yangtze Craton), the FPDs from this study are well aligned with the GPS velocity directions with only limited exceptions. The major disagreement between these two datasets occurs in the Sichuan basin and the Huayingshan thrust and fold belt, where the two directions are nearly perpendicular to each other.
The FPDs from the Pms studies are more scattered than those obtained in this study. Nonetheless, it is clear to find that the FPDs from the two methods are largely consistent in the Songpan-Ganze terrane, the north Chuan-Dian block, and the southern Sichuan basin. Beneath the LMSF and western Sichuan basin, the FPDs from this study and the Pms studies are primarily parallel to the fault system.
Inside the south Chuan-Dian block, the dominant N-S FPDs from this study are systematically different from the GPS velocity at an angle of 30-45°. Considering that FPDs from this study are very well aligned with the large-scale N-S striking faults inside the south Chuan-Dian block, while the GPS directions are parallel to the block boundary of the RRF, the systematic differences may reveal that regional displacements are largely accommodated by the RRF, while the N-S striking faults inside the south Chuan-Dian block are formed by the NW-SE striking regional stresses. This is consistent with the results of the uniaxial compressional experiment: the angle between the generated faults and the maximum compressional stresses is approximately 45°.
The different tectonic processes may lead to different angular differences between the FPDs and the GPS velocity. Both the GPS velocity and the FPDs of anisotropies are related to regional tectonic stresses. In some regions, the flow of materials results in the alignment of planar or linear structures, such as rock layers, fault fabrics, and elongated rocks in the direction of the flow. The surface wave traveling in this direction preferentially travels in the fastest layer/column. In this case, the FPDs are parallel or subparallel to the GPS velocity directions. This scenario may apply to the north Chuan-Dian block and the Songpan-Ganze terrane. In other cases, especially in a thrust and fold belt, where the movement of materials is blocked by strong barriers, the planar or elongated structures most likely align parallel or subparallel to the tectonic boundary, such as the LMSF, but perpendicular to the GPS velocity directions in the far field. The second scenario may explain the large angles between the GPS velocity and the FPDs, such as the NE-SW-aligning FPDs beneath the LMSF and the Huayingshan thrust and fold belt. Figure 13 shows two 3D illustrations of the inverted velocity model. Figure 13a shows the cross section passing the Wenchuan earthquake epicenter and perpendicular to the LMSF (AA' in Figure 1). The Sichuan basin is characterized by thick sediments on top of the strong midlower crust, with a Moho depth of approximately 40 km. A prominent LVZ (with a velocity ranging from 3.2 to 3.5 km/s) is observed at approximately 20-30 km beneath the Songpan-Ganze terrane. Surprisingly, a low-velocity crust is also present to the southeast of the Sichuan basin between the HYSF and the DLF, that is, the Huayingshan thrust and fold belt. Beneath this zone lies a relatively low-velocity upper mantle (<8.0 km/s, e.g., Liang et al., 2004). The Yungui Plateau to the SE of the DLF has a relatively strong crust and mantle.

The Low-Velocity Zone and the Deformation Mechanism
The LMSF is located above the westward extension of the strong basement of the Sichuan basin, and this is probably the tectonic environments responsible for the large earthquake like Wenchuan earthquake . Finer seismic imaging by Liu et al. (2018) reveals that the midlower crust of the seismic gap between the Wenchuan and Lushan earthquakes is occupied by LVZs extending from the interior of the Songpan-Ganze terrane. The LVZ beneath the seismic gap is thought to result from partial melting at the 10.1029/2019TC005747 Tectonics midlower crust (Nelson et al., 1996), and it heats up the faults above and makes them aseismic in nature . Figure 13b shows the cross section cutting across the LRBF, the XSHF, and the XJHF (BB' in Figure 1). The cross section shows the vertical profile of the Songpan-Ganze terrane, the north Chuan-Dian block, and the south Chuan-Dian block from NW to SE. Across the XSHF, the LVZ is also prominent beneath the north Chuan-Dian block. However, to the SE of the XJHF, the LVZ is less prominent and is only confined to the lower crust.

10.1029/2019TC005747
Tectonics To illuminate the LVZ distribution with depth, Figure 14 shows the model beneath 14 km with velocities less than 3.5 km/s and larger than 4.0 km/s. The LVZs are well defined by major tectonic boundaries. The LVZs are obviously absent from the Sichuan basin and the Yungui Plateau. The upper bound of the LVZs in the Songpan-Ganze terrane and the north Chuan-Dian block is shallower than 14 km. Consistent with Figure 13b, the LVZ beneath the south Chuan-Dian block is much thinner and is confined only to the lower crust. Similar to the Songpan-Ganze terrane and the north Chuan-Dian block, a thick and narrow LVZ bounded by the HYSF and the DLF is present beneath the Huayingshan thrust and fold belt between the strong basin to the west and the Yungui Plateau to the east. The LVZ in the middle crust has often been interpreted as partial melting (e.g., Hacker et al., 2014;Liang et al., 2018;Nelson et al., 1996;Xie et al., 2017).
The northwest part of the south Chuan-Dian block is occupied by the strong Emeishan large igneous province, which may have played a role similar to the Sichuan basin in stopping the expansion of the LVZ toward the south Chuan-Dian block. Huang et al. (2010) found strong negative anisotropies associated with the Emeishan large igneous province. Xie et al. (2013) also found that relatively weak anisotropies exist beneath the Emeishan large igneous province in the lower crust compared to the northern and southern Chuan-Dian blocks.
Another point to make is that the locations of the LVZs appear to correlate positively with the topography and this is one of major arguments that supports the lower crust flow model (Royden et al., 1997, Clark & Royden, 2000. For example, the XJHF separates the higher topography of the north Chuan-Dian block from the lower topography of the south Chuan-Dian block; this coincides with the much thicker LVZ beneath the north Chuan-Dian block compared to the south Chuan-Dian block. The correlation is even more prominent in the Songpan-Ganze terrane compared to the Sichuan basin across the LMSF. The LVZ beneath the Huayingshan thrust and fold belt is also correlated with the relatively higher topography across the surface compared to the topography of the Sichuan basin. However, this correlation does not hold for the Yungui Plateau, where the strong crust is associated with relatively higher topography compared to the Huayingshan thrust and fold belt. The locations of the LVZs in the midlower crust are also strongly correlated with the low-velocity layer in the uppermost mantle constrained by Pn wave travel time (Liang et al., 2004;Liang & Song, 2006). The Sichuan basin is found to be strong both in the midlower crust and uppermost mantle, while the region around the basin is found to be relatively weak in the midlower crust and the uppermost mantle (Figures 13 and 14 in this study, Figure 2 in Liang & Song, 2006). The strong correlation between the crust and mantle at a large scale may suggest that the deformations are vertically coupled from the midlower crust to the uppermost mantle, at least at a regional scale. However, Legendre et al. (2015) found that the deformation in the crust is decoupled from the rest of the lithosphere, but their conclusion was based only on the azimuthal anisotropies at T = 20, 40, and 80 s extracted from a limited number of rays without knowing the depth-dependent anisotropies. One may argue that the low velocity in crust and mantle may result from different mechanisms. For example, the upwelling of the mantle may cause a low-velocity belt in the uppermost mantle. If this is the case, a widespread distribution of volcanoes around the Sichuan basin would be expected, and this is not the case. Sol et al. (2007) also pointed out the consistency among the crustal structures, surface strain, and mantle anisotropy. They suggest that deformation in the lithosphere is mechanically coupled across the crust-mantle interface and that the lower crust is sufficiently strong in transmitting stress.
As discussed in section 4.1, the overall consistent FPDs at depths of 0-20 and 20-50 km suggest that the deformation in the upper crust is largely coupled with the deformation in the midlower crust. Integrated with the 3D velocity, we can make the assumption that the deformation in the upper part of the lithosphere, including the whole crust and the uppermost mantle, is largely coupled for a majority of the

10.1029/2019TC005747
Tectonics eastern Tibetan Plateau and the western Yangtze Craton, with the only exception occurring in the western Songpan-Ganze terrane. This result is consistent with the SKS splitting analyses in the eastern Tibetan Plateau by Sol et al. (2007), Wang et al. (2008), León Soto et al. (2012), Chang et al. (2015), and Yang et al. (2018).
For the lower crust flow model (Clark & Royden, 2000;Klemperer, 2006;Royden, 1996;Royden et al., 1997;Royden et al., 2008), the vertical expansion of the LVZ in the midlower crust plays a leading role in creating high topography. This model can explain the positive correlation between the LVZs and high topographies. If the ductile midlower crust flows in the NW-SE direction from the Songpan-Ganze terrane and the north Chuan-Dian block to the south Chuan-Dian block and another channel flows in the NE-SW along the Huayingshan thrust and fold belt, the FPDs in the middle to lower crust (Figure 11b) in those blocks are generally consistent with the flow directions. This is also one of the major arguments used by Kong et al. (2016), Sun et al. (2012), Sun et al. (2015), Legendre et al. (2015), and Zheng et al. (2018) to conclude the existence of a lower crust flow. However, a highly ductile and flowing channel would result in decoupled deformation of the upper crust and lower crust with the uppermost mantle, which is inconsistent with coupled deformation in the lithosphere, as discussed previously.
From the aspect of stress distribution, the vertical expansion yielded from the lower crust flow model would require high-angle maximum compressional stresses. This is, however, not supported by the stress map inverted from focal mechanisms by Yang et al. (2017Yang et al. ( , 2018.
For the pure shearing model (Hubbard et al., 2010;Hubbard & Shaw, 2009;Tapponnier et al., 1982), the LVZ in the middle crust serves as a lubricant to detach the upper crust from the lower crust and to form thrust and fold belts in the upper crust. Laterally, this results in block movement along large strike-slip faults. This model is also consistent with the positive correlation between the LVZs and high topographies and also explain the positive correlation between the FPDs and the strikes of the faults in the north and south Chuan-Dian blocks. Furthermore, the pure shearing model is also consistent with coupled deformations in the whole lithosphere, as discussed previously. The low-angle P axes (the direction of the maximum compressional stress) computed by Yang et al. (2017Yang et al. ( , 2018 are also consistent with the deformations accommodated by thrusting, folding, and strike slipping. Based on the above assessments, we conclude that the pure shearing model that results in vertical thrusting and folding and lateral strike slipping is the dominant deformation in large scales. However, we cannot rule out the existence of a midlower crust flow at local scales or in the interior of the Tibetan Plateau (e.g., west of the LRBF). Figure 15 shows a schematic model that reflects the general features of the velocity and anisotropy obtained by this study. The pure shearing model is applied to explain coupling between the lower crust and uppermost mantle and the lateral displacement of the tectonic blocks in the southeast Tibetan Plateau. The compressional stress incurred by the northward and eastward subduction of the India plate is transferred to the Songpan-Ganze terrane and the north Chuan-Dian block. Under increased pressure and temperature, the midlower crust is weakened; possibly catalyzed by the decreased solidus of rocks due to the existence of water, partial melting may take place and result in LVZs in the midlower crust beneath the Songpan-Ganze terrane and the north Chuan-Dian block. The LVZs serve as a lubricant to detach the upper crust from the lower crust to form a thrust and fold belt at the western edge of the Sichuan basin along the LMSF. The Emeishan large igneous province played a leading role in limiting the southeastward expansion of the LVZ into the south Chuan-Dian block. The basalt widely outcropped in the Emeishan may be too strong to be weakened, and the layers beneath the basalt layer

10.1029/2019TC005747
Tectonics can still undergo partial melting and therefore result in a less prominent LVZ at the base of the crust in the south Chuan-Dian block (Figures 13b and 14).
The existence of the LVZ belt beneath the Huayingshan thrust and fold belt is surprising. It may indicate that a strong unit, such as the Yangtze Craton, may still be broken up under compressional stress resulting from the India-Eurasia collision or earlier orogenic event.

Conclusions
In this study, we used a simple and straightforward azimuth-dependent dispersion curve inversion method, that is, the ADDCI method, to invert a 3D velocity and anisotropy mode in the eastern Tibetan Plateau and the western Yangtze Craton. Significant variations are observed both laterally and vertically. The sharp contrast of MOAs and FPDs at depths may help to identify compositional boundaries.
The FPDs and the MOAs present strong block variations. The positive correlation between the FPDs and the strikes of major faults may indicate that the preferred orientation of the planar or elongated fault fabrics is the major mechanism to form anisotropies in those blocks.
The MOAs are found to be particularly large at the depths of 0-10, 20-30, and 50-60 km. They are most likely due to the strong deformation near surface, low-velocity zone at middle crust and the preferential alignment of the Olivine in the uppermost mantle, respectively.
The average MOAs at 0-50 km depths are between 2% and 4%. This is much smaller than the MOAs at any single depth ranges because anisotropies at different depths may cancel out. Therefore, the anisotropies measured from Pms phases are useful to study the averaged anisotropic features of the crust, but they are often too simple to reveal the vertical variations of compositions, stresses, and temperatures.
The LVZs are geographically correlated with the LVZs in the uppermost mantle and also positively correlated with the high topography. The FPDs are overall consistent between the depth ranges of 0-20 and 20-50 km. These observations may suggest that the large-scale deformation is coupled vertically from the surface to the uppermost mantle. Integrating with other evidences in stress analysis and seismic imaging, we suggest that the pure shearing process is the major mechanism to account for the growth of the Eastern Tibetan plateau.