Distribution of the crustal low velocity zones beneath the central and northeastern Tibetan Plateau: Insights from joint analysis of receiver functions and surface wave dispersion observations

Accepted Manuscript Distribution of the crustal low velocity zones beneath the central and northeastern Tibetan Plateau: insights from joint analysis of receiver functions and surface wave dispersion observations Mengkui Li, Tengfei Wu, Xu Lin, Yujin Hua PII: DOI: Reference:

S0031-9201(18)30296-6 https://doi.org/10.1016/j.pepi.2018.12.003 PEPI 6216

To appear in:

Physics of the Earth and Planetary Interiors

Received Date: Revised Date: Accepted Date:

17 October 2018 20 December 2018 20 December 2018

Please cite this article as: Li, M., Wu, T., Lin, X., Hua, Y., Distribution of the crustal low velocity zones beneath the central and northeastern Tibetan Plateau: insights from joint analysis of receiver functions and surface wave dispersion observations, Physics of the Earth and Planetary Interiors (2018), doi: https://doi.org/10.1016/j.pepi. 2018.12.003

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Distribution of the crustal low velocity zones beneath the central and northeastern Tibetan Plateau: insights from joint analysis of receiver functions and surface wave dispersion observations

Mengkui Lia,b,*,Tengfei Wuc,*, Xu Lind, Yujin Huaa a

Department of Geophysics, School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China;

b

Key Laboratory of Geospace Environment and Geodesy of Ministry of Education, Wuhan University, Wuhan

430079, China; c

State Key Laboratory of Water Resources and Hydropower Engineering Science, Institute of Engineering Risk

and Disaster Prevention, Wuhan University, Wuhan 430072, China; d

College of Earth Sciences, Chengdu University of Technology, Chengdu 610059, China

*Corresponding author: Mengkui Li ([email protected]); Tengfei Wu ([email protected])

Abstract A crustal and uppermost mantle shear-wave velocity model beneath the central and northeastern Tibetan Plateau is obtained from the joint analysis of Rayleigh wave phase dispersions and teleseismic P-wave receiver functions at 188 broadband seismic stations deployed in this region. This newly presented model reveals wide-spread low velocity zones in the middle crust, and these zones are mainly confined to the region south of the Kunlun Fault. The velocities in the low velocity zones are extremely low (<~3.3 km/s), suggesting the presence of actively melting or previously melted materials in the crust, and further implying the possibly large contributions of the crustal channel flow model to the deformation in the Tibetan Plateau. Additionally, a continuous low velocity anomaly in the uppermost mantle beneath the eastern Kunlun Fault, the Songpan-Ganzi terrane and part of the Qiangtang terrane is observed, which may suggest an induced local mantle upwelling in this region.

Keywords: Tibetan Plateau; Joint inversion; Receiver functions; Surface wave dispersions; Partial melt

1. Introduction The ongoing collision beginning ~50 Ma between the Indian and Eurasian plates resulted in the uplift of the Tibetan Plateau (TP), a highly deformed zone with thickened crustal thickness and elevated topography (Rowley, 1996; Yin and Harrison, 2000). Deformation has persisted since the early stages of collision. At present, GPS observations reveal both inter-plate convergence between the two plates and intra-plate deformation within TP (Bilham et al., 1997; Gan et al., 2007; Pan and Shen, 2017; Pan et al., 2018). For decades, the evolution of TP has been a fascinating topic for geoscientists (e.g., Ceylan et al., 2012; Jiang et al., 2014; Kind et al., 1999; Li et al., 2008, 2016, 2017; Liang et al., 2012; Nunn et al., 2014; Royden, 1997; Royden et al., 2008; Tapponnier, 2001; Tilmann et al., 2003; Wu et al., 2016; Xu et al., 2018; Yang et al., 2012; Yin and Harrison, 2000; Zheng et al., 2016), especially the mechanism responsible for the Tibetan growth and uplift. Several models have been proposed to explain the crustal deformation of TP, such as the lateral extrusion of rigid block model (Avouac and Tapponnier, 1993; Tapponnier, 2001), uniform shortening and thickening of the lithosphere within TP (England and Houseman, 1986; Yang and Liu, 2013), and the crustal channel flow model (Clark and Royden, 2000; Royden, 1997; Royden et al., 2008). A recent study in the eastern TP from the joint inversion of surface-wave dispersion (SWD) data and receiver functions (RF) reveals that both the rigid block model and the crustal channel flow model may contribute significantly to the plateau expansion (Liu et al., 2014). If the crustal channel flow model contributes to the deformation of TP, it would predict crustal partial melts, which are directly indicated by the crustal low velocity zones (LVZs) in seismic studies. Therefore, to better understand the deformation mechanism in TP, additional information on the lithospheric structure is still needed,

