The lithospheric-scale deformation in NE Tibet from joint inversion of receiver function and surface wave dispersion

As a region far away from the Indian-Asian collision front, Northeastern (NE) Tibet has attracted great attention due to its implications for the plateau growth and the devastating earthquakes. Several geodynamic models have been proposed for the deformation mechanism of NE Tibet, but it remains controversial. Here we introduce the crustal Vp/Vs ratios from the H-k stacking to the joint inversion of receiver function and surface wave dispersion, which is based on the Neighborhood Algorithm. The crustal and uppermost mantle velocities are obtained beneath 203 stations in NE Tibet. The results show (1) deep Moho in Songpan-Ganzi and Qiangtang terranes, and distinct Moho depth between northwest and southeast Qilian; (2) positive correlation between the average velocity in the crust and in the uppermost mantle, suggesting a coupled lithospheric-scale deformation; (3) diverse high upper mantle velocities distributing in northern part of the study region, which may indicate the existing of local lithosphere or a non-uniform southward subduction of Asian lithosphere. Article history: Received 18 January 2018 Revised 20 December 2018 Accepted 18 January 2019


INTroducTIoN
Several geodynamic models have been proposed to explain the uplift and outgrowth of the Tibetan Plateau (TP), such as rigid block model, continuum model, and crustal flow model.Block model suggests that lithospheric-scale faults exhibit as boundaries separating the continent interior into a limited number of tectonic blocks, and the deformation takes place mainly along the strike-slip faults (Tapponnier et al. 2001).Continuum model allows stress to transfer within and across blocks and the deformation is broadly distributed (e.g., Neil and Houseman 1997;Zhang et al. 2004).In addition, the crustal flow model, in which the mid-lower crust of the plateau interior moves eastward and is diverted to the northeast and southeast around the Sichuan basin, is commonly used to explain the topography gradient around the TP (Royden et al. 1997;Clark and Royden 2000).
Due to its location and devastating earthquakes, Northeastern (NE) Tibet has long been an ideal site to study the deformation and uplift of the TP.Widespread mid-crustal low-velocity zones (LVZs) (Bao et al. 2013;Li et al. 2014) and low resistivity zones (Le Pape et al. 2012) were observed in the Songpan-Ganzi terrane (SGT) that penetrate across Kunlun fault (KF), consistent with the crustal flow model.However, isolated LVZs were observed under Qilian from ambient noise (Bao et al. 2013;Li et al. 2014), deep seismic sounding (Zhang et al. 2011(Zhang et al. , 2013)), and deep seismic reflection results (Gao et al. 2013;Guo et al. 2015), which favor continuous shortening as a deformation mechanism.Furthermore, observed Moho offsets (Zhu and Helmberger 1998;Vergne et al. 2002;Zhang et al. 2011;Ye et al. 2015) and the correlation between GPS and SKS splitting (Wang et al. 2008;León Soto et al. 2012;Chang et al. 2017) support a coupled lithospheric scale deformation (Tapponnier et al. 2001;Yin et al. 2008;Ye et al. 2015).The images from joint inversion of receiver function (RF) and surface wave dispersion argued that the rigid block model and crustal flow model could be concomitant (Liu et al. 2014).In summary, though abundant geophysical investigations have been conducted in this region, the deformation mechanism in this region is still controversial, which is the key to understand the uplift and growth of the TP.
Since more geophysical data have become available, joint inversions have been commonly used to resolve the ambiguities inherent in using individual methods separately.A common method is the joint inversion of surface wave dispersions and RFs (Julià et al. 2000).However, the traditional joint inversion usually sets the crustal Vp/Vs to the global average of 1.75, which would bias the crustal structure (Li et al. 2017a).Moreover, the results from the traditional linear inversion may also be affected by the starting model.In order to better illuminate the lithospheric structure in NE Tibet, here we utilize the data recorded by 230 stations and apply the nonlinear joint inversion of RFs and surface wave dispersions with independently constrained Vp/Vs ratios (from H-kappa stacking).The approach allows us to resolve lithospheric S velocity, surface sedimentary layers, and Moho depth.Our results suggest a coupled lithospheric-scale deformation in NE Tibet, which could provide a possible common process of the plateau growth.