encompassing the entire or local parts of TP. Numerous previous studies based on different tomographic methods have been undertaken to investigate the lithospheric structures during the past decades (e.g., Bao et al., 2015b; Ceylan et al., 2012; Hung et al., 2010; Jiang et al., 2016; Kind et al., 1999; Liang et al., 2012; Nunn et al., 2014; Pandey et al., 2015; Priestley et al., 2006; Tilmann et al., 2003; Wei et al., 2017; Xie et al., 2013; Yang et al., 2012; Yao et al., 2008; Zhang et al., 2011). These studies have revealed some common first-order features of the crustal and uppermost mantle structures beneath TP, such as significantly heterogeneous lithospheric structure and widely distributed crustal LVZs. Most of these studies involve only one geophysical observation, e.g., RFs (e.g., Kind et al., 1999; Murodov et al., 2018; Owens and Zandt, 1997; Xu et al., 2007; Zhu et al., 1995), P-wave traveltime (e.g., Hung et al., 2010; Li et al., 2008; Liang et al., 2012) and surface-wave group/phase velocities (e.g., Ceylan et al., 2012; Chen et al., 2016; Jiang et al., 2014; Wei et al., 2017; Yang et al., 2012; Yao et al., 2008; Zheng et al., 2016), which may suffer the intrinsic limitations of the geophysical techniques, such as the low vertical resolution of P-wave tomography, insensitivity to the absolute shear-wave speed of RF and insensitivity to the conversion interfaces of SWD. Joint analysis of multiple datasets can provide tighter constraints on velocity structures than inverting each data set independently (e.g., Bodin et al., 2012; Julià et al., 2000; Li et al., 2018a; Sosa et al., 2014). In this study, we constructed a 3-D shear-wave velocity model of the crust and uppermost mantle beneath the central and northeastern TP through the joint inversion of Rayleigh-wave phase velocity dispersion from 10 to 70 s (Bao et al., 2015a) and teleseismic P-wave RFs using seismic data recorded at the broadband seismic stations deployed in this region. Our results provide new

insights into the detailed deep structure under the central and northeastern TP, especially the distribution of crustal LVZs, which are possibly related to the presence of actively melting or previously melted materials in the crust and crustal channel flow.

2 Data and Method 2.1 Data The seismic data used for calculating radial P-wave RFs were downloaded from the IRIS Data Management Center, which were recorded by the temporal and permanent digital broadband stations deployed in the central and northeastern TP (Fig. 1). The total number of the stations is 188, including: 4 from the China National Seismic Network (Zheng et al., 2010), 51 from Namcha Barwa (Sol et al., 2007), 89 from INDEPTH II/ASCENT (Sandvol et al., 2008), 19 from INDEPTH IV (UK) (Zhao et al., 2008), 14 from the Broadband Seismic Experiment in Various Regions in China (Central China), and 11 from the Northeast Tibet Plateau Seismic Experiment (NETS) (Li et al., 2013). We show the network code associated to each network in Table 1. The locations of all the seismic stations are shown in Table S1. We selected events with moment magnitudes between 5.5 and 7.5, and within epicentral distances between 30° and 90°. The Rayleigh-wave phase velocity dispersion data were extracted from a recent surface wave tomography of mainland China (Bao et al., 2015a), in which the precise Rayleigh-wave phase velocity up to 70 s were determined using cross-correlation of ambient noise and earthquake data. In this study, the phase velocity dispersion data ranging from 10 to 70 s at each station were interpolated from the surface wave tomography model and used as the SWD dataset for the joint inversion.

2.2 Retrieval of radial P-wave RFs The processing strategies adopted by Li et al. (2016, 2018b) were used to process the waveforms. We first removed the trend and the mean from the raw seismograms, resampled the seismograms to 10 sps (samples per second) and performed a causal bandpass Butterworth filter with corner values of 0.05 and 2.0 Hz to remove noises with low frequencies. Subsequently, we rotated them from north and east components to corresponding radial and transverse components, manually picked the P-wave arrival times, and cut the time window to 20 s before and 120 s after P-wave arrival times. We then applied a time-domain iterative deconvolution technique (Ligorría and Ammon, 1999), which is based on the least square minimization of the difference between the observed and predicted horizontal component waveforms, to extract radial RFs from observed seismograms. A Gaussian filter width of 1.5 was selected to attenuate signals with frequencies greater than ~0.7 Hz. We discarded the RFs with waveform fits less than 85%. This procedure is hereinafter referred as the preliminary quality control of RFs. Fig. 2 shows the distribution of selected events after preliminary quality control. Finally, for each station, we adopted the following strategies to obtain the stacked average RF (Fontaine et al., 2015, 2013): (1) Bin the RFs into four quadrants with back-azimuths between N0° and N90°, N90° and N180°, N180° and N270°, and N270° and N360°, and then select the quadrant that has most RFs. (2) Compute the median pmedian of the ray parameters of all seismic events in the selected quadrant in step (1). (3) Select events with a ray parameter= pmedian ±0.004 (s/km). According to the experience of

the previous studies (Fontaine et al., 2015, 2013; Lamarque et al., 2015), this range of parameter is narrow enough to avoid adopting moveout corrections to the RFs. (4) Select the most mutually coherent RFs from the selected RFs in the previous step. Since most stations have available RFs (after the previous step) less than 50, we directly selected the coherent RFs by visual inspection instead of using the Cross-Correlation Matrix method (Tkalčić et al., 2006, 2011; Fontaine et al., 2015, 2013). (5) Stack the most mutually coherent RFs selected in the previous step to obtain the average RF, which will be used for the joint inversion. We show the individual RF used for the stack compared to the stacked RF for four selected stations, which locate at different geological units (See Fig. 1 for the detail locations), in Fig.3. They are XE_ES20 in the Lhasa terrane, X4_D14 in the Qiangtang terrane, X4_D10 in the Songpan-Ganzi terrane and X4_H16 in the Tarim Basin, respectively.

2.3 Joint Inversion The combination of RF and SWD intuitively provides tighter constraints on shear-wave velocity structure than inverting each data set independently (e.g., Bodin et al., 2012; Julià et al., 2000; Li et al., 2016). We used the iterative damped least squares method (Julià et al., 2000) to invert the RF and SWD for the crustal and uppermost mantle shear-wave velocity structure beneath each station. We briefly describe the method here. Accounting for the differences in physical units and the number of data points of the two datasets, the joint prediction error is defined as (Julià et al. 2000): S

1    Nr

 Ori  Pri    ri i 1  Nr

    Ns

 Osj  Psj    sj  j 1  Ns

 

(1)

where

is the influence factor, which controls the relative contributions of RF and SWD. It ranges

from 0 to 1, and a smaller value of

corresponds to larger contribution of RF. Nr and Ns are the

length of RF and SWD, respectively. Ori and Pri are the observed and predicted RF at time ti, and the standard deviation is denoted by σri. Osj and Psj are the observed and predicted phase velocities at the jth point, respectively, and the standard deviation is denoted by σrj. We used the programs joint96, surf96 and rftn96 to perform the joint inversion, compute predicted SWD and compute predicted RF, respectively (Herrmann and Ammon, 2004). In our joint inversion, the initial model beneath each station was parameterized as a stack of horizontal layers with layer thicknesses 2 km in the top 50 km, 5 km at depths from 50 km to 100 km, and 10 km below a depth of 100 km. We set the shear-wave velocity in the top 150 km to a constant 4.48 km/s to avoid the assumption of a prefixed Moho interface (e.g., Bao et al., 2015b; Chen et al., 2014; Li et al., 2016). The deeper part of the initial model was taken from the AK135 model (Kennett et al., 1995). The shear-wave velocity was linked to the P-wave velocity using a Vp/Vs ratio of 1.75 and was scaled to obtain density by using the empirical relationship of (Brocher, 2005). The maximum iteration step for each inversion was set to be 30 to make sure the optimal model can be obtained. We determined the Moho depth by searching for the local maximum of the vertical velocity gradient around the velocity range between the average lower crustal and uppermost mantle shear-wave velocities of ~3.8 km/s and ~4.3 km/s, respectively (e.g., Bao et al., 2015b; Guo et al., 2018; Li et al., 2016). In Fig. 4, we show the joint inversion results, as well as the data fitting, for the four selected stations in section 2.2. As shown in Fig. 4, the Moho discontinuities beneath the four stations are clearly identified. Meanwhile, there are clearly low velocity zones shown in the crust.

The predictions fit the observed RFs and SWDs well. We used a bootstrap method (Efron and Tibshirani, 1991) to evaluate the error of our velocity model, like Bao et al. (2015b). Suppose that one station has N RFs for stacking (after selecting based on the procedure described in section 2.2), we randomly selected 1.5*N RFs from the data to produce a stacked RF. We repeated this resampling process 500 times to generate 500 stacked RFs, and then each stacked RF was used individually for joint inversion to obtain 500 velocity models. The inversion results for four selected stations are shown in Fig. S1. It was found that the general features of the resulting velocity structure remained stable and the uncertainty of our inverted velocity models is ~0.05-0.15 km/s depending on depth. The tests demonstrated the stability of the inversion method we used.

3.Result We constructed a 3-D shear-wave velocity model of the crust and uppermost mantle beneath the central and northeastern TP by interpolating the 1-D shear-wave velocity profiles at each station. The horizontal shear-wave velocity slices at 5, 15, 30, 60, 90, and 120 km depths are depicted in Fig. 5. The results clearly show the heterogeneity of the crustal and uppermost mantle shear-wave velocities. As shown in Fig. 5a, at 5 km depth, the Qaidam basin, the Qiangtang terrane, and the northern Lhasa terrane are delineated with low velocities, while other regions show relatively high velocities. The shear-wave velocities of Qiangtang terrane are lower than those of other parts of TP at this depth, which is due to that the Qiangtang terrane is a region of TP with a relatively large number of irresolvable sedimentary basins that tend to be narrow in the north-south direction and elongated

east-west (Yang et al., 2012). At 15 km depth (Fig. 5b), most part of the study region shows relatively low velocities, except the northwestern Qaidam basin and western Qilian orogen. Since there are few stations in the western Qilian orogen (Fig. 1), we did not pay much attention to the structure in this region. At 30 km depth (Fig. 5c), the Qiangtang terrane, the Songpan-Ganzi terrane and the eastern Qinling orogen, and the eastern Qilian orogen are delineated with low velocities. For the Lhasa terrane, most areas show low velocity anomalies, with pervasive high velocity zones in-between. We also observed a LVZ in the eastern Qilian orogen. The Qaidam basin is featured with high velocities. At 60 km depth (Fig. 5d), a notable observation is the significant velocity contrast besides the Kunlun Fault (KLF). The south region to the KLF is featured with widely distributed low velocities, while the north region to the KLF shows high velocities. The velocity contrast is due to the difference in the crustal thickness besides the KLF. The crustal thickness of the south region to the KLF is thicker (> 60 km) than that of the north region to the KLF (< 60 km; see velocity profiles in Fig. 6). The velocities of the south region to the KLF are still crustal velocities; while velocities in surrounding areas already reflect the seismic velocities of the uppermost mantle. At 90 km depth (Fig. 5e), the Qiangtang terrane, the Songpan-Ganzi terrane and the Qinling orogen show low velocities. The Lhasa terrane and Qaidam basin are characterized as relatively high velocities. At 120 km depth (Fig. 5f), the Qinling orogen, the eastern Qilian orogen and the central Songpan-Ganzi terrane show low velocities. The other regions are featured with relatively high velocities with pervasive low velocity zones distribute around and in-between. We show the shear-wave velocity structures along several vertical profiles in Figs. 6 to 8, which can highlight some important observations from our preferred model. The locations of these profiles