dATA
Because we introduce the specific crustal Vp/Vs on the joint inversion of RF, surface wave dispersion, so the data used in this study include three types, i.e., RFs, dispersion curves of surface wave group and phase velocities, and Vp/ Vs ratio from H-k stacking.
The RF data came from 230 broadband stations archived from 2007 -2011, including 90 stations from Incorporated Research Institutions for Seismology (IRIS), 118 stations from China Earthquake Administration (CEA) recorded during 2007-2009(Data Management Centre of China National Seismic Network 2007;Shen et al. 2008;Zheng et al. 2010), and 22 stations from Institute of Geology and Geophysics, Chinese Academy of Sciences (IGGCAS) recorded between November 2010 and June 2011 (Xu et al. 2014;Wu et al. 2015) (Fig. 1).To calculate the RFs, we selected teleseismic P waveforms from earthquakes with magnitude Mw ≥ 5.5 and epicentral distance range from 30 -90°.Before using the time-domain iterative deconvolution method (Ligorría and Ammon 1999) to obtain the RFs, we applied the bandpass filter of 0.1 ~ 1.0 Hz to the raw data and Gaussian filter with a width of 2.5 to reduce high-frequency noise.Afterwards, we performed a quality control process based on the fitting ratio, amplitude, and harmonic stripping (Shen et al. 2013;Deng et al. 2015) on the RFs of all the stations (~23000).Finally, we obtained 8354 good observed RFs over 203 stations.
The dispersion data were extracted from the recent surface wave tomography of China (Bao et al. 2015), who used both ambient noise and earthquake data beneath more than 1300 stations, including 864 CEA permanent stations from 2008 -2011, 401 temporary PASCCAL stations, and 51 permanent stations from the IRIS Data Management Center.More than 700000 dispersion curves were measured to generate group and phase velocity maps at periods of 10 -70 s (group velocity maps extend to 140 s).In this study, according to the inverse-weighted distance algorithm, we extracted group and phase velocity dispersion curves at each of our stations from the dispersion maps of 10 -70 s, which are better constrained with both the ambient noise and earthquake data.
80% of Vp/Vs data in this study were from Zheng et al. (2016), who used 847 teleseismic earthquakes and obtained Vp/Vs ratios beneath 143 stations by using the H-k stacking method.In addition, for IGGCAS stations, we used the Vp/Vs ratios from Xu et al. (2014).For the remaining stations, we obtained the Vp/Vs ratios by applying the H-k stacking method (Zhu and Kanamori 2000).Two examples of the estimation are shown in Figs.2b and c.Although the measurements have uncertainties, they generally have a consistent pattern among stations (Fig. 2), which shows the relatively higher Vp/Vs in Lhasa terrane (LT), Qiangtang terrane (QT), and SGT, but lower in Qilian.