are indicated in Fig. 5f. A number of common features are observed, such as the extremely low velocity anomalies in the uppermost crust beneath our study region, especially under the Qaidam basin, and the wide-spread mid-crustal LVZs. As presented in these profiles, LVZs are observed throughout the middle crust of the central TP, although the minimum velocities vary slightly, suggesting the difference in the shear strength of crustal materials. The LVZs are mainly presented at ~15-40 km depths. As shown in Fig. 6, the mid-crustal LVZs are bounded to the south to the KLF and do not appear beneath the Qaidam basin. A notable observation in Fig. 6 is the identification of a continuous low velocity anomaly in the uppermost mantle beneath the eastern KLF, the Songpan-Ganzi terrane and part of the Qiangtang terrane. The low velocity anomaly in the uppermost mantle was also identified by previous a surface wave tomographic study (Zheng et al., 2016). The eastern extent of the mid-crustal LVZs is relatively more complex than their northern extent, as shown in Fig. 7 and Fig. 8. Within our study region, the low velocity zones are confined to the western of Lhasa terrane, while extending along the entire Qiangtang and Songpan-Ganzi terranes. We also observed a small-scale low velocity zone beneath the Qinling orogen around the latitude 35ºN, indicating a certain connection to the low velocity zone in the Songpan-Ganzi terrane at a similar depth range (Fig. 8; profile LL1).

4. Discussion The 3-D shear-wave velocity model discussed in the previous section exhibits prominent low shear wave velocities in the middle crust beneath most part of our study from depths of ~15 km to 40 km, especially beneath the central TP. As shown by the south-north trending vertical profiles in Fig.

6, the northern extension of the LVZs clearly terminate in the front of the KLF, suggesting the KLF to be an important northern boundary of TP. We found no mid-crustal LVZs beneath the Qaidam basin. Instead, a strong crust with relatively high velocities is observed beneath this basin, which may block the northward penetration of the LVZs in TP. In the west-east trending profile LL1 (Fig. 8), we observed a small-scale LVZ beneath the Qinling orogen, showing a certain connection at the similar depth range with the LVZs in the Songpan-Ganzi terrane. This LVZ may be the northward extrusion of the mid-crustal LVZs beneath Songpan-Ganzi terrane across the KLF. A similar finding was also revealed by a previous surface wave tomographic study based on ambient noise (Jiang et al., 2014). The presence of the low velocity zone indicates a weaker crust in the Qinling orogen, which is due in part to the facts that its chemical composition is felsic, as constrained by refraction and reflection explorations (Liu et al., 2006; Pan and Niu, 2011). The felsic rocks are usually weaker than intermediate/mafic rocks under the same temperature and pressure conditions (Jiang et al., 2014). Various previous studies indicate that the middle crust of TP is warm, probably ductile, and that it may contain partial melt (e.g., Bai et al., 2010; Jiang et al., 2014; Le Pape et al., 2012; Xie et al., 2013; Yang et al., 2012; Yao et al., 2008). Basing on the assumption of the pervasive partial melt exits in the middle crust across TP, Hacker et al. (2014) calculated wave speeds within the Tibetan crust based on mineral physical properties and petrological constraints, and then compared these calculations to the crustal shear-wave and radial anisotropy models from surface wave tomography (Xie et al., 2013; Yang et al., 2012). In Hacker et al. (2014), they took the crustal shear-wave speeds and radial anisotropies as observed data, and then set several models to calculate the theoretical wave speeds. They found that the model that contains both horizontally oriented mica and partial melt in

crustal section shows better performance in fitting the observed data. The percent of partial melt is ~2%. They concluded that the middle crust beneath TP is either melting or has melted. The point used in their study for the calculation is located at (93°E, 34°N), which was treated as a typical point for the Tibetan Plateau. This point (green star in Fig. 1) is in the Qiangtang terrane and within our study region. Our model shows that the velocities of the LVZs in the central Tibetan are as low or less than ~3.3 km/s (Figs. 6 to 8), which are consistent with the previous surface wave tomographic studies (e.g., Jiang et al., 2014; Xie et al., 2013; Yang et al., 2012) and are comparable to the velocities at the typical point in Hacker et al. (2014). Assuming a similar thermal status and rock compositions in the central and northeastern TP, the mid-crustal low velocities in our preferred model may suggest that the middle crust in our study region should be actively melting or has undergone melting in the past. The percent of partial melt is inferred to ~2%, as was suggested by Hacker et al. (2014). A recent study (Wang et al., 2012) found the tourmaline-bearing mica and biotite rhyolites in the Bukadaban-Malanshan area, near the northern margin of TP. The authors suggested that the Bukadaban-Malanshan rhyolites were generated by dehydration melting of metasedimentary rocks under certain conditions (Wang et al., 2012). Their data not only confirm the occurrence of a partially molten zone in the mid-to-lower crust beneath northern TP, but also constrain the crustal melting to have existed from middle Miocene to Quaternary times. This is the direct petrologic evidence of the existence of partial melt in the Tibetan middle crust, as well as other evidences from magnetotelluric studies (e.g., Bai et al., 2010; Le Pape et al., 2012) and seismic tomography studies (e.g., Jiang et al., 2014; Yang et al., 2012; Zheng et al., 2016). The above discussions show that our model strongly supports the existence of partial melt in the

middle crust beneath the central and northeastern TP, implying the possible contribution of the crustal channel flow model to the deformation of TP. Additionally, as shown in Fig. 6, we also identified a continuous low velocity anomaly in the uppermost mantle beneath the eastern KLF, the Songpan-Ganzi terrane and part of the Qiangtang terrane, which was also clearly identified by Zheng et al. (2016) using ambient noise tomography. The low velocity anomaly is consistent with the detachment of the lithosphere root that has completely or partially removed the mantle lithosphere, suggesting an induced local mantle upwelling after the detachment of the lithosphere (Mo et al., 2005; Zheng et al., 2016).

5. Conclusion Based on the joint analysis of RF and SWD, we derived a 3-D lithospheric shear wave velocity model in the central and northeastern TP. Our new model provides new insights into the crust and uppermost mantle structure, especially the distribution of possible partial melts under the central and northeastern TP. The present model indicates wide-spread mid-crustal LVZs beneath this region, which are mainly confined to the south region to the KLF and at the depths ranging from ~15 km to 40 km. The shear velocities in the low velocity zones are less than ~3.3 km/s. By referring to the conclusions of Hacker et al. (2014), the mid-crustal LVZs in our preferred model may suggest that the crust beneath the central and northeastern TP should be actively melting or has undergone melting in the past and further imply the possible contribution of the crustal channel flow model to the deformation of TP.

Acknowledgement The seismic waveform data used in this study were downloaded from IRIS DMC (X4, doi.org/10.7914/SN/X4_2007; XE, doi.org/10.7914/SN/XE_2003; XD doi.org/10.7914/SN/XD_2002; ZV, doi.org/10.7914/SN/ZV_2008; CB, doi.org/10.7914/SN/CB). We greatly thank Professor Vernon Cormier (the editor) and two anonymous reviewers for their insightful comments and suggestions that helped to improve this manuscript. We thank Herrmann R.B. for making the CPS (Computer Programs in Seismology) package (Herrmann and Ammon, 2004) freely available. We also thank Dr. X.W. Bao for providing us the surface wave dispersion data. We thank the Wiley’s language service for language editing. This study is jointly supported by the National Natural Science Foundation of China (No. 41704046; 41801389), the China Postdoctoral Science Foundation (No. 2016M602343). The calculation in this paper have been done on the superconducting system in the Supercomputing Center of Wuhan University.

References Avouac, J.-P., Tapponnier, P., 1993. Kinematic Model of the Active Deformation in Central Asia. Geophys. Res. Lett. 20, 895-898. Bai, D., Unsworth, M.J., Meju, M.A., Ma, X., Teng, J., Kong, X., Sun, Y., Sun, J., Wang, L., Jiang, C., Zhao, C., Xiao, P., Liu, M., 2010. Crustal deformation of the eastern Tibetan plateau revealed by magnetotelluric imaging. Nat. Geosci. 3, 358-362. Bao, X., Song, X., Li, J., 2015a. High-resolution lithospheric structure beneath Mainland China from ambient noise and earthquake surface-wave tomography. Earth Planet. Sci. Lett. 417, 132-141.

Bao, X., Sun, X., Xu, M., Eaton, D.W., Song, X., Wang, L., Ding, Z., Mi, N., Li, H., Yu, D., Huang, Z., Wang, P., 2015b. Two crustal low-velocity channels beneath SE Tibet revealed by joint inversion of Rayleigh wave dispersion and receiver functions. Earth Planet. Sci. Lett. 415, 16-24. Bilham, R., Larson, K., Freymueller, J., Jouanne, F., Le Fort, P., Leturmy, P., Mugnier, J.L., Gamond, J.F., Glot, J.P., Martinod, J., Chaudury, N.L., Chitrakar, G.R., Gautam, U.P., Koirala, B.P., Pandey, M.R., Ranabhat, R., Sapkota, S.N., Shrestha, P.L., Thakuri, M.C., Timilsina, U.R., Tiwari, D.R., Vidal, G., Vigny, C., Galy, A., De Voogd, B., 1997. GPS measurements of present-day convergence across the Nepal Himalaya. Nature 386, 61-64. Bodin, T., Sambridge, M., Tkalćić, H., Arroucau, P., Gallagher, K., Rawlinson, N., 2012. Transdimensional inversion of receiver functions and surface wave dispersion. J. Geophys. Res. Solid Earth 117, B02301. Brocher, T.M., 2005. Empirical relations between elastic wavespeeds and density in the Earth’s crust. Bull. Seismol. Soc. Am. 95, 2081-2092. Ceylan, S., Ni, J., Chen, J.Y., Zhang, Q., Tilmann, F., Sandvol, E., 2012. Fragmented Indian plate and vertically coherent deformation beneath eastern Tibet. J. Geophys. Res. B Solid Earth 117, 1-11. Chen, H. peng, Zhu, L. bao, Wang, Q. dong, Zhang, P., Yang, Y. hang, 2014. S-wave velocity structure of the North China from inversion of Rayleigh wave phase velocity. J. Asian Earth Sci. 88, 178-191. Chen, H., Zhu, L., Su, Y., 2016. Low velocity crustal flow and crust–mantle coupling mechanism in Yunnan, SE Tibet, revealed by 3D S-wave velocity and azimuthal anisotropy. Tectonophysics 685, 8-20.