METhod
The teleseismic P-wave RF is sensitive to shear velocity contrast and depth-velocity product, instead of velocity alone, while the surface wave dispersion is sensitive to vertical shear-velocity averages but insensitive to velocity discontinuities.Thus, joint inversion of RF and surface wave dispersion has been widely used to reduce parameter ambiguity in the inversion (Julià et al. 2000).However, the joint inversion is greatly influenced by the Vp/Vs ratio in the crust, which is often unconstrained and normally set to the global average (1.75).The wrong Vp/Vs ratio would bias the crustal structure (Li et al. 2017a).
Here, we use the crustal Vp/Vs ratio beneath every station obtained from H-k stacking and perform the joint inversion of receiver function and dispersion, which is based on the fast global search method, the Neighborhood Algorithm (Sambridge 1999;Xu et al. 2013;Li et al. 2017a).The model parameterization is similar as Li et al. (2017a).The crust and the mantle are parameterized as cubic splines, with fixed numbers of nodes of 12 and 9, respectively.The depth of the last (bottom) node in the mantle is fixed at 150 km.Specifically, we introduced a separate sedimentary layer at the surface, which has a great influence on the waveform of the RF, particularly the beginning portion.We used a linear gradient for the sedimentary layer, but allowed the thickness to vary.We fixed the Vp/Vs ratio for the sedimentary layer at 2.0 as it trades off with the sedimentary layer thickness (not of interest in this study), the Vp/Vs ratio for the rest of the crust beneath the station from H-k stacking, and the ratio for the mantle at the global average of 1.795.The velocities below 95 km from Bao et al. (2015) are also used as reference in the joint inversion, where the resolution of dispersion curves degrades.The surface wave inversion alone is applied first to get the initial model for the joint inversion, and 200200 total models are searched to get the final best fitting model.Li et al. (2017a) has done several synthetic tests to illustrate the robustness of this method.Here we perform two synthetic tests to emphasize that the joint inversion could recover the correlated or anti-correlated velocity structure between crust and mantle (Fig. 3).For test 1, the input model has 2 km sediment, three layers in the crust, which includes a low velocity zone in the middle crust over a half-space of low velocity in the mantle.For test 2, the input model is the same as test 1 in the crust, but has high velocity in the mantle.Even though the surface wave inversion could recover the dispersion data (group and phase), the recovered velocity model is smoother than the input model.However, the joint inversion can recover well both input models, including the input Vs values and discontinuities.The predicted RF and dispersion curves for the inverted final model also agree well with the input ("observed") data.These results suggest that the correlation between crustal and uppermost mantle Vs observed in our data is robust, rather than biased from the trade-off between crustal and mantle structures.

Joint Inversion using real data
Two examples using real data (for station JCDC and AMU) are shown in Fig. 4. Obviously, the Vs model from joint inversion reveals more detailed features in the crust and has higher resolution in the mantle, compared with the Vs model from surface wave inversion alone.Note the inclusion of the sedimentary layer improves the fit to the direct P arrival near zero time in the RF, in comparison with the previous study (Xu et al. 2013).Because of the depth trade-off in surface wave dispersion, the improved resolution for the crust will also improve the resolution for the mantle, which is mostly constrained by the dispersion data and the background reference model.That is another advantage of this study compared with that from Zheng et al. (2016).Specifically, the Moho depth beneath JCDC is ~46 km, and there is a distinct high velocity layer in the crust.Comparatively, the Moho depth beneath AMU is ~45 km, and the velocity is smooth in the crust but high in the uppermost mantle.

Velocity Maps at different depths
The sedimentary layer has strong effect on the RFs and may bias the final results.The inclusion of searching the thickness and velocities of this layer not only make this method able to recover well the real sedimentary layer, but also prevent its contamination on the robustness of the crustal LVZ distribution and deep crustal structure (Li et al. 2017a).The velocity map at 2 km shows distinct low velocity beneath basins, such as the Qaidam basin, Gonghe basin, and Ordos basin (Fig. 5a).These features highlight the effectiveness of introducing sedimentary layer in the inversion.
The velocity map at 20 km (Fig. 5b) shows the uppermiddle crustal characteristics.SGT and QT are typified by low velocities, which may have penetrated across the KF.This result is consistent with previous studies (Le Pape et al. 2012;Jiang et al. 2014;Li et al. 2014Li et al. , 2017b)).We can also find isolated low velocity in northwestern Qilian, which is considered as the initial stage of deformation (Bao et al. 2013).
The velocity map at 60 km (Fig. 5c) indicates Moho depth variation.In general, the lower velocity indicates a deeper Moho.For example, SGT, QT, and LT terranes present low velocity at this depth, which is consistent with the deeper Moho.The same fea ture could also be observed in the northwestern Qilian.The Qaidam basin has a shallower Moho, sandwiched between Qilian and SGT.
The velocity at 80 km (Fig. 5d) illustrates the property of uppermost mantle.Most of the SGT and QT exhibit low velocity, which reveals a close correlation between the thick (> 60 km) crust and the presence of an upper-mantle lowvelocity zone beneath the QT and SGT.