Clark, M.K., Royden, L.H., 2000. Topographic ooze: Building the eastern margin of Tibet by lower crustal flow. Geology 28, 703-706. Efron, B., Tibshirani, R., 1991. Statistical data analysis in the computer age. Science 253, 390-395. England, P., Houseman, G., 1986. Finite Strain Calculations of Continental Deformation 2. Comparison With the India_Asia Collision Zone. J. Geophys. Res. 91, 3664-3676. Fontaine, F.R., Barruol, G., Tkalčić, H., Wölbern, I., Rümpker, G., Bodin, T., Haugmard, M., 2015. Crustal and uppermost mantle structure variation beneath La Réunion hotspot track. Geophys. J. Int. 203, 107-126. Fontaine, F.R., Tkalčić, H., Kennett, B.L.N., 2013. Imaging crustal structure variation across southeastern Australia. Tectonophysics 582, 112-125. Gan, W., Zhang, P., Shen, Z.K., Niu, Z., Wang, M., Wan, Y., Zhou, D., Cheng, J., 2007. Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements. J. Geophys. Res. Solid Earth 112, 1-14. Guo, Z., Gao, X., Li, T., Wang, W., 2018. Crustal and uppermost mantle structures of the South China from joint analysis of receiver functions and Rayleigh wave dispersions. Phys. Earth Planet. Inter. 278, 16-25. Hacker, B.R., Ritzwoller, M.H., Xie, J., 2014. Partially melted, mica-bearing crust in Central Tibet. Tectonics 33, 1408-1424. Herrmann, R. B., and C. J. Ammon (2004), Surface waves, receiver functions and crustal structure, in Computer Programs in Seismology, Version 3.30, Saint Louis University. [Available at http://www.eas.slu.edu/People/RBHerrmann/CPS330.html.]

Hung, S.H., Chen, W.P., Chiao, L.Y., Tseng, T.L., 2010. First multi-scale, finite-frequency tomography illuminates 3-D anatomy of the Tibetan Plateau. Geophys. Res. Lett. 37, 1-5. Jiang, C., Yang, Y., Rawlinson, N., Griffin, W.L., 2016. Crustal structure of the Newer Volcanics Province, SE Australia, from ambient noise tomography. Tectonophysics 683, 382-392. Jiang, C., Yang, Y., Zheng, Y., 2014. Penetration of mid-crustal low velocity zone across the Kunlun Fault in the NE Tibetan Plateau revealed by ambient noise tomography. Earth Planet. Sci. Lett. 406, 81-92. Julià, J., Ammon, C.J., Herrmann, R.B., Correig, A.M., 2000. Joint inversion of receiver function and surface wave dispersion observations. Geophys. J. Int. 143, 99-112. Kennett, B.L.N., Engdahl, E.R., Buland, R., 1995. Constraints on seismic velocities in the Earth from travel times. Geophys. J. Int. 122, 108-124. Kind, R., Sobolev, S. V, Yuan, X., 1999. Seismic Evidence for a Detached Indian Lithospheric Mantle Beneath Tibet . Seismic Evidence for a Detached Indian Lithospheric Mantle Beneath Tibet. Science (80-. ). 283, 1306-1310. Lamarque, G., Barruol, G., Fontaine, F.R., Bascou, J., Menot, R.-P., 2015. Crustal and mantle structure beneath the Terre Adelie Craton, East Antarctica: insights from receiver function and seismic anisotropy measurements. Geophys. J. Int. 200, 807–821. Le Pape, F., Jones, A.G., Vozar, J., Wenbo, W., 2012. Penetration of crustal melt beyond the Kunlun Fault into northern Tibet. Nat. Geosci. 5, 330-335. Li, C., van der Hilst, R.D., Meltzer, A.S., Engdahl, E.R., 2008. Subduction of the Indian lithosphere beneath the Tibetan Plateau and Burma. Earth Planet. Sci. Lett. 274, 157-168.

Li, L., Li, A., Shen, Y., Sandvol, E.A., Shi, D., Li, H., Li, X., 2013. Shear wave structure in the northeastern Tibetan Plateau from Rayleigh wave tomography. J. Geophys. Res. Solid Earth 118, 4170-4183, doi:10.1002/jgrb.50292. Li, M., Zhang, S., Bodin, T., Lin, X., Wu, T., 2018a. Transdimensional inversion of scattered body waves for 1D S-wave velocity structure - Application to the Tengchong volcanic area, Southwestern China. J. Asian Earth Sci. 159, 60-68. Li, M., Zhang, S., Wang, F., Wu, T., Qin, W., 2016. Crustal and upper-mantle structure of the southeastern Tibetan Plateau from joint analysis of surface wave dispersion and receiver functions. J. Asian Earth Sci. 117, 52-63. Li, M., Zhang, S., Wu, T., Hua, Y., Zhang, B., 2018b. Fine crustal and uppermost mantle S-wave velocity structure beneath the Tengchong volcanic area inferred from receiver function and surface-wave dispersion : constraints on magma chamber distribution. Bull. Volcanol. 80. Li, Y., Pan, J., Wu, Q., Ding, Z., 2017. Lithospheric structure beneath the northeastern Tibetan Plateau and the western Sino-Korea Craton revealed by Rayleigh wave tomography. Geophys. J. Int. 210, 570-584. Liang, X., Sandvol, E., Chen, Y.J., Hearn, T., Ni, J., Klemperer, S., Shen, Y., Tilmann, F., 2012. A complex Tibetan upper mantle: A fragmented Indian slab and no south-verging subduction of Eurasian lithosphere. Earth Planet. Sci. Lett. 333–334, 101-111. Ligorría, J.P., Ammon, C.J., 1999. Iterative deconvolution and receiver function estimation. Bull. Seismol. Soc. Am. 89, 1395-1400. Liu, M., Mooney, W.D., Li, S., Okaya, N., Detweiler, S., 2006. Crustal structure of the northeastern

margin of the Tibetan plateau from the Songpan-Ganzi terrane to the Ordos basin. Tectonophysics 420, 253-266. Liu, Q.Y., Van Der Hilst, R.D., Li, Y., Yao, H.J., Chen, J.H., Guo, B., Qi, S.H., Wang, J., Huang, H., Li, S.C., 2014. Eastward expansion of the Tibetan Plateau by crustal flow and strain partitioning across faults. Nat. Geosci. 7, 361-365. Mo, X., Dong, G., Zhao, Z., Zhou, S., Wang, L., Qiu, R., Zhang, F., 2005. Spatial and temporal distribution and characteristics of Granitoids in the Gangdese, Tibet and implication for crustal growth and evolution. Geol. J. China Univ., 11(3), 281–290. Murodov, D., Zhao, J., Xu, Q., Liu, H., Pei, S., 2018. Complex N–S variations in Moho depth and Vp/Vs ratio beneath the western Tibetan Plateau as revealed by receiver function analysis. Geophys. J. Int. 214, 895-906. Nunn, C., Roecker, S.W., Priestley, K.F., Liang, X., Gilligan, A., 2014. Joint inversion of surface waves and teleseismic body waves across the Tibetan collision zone: The fate of subducted Indian lithosphere. Geophys. J. Int. 198, 1526-1542. Owens, T.J., Zandt, G., 1997. Implications of crustal property variations for models of Tibetan plateau evolution. Nature. 387, 37-43. https://doi.org/10.1038/387037a0 Pan, S., Niu, F., 2011. Large contrasts in crustal structure and composition between the Ordos plateau and the NE Tibetan plateau from receiver function analysis. Earth Planet. Sci. Lett. 303, 291-298. Pan, Y., Shen, W. Bin, 2017. Contemporary crustal movement of southeastern Tibet: Constraints from dense GPS measurements. Sci. Rep. 7, 1-7. Pan, Y., Shen, W.-B., Shum, C.K., Chen, R., 2018. Spatially varying surface seasonal oscillations and

3-D crustal deformation of the Tibetan Plateau derived from GPS and GRACE data. Earth Planet. Sci. Lett. 502, 12-22. Pandey, S., Yuan, X., Debayle, E., Tilmann, F., Priestley, K., Li, X., 2015. Depth-variant azimuthal anisotropy in Tibet revealed by surface wave tomography. Geophys. Res. Lett. 42, 4326-4334. Priestley, K., Debayle, E., McKenzie, D., Pilidou, S., 2006. Upper mantle structure of eastern Asia from multimode surface waveform tomography. J. Geophys. Res. Solid Earth 111, 1-20. Rowley, D.B., 1996. Age of initiation of collision between India and Asia: A review of stratigraphic data. Earth Planet. Sci. Lett. 145, 1-13. Royden, L.H., 1997. Surface Deformation and Lower Crustal Flow in Eastern Tibet. Science (80-. ). 276, 788-790. Royden, L.H., Burchfiel, B.C., Hilst, R.D. Van Der, 2008. The Geological Evolution of the Tibetan Plateau. Science (80-. ). 321, 1054-1058. Sandvol, E., Chen, J., Ni, J., Zhou, S., Ma, Y., Zhang, X., Yue, H., Ceylan, S., Bao, X.C., Hearn, T., 2008. Lithospheric Seismic Velocity Structure of the Northern Tibetan Plateau: The ASCENT Seismic Experiment. EOS Trans AGU 89, Fall Meeting Supplement (abstract S14A-01). Sol, S., Meltzer, A., Burgmann, R., Hilst, R.D, v.d., King, R., Chen, Z., Koons, P.O., Lev, E., Liu, Y.P., Zeitler, P.K., Zhang, X., Zhang, J., Zurek, B., 2007. Geodynamics of the southern Tibetan Plateau from seismic anisotropy and geodesy. Geology 35, 563–566. Sosa, A., Thompson, L., Velasco, A.A., Romero, R., Herrmann, R.B., 2014. 3-D structure of the Rio Grande Rift from 1-D constrained joint inversion of receiver functions and surface wave dispersion. Earth Planet. Sci. Lett. 402, 127-137.

Tapponnier, P., 2001. Oblique Stepwise Rise and Growth of the Tibet Plateau. Science (80-. ). 294, 1671-1677. Tilmann, F., Ni, J., Team, I.I.S., 2003. Seismic imaging of the Indian continental lithosphere. Science (80-. ). 300, 1424-1427. Tkalčić, H., Chen, Y., Liu, R., Zhibin, H., Sun, L., Chan, W., 2011. Multistep modelling of teleseismic receiver functions combined with constraints from seismic tomography: Crustal structure beneath southeast China. Geophys. J. Int. 187, 303-326. Tkalčić, H., Pasyanos, M.E., Rodgers, A.J., Gök, R., Walter, W.R. and Al-Amri, A., 2006. A multistep approach for joint modeling of surface wave dispersion and teleseismic receiver functions: Implications for lithospheric structure of the Arabian Peninsula. J. Geophys. Res. Solid Earth 111, 1-25, doi: 10.1029/2005JB004130. Wang, Q., Chung, S.L., Li, X.H., Wyman, D., Li, Z.X., Sun, W.D., Qiu, H.N., Liu, Y.S., Zhu, Y.T., 2012. Crustal melting and flow beneath northern Tibet: Evidence from mid-miocene to quaternary strongly peraluminous rhyolites in the Southern Kunlun Range. J. Petrol. 53, 2523-2566. Wei, X., Jiang, M., Liang, X., Chen, L., Ai, Y., 2017. Limited southward underthrusting of the Asian lithosphere and material extrusion beneath the northeastern margin of Tibet, inferred from teleseismic Rayleigh wave tomography. J. Geophys. Res. Solid Earth 122, 7172-7189. Wu, T., Zhang, S., Li, M., Qin, W., Zhang, C., 2016. Two crustal flowing channels and volcanic magma migration underneath the SE margin of the Tibetan Plateau as revealed by surface wave tomography. J. Asian Earth Sci. 132, 25-39.

Xie, J., Ritzwoller, M.H., Shen, W., Yang, Y., Zheng, Y., Zhou, L., 2013. Crustal radial anisotropy across eastern Tibet and the Western Yangtze Craton. J. Geophys. Res. Solid Earth 118, 4226-4252. Xu, C., Luo, Z., Sun, R., Zhou, H., Wu, Y., 2018. Multilayer densities using a wavelet-based gravity method and their tectonic implications beneath the Tibetan Plateau. Geophys. J. Int. 213, 2085-2095. Xu, L., Rondenay, S., van der Hilst, R.D., 2007. Structure of the crust beneath the southeastern Tibetan Plateau from teleseismic receiver functions. Phys. Earth Planet. Inter. 165, 176-193. Yang, Y., Liu, M., 2013. The indo-asian continental collision: A 3-D viscous model. Tectonophysics 606, 198-211. Yang, Y., Ritzwoller, M.H., Zheng, Y., Shen, W., Levshin, A.L., Xie, Z., 2012. A synoptic view of the distribution and connectivity of the mid-crustal low velocity zone beneath Tibet. J. Geophys. Res. Solid Earth 117, 1-20. Yao, H., Beghein, C., Van Der Hilst, R.D., 2008. Surface wave array tomography in SE Tibet from ambient seismic noise and two-station analysis - II. Crustal and upper-mantle structure. Geophys. J. Int. 173, 205-219. Yin, A., Harrison, M., 2000. Geological evolution of the Himalayan- Tibetan orogen. Annu. Rev. Earth planet. Sci., 28, 211-280. Zhang, Q., Sandvol, E., Ni, J., Yang, Y., Chen, Y.J., 2011. Rayleigh wave tomography of the northeastern margin of the Tibetan Plateau. Earth Planet. Sci. Lett. 304, 103-112. Zhao, W., Brown, L., Wu, Z., Klemperer, S.L., Shi, D., Mechie, J., Su, H., Tilmann, F., Karplus, M.S.,

Makovsky, Y., 2008. Seismology across the northeastern edge of the Tibetan Plateau. EOS Trans. AGU 89, 487. Zheng, D., Li, H., Shen, Y., Tan, J., Ouyang, L., Li, X., 2016. Crustal and upper mantle structure beneath the northeastern Tibetan Plateau from joint analysis of receiver functions and Rayleigh wave dispersions. Geophys. J. Int. 204, 583-590. Zheng X.F., Yao, Z.X., Liang, J.H., Zheng, J. The role played and opportunities provided by IGP DMC of China National Seismic Network in Wenchuan earthquake disaster relief and researches. Bull. Seismol. Soc. Am., 2010, 100(5B): 2866-2872. Zhu, L., Owens, T., Randall, G.E., 1995. Lateral Variation in Crustal Structure of the Northern Tibetan Plateau Inferred from Teleseismic Receiver Functions. Bull. Seismol. 85, 1531-1540.

Table Captions Table 1 Broadband seismic networks used in this study

Figure Captions Fig. 1. Topographic map of the study area showing the main tectonic units in the central and northeastern Tibetan Plateau, and the distributions of seismic stations. Symbols for each seismic network are defined in Table 1. The gray dashed lines represent faults and the boundaries of geological units, respectively. The green star represents the location of the selected point within Qiangtang terrane in Hacker et al. (2014). The abbreviations are: Indus-Yarlung suture (IYS), Bangong-Nujiang suture (BNS), Jinsha suture (JS), Kunlun Fault (KLF), Haiyuan Fault (HF), Altyn Tagh Fault (ATF).

Fig. 2. Epicenter locations of the earthquakes after the preliminary quality control of RFs (blue circles) with moment magnitudes between 5.5 and 7.5 and epicentral distances within 30° and 90° recorded. The red triangles denote the locations of the seismic stations. See the text for the explanation for the preliminary quality control of RFs.

Fig. 3. The individual RF (thin gray curves) used for the stack compared to the stacked RF (thick black curves) at four selected stations. The short dashed red lines mark the arrival times of the Pms phase. (a) XE_ES20 in the Lhasa terrane; (b) X4_D14 in the Qiangtang terrane; (c) X4_D10 in the Songpan-Ganzi terrane; (d) X4_H16 in the Tarim Basin.

Fig. 4. The joint inversion results for four selected stations. The panels for each station are the preferred

1-D shear-wave velocity model, the observed (blue curve) and predicted (red curve) RFs, and the observed (red circle) and predicted (black cure) phase dispersion curves, respectively. The short dashed red lines mark the arrival times of the Pms phase. (a) XE_ES20 in the Lhasa terrane; (b) X4_D14 in the Qiangtang terrane; (c) X4_D10 in the Songpan-Ganzi terrane; (d) X4_H16 in the Tarim Basin.

Fig. 5. Preferred shear-wave velocity model at 5, 15, 30, 60, 90 and 120 km depths. The abbreviations in a) are the same as that in Fig. 1. The blue dashed lines in f) represent the locations of the vertical velocity profiles in Figs. 6-8.

Fig. 6. Six south-north trending vertical profiles from across our study region (along the longitudes from 92°Eto 97°E). The locations are denoted in Fig. 5f. The thick white curves represent the estimated Moho interfaces. Topography is plotted above each profile. The seismic stations on each topographic profile are marked as black inverted triangles.

Fig. 7. Four west-east trending vertical profiles from across our study region (along the latitudes from 30°Nto 33°N). The locations are denoted in Fig. 5f. The thick white curves represent the estimated Moho interfaces. Topography is plotted above each profile. The seismic stations on each topographic profile are marked as black inverted triangles.

Fig. 8. Three west-east trending vertical profiles from across our study region (along the latitudes

from 34°Nto 36°N). The locations are denoted in Fig. 5f. The thick white curves represent the estimated Moho interfaces. Topography is plotted above each profile. The seismic stations on each topographic profile are marked as black inverted triangles.

Figures

Fig. 1.

Fig. 2.

Fig. 3.

Fig. 4.

Fig. 5.

Fig. 6.

Fig. 7.

Fig. 8.

Table Table 1 Broadband seismic networks used in this study Network

Network Code

Symbol in Fig. 1

Namche Barwa

XE

red triangle

INDEPTH IV/ASCENT

X4

blue triangle

Broadband Seismic Experiment in Various Regions in China (Central China)

XD

red square

INDEPTH IV/UK

XO

yellow square

Northeast Tibet Plateau Seismic Experiment (NETS)

ZV

blue square

China National Seismic Network

CB

red star

Highlights: 

Lithospheric model in the central and northeastern Tibetan is proposed by joint inversion



Widely spread mid-crustal partial melts are observed in the central and northeastern Tibetan



Crustal channel flow model may contribute largely to the Tibetan deformation