Moho depth
With the inclusion of constrained Vp/Vs in the joint inversion, the Moho depth is more robust than the traditional joint inversion (Li et al. 2017a).The deepest Moho (~75 km) in this study area is observed beneath QT (Fig. 6a).The Moho lies at the depth of 60 -70 km in the SGT.In the Qaidam basin, however, the Moho becomes shallower to ~50 -55 km depth.The NW Qilian has a relatively deeper Moho (~55 -65 km depth) than that of the SE Qilian, Ordos basin and West Qinling (~40 -55 km depth).Compared with the Moho depth from H-k stacking (Xu et al. 2014;Zheng et al. 2016), our Moho depth is more robust because the multiple phases are not generally clear and the sediment would also bias the Moho depth (Yeck et al. 2013).
Figure 6b presents the standard deviation of the 20000 models that have smaller misfits with the observations in our search-based inversion scheme.The standard deviation  is quite small compared to the Moho variations, indicating that our inversion results are robust.

does the Asian lithosphere subduct beneath NE TP?
It is well accepted that the Indian lithosphere has subducted beneath southern TP, even though the frontier is still debated.However, whether the Asian lithosphere has subducted beneath TP remains controversial.Meyer et al. (1998) suggested the south-directed subduction in the NE Tibet based on surface geology and geomorphology.Tapponnier et al. (2001) proposed a model of oblique southward lithosphere subducting beneath NE TP.Regional seismic-reflection profiles and strata data (Yin et al. 2008), as well as receiver function studies of local temporary seismic arrays have reported evidence for possible southward subducting of Asian lithosphere beneath central and northern Tibet (Kind et al. 2002;Zhao et al. 2010Zhao et al. , 2011;;Ye et al. 2015).However, a seismic tomography study (Liang et al. 2012) observed a low velocity upper mantle beneath northern Tibet, which argues against subducting of the Asian lithosphere.Wei et al. (2017) used the Rayleigh wave tomography and further proposed a limited southward underthrusting of the Asian lithosphere.Furthermore, the recent S-RF images are inconsistent with the model of Asian mantle lithosphere subducting beneath the Tibet Plateau (Shen et al. 2015).
In order to better illustrate the velocity structures at different terranes, we slice four profiles from SW to NE cutting across the region.Profile (a) spans from the boundary between LT and QT to the North China Craton (NCC).The Moho depth decreases gradually from ~70 km in QT to ~65 km in SGT, to ~ 60 km in Kunlun-Qaidam terrane (KQT), and then increases to ~65 km in Qilian and finally decreases to ~50 km in NCC.Prominent low velocities are present in QT and SGT, and have penetrated across the KF.Qilian also has low velocity in the middle crust.Profile (b) extends from KQT to the NCC.The Moho is deeper in Qilian, but shallower in KQT and NCC.A low velocity zone is also present in the middle crust of Qilian.The profiles (c) and (d) have relatively flat Moho.The low velocities still exist in the crust beneath SGT and KQT along Profile (c), but not Qilian.Nevertheless, there is no apparent low velocity existing along profile (d).In the northern part of these four slices, the uppermost mantle exhibits high velocity.However, the high velocities in profiles (a) and (b) of the upper mantle do not connect with the southern part.
In summary, the velocity structure is heterogeneous in different tectonic units in this region.Even though the connected high velocity of the upper mantle in the northern part of profiles (c) and (d) seems to suggest a southward subduction of the Asian lithosphere, we do not see this fea-ture in profiles (a) and (b).Here we propose two alternative explanations.First, the Asian lithosphere does subduct beneath TP, but the subduction is not uniform.Second, there is no Asian southward subduction, and the high velocity just reflects the variation of local lithospheric structure.The heterogeneity may explain the different models proposed in the previous studies (discussed above).

The coupled deformation in NE Tibet
Moho offsets across major faults are widely reported in this region (Zhu and Helmberger 1998;Vergne et al. 2002;Zhang et al. 2011;Tian et al. 2014;Ye et al. 2015;Deng et al. 2018).The correlation between GPS measurements, crustal anisotropy and the mantle deformation field measured by shear-wave splitting suggests that the deformation of lithosphere is vertically coherent (Wang et al. 2008(Wang et al. , 2016;;Chang et al. 2017).Azimuthal anisotropy derived from multimode Rayleigh wave tomography implied coupled deformation throughout the whole lithosphere (Pandey et al. 2015).
From Fig. 7, we have found another interesting feature that the low velocity in the crust corresponds to the low velocity in the mantle.In order to better confirm this hypothesis, we calculate the average velocity in the crust (from 4 km below sediment to Moho) and uppermost mantle (top 20 km). Figure 8 shows the distribution of crustal and uppermost mantle velocities and there was a significantly positive correlation between them.In other words, the stronger crust corresponds to the stronger mantle, which is consistent with the coupled lithospheric-scale deformation in this region (Deng et al. 2018;Xu et al. 2018).It should be noted that, the Vp/Vs have uncertainties from the H-k stacking method, while the different Vp/Vs would mainly have trade-off with the Moho depth but not the velocity (Li et al. 2017a;Deng et al. 2018), which herein does not change our final conclusion.

coNclusIoN
In an effort to understand the deformation in NE Tibet, we perform a nonlinear joint inversion of receiver function and surface wave dispersion with specified Vp/Vs beneath 230 stations.The crustal and uppermost mantle velocities are obtained beneath the 203 stations.The main conclusions from the detailed velocity structure include: (1) The crustal thickness is heterogeneous in different tectonic units.
(2) Positive correlation exists between crustal velocity and uppermost mantle velocity, which suggests a coupled lithospheric-scale deformation in this region.(3) Two alternative explanations are proposed for the diverse high velocities in the upper mantle distributing in northern part of the region, which may indicate a non-uniform southward subduction of the Asian lithosphere or just reflect local lithosphere variations.
Fig. 2. The Vp/Vs ratio used in this study.The tectonic units are the same as Fig. 1.The circles with black border indicate the results from the H-k stacking method in this study.Others refer to Xu et al. (2014) and Zheng et al. (2016).Two stations for (b) and (c) are labeled by red characters in (a).The average crustal Vp/Vs ratio and Moho depth of the two stations are marked with white circles in (b) and (c).

Fig. 3 .Fig. 4 .
Fig. 3.The synthetic tests of the joint inversion of RF and dispersion with specified Vp/Vs (1.75 in the crust).(a) The input model has a 2-km sedimentary layer, three layers in the crust (18 km upper crust, 20 km middle crust, and 20 km lower crust) and one layer in the mantle.(b) The input model is the same as (a) except the high velocity in the mantle.The left panel: the best fitting model of the joint inversion (JI) (blue) and the best fitting model of the surface wave inversion alone (SI) (red).The arrow marks the Moho depth from the JI.The right topo panel: the observed RF (ORF, black envelope) and the predicated (calculated) RF from the best-fitting model of the JI (PRF) (blue); the right bottom panel: the observed Rayleigh wave phase and group dispersion (OD) curves (black), the predicated Rayleigh wave phase and group dispersion (PDJ) curves from the JI (blue), and the predicated Rayleigh wave phase and group dispersion (PDS) curves from the SI (red).The inverted sedimentary and crustal thicknesses from joint inversion are similar with the input model, as well the velocity value.

Fig. 7 .
Fig. 7.The topography and velocity structure beneath the four profiles as marked in Fig. 6.The tectonic units are the same as Fig. 1.