Lithospheric structures of and tectonic implications for the central–east Tibetan plateau inferred from joint tomography of receiver functions and surface waves

SUMMARY

1 INTRODUCTION AND GEOLOGY

The Qinghai–Tibetan plateau (or Tibetan plateau) is the largest and highest plateau on Earth, with a surface elevation primarily >3000 m above sea level, except for the Qaidam Basin (Fig. 1). The plateau formed via the successive amalgamation of continental terranes (Fig. 1) during the Palaeozoic and Mesozoic, and has undergone extensive shortening, crustal thickening, and uplift since the Mesozoic (e.g. Dewey & Bird 1970; Molnar & Tapponnier 1975; Dewey et al. 1988), following collision with the Indian Plate to the south. Within the plateau, a series of E–W striking tectonic sutures/faults, such as the Main Frontal Thrust, the Yarlung–Zangbo suture (YZS), the Bangong–Nujiang suture (BNS), the Jinsha River suture (JRS) and the Kunlun Fault (KF), separate the plateau into several tectonic terranes, including the Himalaya, Lhasa, Qiangtang and Songpan–Ganzi terranes and the Kunlun–Qilian Orogen from south to north. Outside the plateau, several stable tectonic blocks surround it, including the Indian Plate in the south, the South China Craton (SCC) in the east, the North China Craton (NCC) in the northeast and the Tarim Craton (or Tarim basin) in the north. Between the plateau and the surrounding stable blocks are several large boundary faults or mountain ranges, including the Main Frontal Thrust (MFT) or Himalayan Mountains in the south, the Hengduan Mountains in the southeast, the Longmen Shan or Longmen Shan Fault (LMSF) in the east, the Qilian Shan or Qilian Shan Frontal Thrust (QSFT) in the northeast and the Altyn mountains or Altyn Tagh Fault (ATF) in the north. The ongoing Indian–Eurasian collision has led to the formation or deformation of these sutures/faults within the Tibetan plateau or the boundary faults and mountain ranges around the plateau (Yin & Harrison 2000). Imaging 3-D lithospheric structures around the central–east Tibetan plateau will hopefully promote our understanding of regional tectonic evolution (e.g. Guillot & Replumaz 2013; Wang 2013).

Figure 1.

Topography and tectonic setting of the study region. The surface topography is from the ETOPO2 data set (National Geophysical Data Center 2006). Plate boundaries (indented lines) are from www-udc.ig.utexas.edu. Boundaries between major crustal blocks (thick grey dashed lines or thick black lines), and important regional faults (thin black lines) are simplified from previous regional studies (Ren et al. 1999; Ren & Liu 2002; Sone & Metcalfe 2008; Pan et al. 2012). Thick white arrow denotes the current movement direction of the Indian Plate at 90°E, 26°N relative to Eurasia (Kreemer et al. 2014). Inset map shows the study region in Asia. Abbreviations are as follows: AS = Ailaoshan Suture, ATF = Altyn Tagh Fault, BM = Burma Microplate, BNS = Bangong–Nujiang suture, CAOB = Central Asian Orogenic Belt, IT = Indochina Terrane, JRS = Jinsha River Suture, KF = Kunlun Fault, LMSF = Longmenshan Fault, LS = Lanchangjiang Suture, MFT = Main Frontal Thrust, NCC = North China Craton, NQT = North Qaidam Thrust, NT = Naga Thrust, QSFT = Qilian Shan Frontal Thrust, RRF = Red River Fault, SCC = South China Craton, SF = Sagaing Fault, SGT = Songpan–Ganzi Terrane, YZS = Yarlung Zangbo Suture.

Open in new tabDownload slide

Broad-band seismic data are a good choice to detect deep earth structures. In the region of the Tibetan plateau, other than permanent seismic stations primarily maintained by the provincial and national governments of China (Zheng et al. 2010), extensive seismological studies have been undertaken since the 1990s by deploying numerous temporary broad-band seismic arrays (e.g. INDEPTH II, INDEPTH III, INDEPTH IV/ASCENT, HI-CLIMB and HIMNT, e.g. Zhao et al. 1993; Nelson et al. 1996; Nabelek et al. 2005; Zhao et al. 2010; Yue et al. 2012; Feng et al. 2014; Wang et al. 2018b; Karplus et al. 2019). Seismic arrays on the southern edge of the plateau along the Himalayan Orogen have primarily identified the deformed structures of the collisional belt between the Indian subcontinent and the Tibetan plateau (e.g. Zhao et al. 1993; Brown et al. 1996; Nelson et al. 1996; Nabelek et al. 2009). Further seismic array studies have extended to the northern and eastern edges of the Tibetan plateau, providing key details on the deep-seated structures and tectonic evolution of Tibet (e.g. Kind et al. 2002; Mechie et al. 2011; Liang et al. 2012; Mechie et al. 2012; Yue et al. 2012; Agius & Lebedev 2013; Feng et al. 2014; Huang et al. 2014; Feng et al. 2017). However, non-uniform data coverage and the diversity of seismological techniques applied in Tibet have led to uncertainties and discrepancies regarding the tectonic evolution of the region, including: (i) the geometry of the underthrusting Indian lithosphere beneath the Tibetan plateau; (ii) whether or not and if so, how far has the Asian Plate subducted southward/southwestward beneath the Tibetan plateau (Kind et al. 2002; Liang et al. 2012; Feng et al. 2014) and (iii) how the interplay of the Tibetan plateau with its surrounding blocks affects the mode of deformation of the plateau (Burchfiel et al. 1995; Tapponnier et al. 2001; Royden et al. 2008). While most of the portable seismic stations have been deployed in linear arrays to address specific research questions, the amalgamation of all these data sets yields improved data coverage over the central–eastern plateau (Fig. 2). Regional 3-D seismic tomography of a broad swath of seismic observations can better address the above questions than 2-D linear array studies.

Figure 2.

Seismic stations used in this study. Open squares are for Rayleigh wave data, filled triangles are for receiver function data. Magenta circles mark the Ps piercing points of receiver-function observations.

Open in new tabDownload slide

Body wave receiver functions (RF) and surface wave dispersions (SWD) are two primary seismic measurements taking different Earth structure information. RF are sensitive to S-wave velocity (Vs) contrasts between neighbouring layers, and are thus suited to detecting interfaces between layers, such as the Moho (e.g. Yue et al. 2012; Feng et al. 2014), while SWD are sensitive to average Vs down to a certain depth and are thus suited for imaging regional absolute Vs distribution (e.g. Yao et al. 2008; Feng et al. 2011; Wei et al. 2017). RF primarily recover structural information in the immediate vicinity of seismic stations, and they are commonly applied to 1-D profile studies or 2-D migration imaging, whereas SWD are ideally suited to 3-D tomographic studies. The intrinsic trade-off between crustal velocity and crustal thickness always influences the reliability of structures in the upper mantle. Joint inversion of both RF and SWD can hopefully better define both Vs contrasts between interfaces and Vs distributions beneath the region, as has been done early in many previous studies that used joint inversions to derive 1-D Vs models (Özalaybey et al. 1997; Julià et al. 2000; Du et al. 2002; An & Assumpção 2004; Chang et al. 2004; Lawrence & Wiens 2004). On the basis of the method of Feng & An (2010), we design a joint tomographic technique to simultaneously and directly invert for a 3-D Vs model from RF waveform and originally measured SWD.

Here, we apply the joint tomographic technique to a comprehensive seismic data set covering the central–east Tibetan plateau to better map its deep seismic structure, from which more informed inferences regarding the tectonic evolution of the Tibetan plateau and its surrounding regions can be made.

2 DATA AND PROCESSING

Rayleigh-wave dispersions and P-wave receiver functions are the observations in the 3-D inversion. Rayleigh wave dispersions are retrieved from the vertical component of both earthquake surface waves and cross-correlated interstation Green's functions from ambient noise. All the data are recorded from ∼900 permanent and portable seismic stations (squares in Fig. 2). The data of the Chinese permanent stations are provided by the Data Management Centre of the China National Seismic Network (Data et al. 2007; Zheng et al. 2010). Most data are from portable seismic stations of temporary experiments in the region, including the INDEPTH-III (Kind et al. 2002), INDEPTH-IV/ASCENT (Zhao et al. 2011; Yue et al. 2012) and Qilian Shan seismic arrays (Feng et al. 2014), and a seismic array in the western North China Craton (Feng et al. 2017), among others available from the Data Management Center of the Incorporated Research Institutions for Seismology.

For earthquake data, only waveforms for events where the theoretical amplitudes estimated from magnitude and epicentral distance were larger than the general noise level were preliminarily chosen for the analysis, but many waveforms with low signal-to-noise ratio were discarded in the later processes. For noise data, desampled (1 sps) continuous data available for provincial seismic stations and the INDEPTH-IV and Qilian Shan seismic arrays were processed. Rayleigh wave group-velocities were measured by the Multiple Filtering Technique (MFT, Dziewonski et al. 1969) with the aid of the program do_mft from the package of Herrmann (2013). Previously processed surface wave group-velocity measurements in and around central–east Tibet (Feng & An 2010; Feng et al. 2011; Li et al. 2014a) were also included in this study. Figs 3(a) and S1(a) show, respectively, the number and the average epicentral distance of group-velocity measurements at different periods. The total number of measurements reaches ∼40 000 at the 20 s period and gradually decreases for shorter and longer periods (Fig. 3a). Figs 3(b) and S1(b) show, respectively, ray azimuthal coverage and density of group-velocity measurements at the period of 50 s. There is a good lateral azimuthal coverage of surface wave data in eastern Tibet (Fig. 3b). The ray density (Fig. S1b) of the surface wave data at the 50 s period, defined as the number of group-velocity measurements in each 0.5° × 0.5° area, is >200 for most of the study region, and decreases dramatically to ∼10 along the western, southern and northern margins of the region.

Figure 3.

Rayleigh wave group-velocity dispersions used in this study. (a) Number of Rayleigh-wave group-velocity measurements for each period. The terms ‘event regional’ and ‘noise regional’ indicate regional measurements with paths completely within the study region from earthquake data and interstation Green's functions, respectively. The terms ‘event teleseismic’ and ‘noise teleseismic’ indicate measurements with only a portion of their paths travelling through the study region from earthquake data and interstation Green's functions, respectively. The term ‘total’ denotes the summation of the four types of measurements. (b) Path azimuthal coverage in 10° azimuthal bins in 0.5° × 0.5° cells for Rayleigh waves at the period of 50 s. The length of the red line represents the number of paths in the 10° azimuthal interval. The other symbols are the same as in Fig. 1.

Open in new tabDownload slide

Teleseismic events with epicentral distances of 30°–95° and magnitudes greater than 5.5 are considered for RF analysis. After rotating the data from the vertical–north–east (ZNE) to the vertical–radial–transverse (ZRT) coordinate system, an iterative time domain deconvolution (Ligorría & Ammon 1999) of the vertical component from the radial component is used to retrieve receiver functions. The Computer Programs in Seismology package by Herrmann (2013) was utilized to process the RF. A larger width parameter of the Gaussian filter in the deconvolution produces higher-frequency RF, which may not only decrease the amplitude of multiples from the Moho (Julià 2007) but also contain minor converted waves and multiples from minor layers, which then produce localized strongly oscillating Vs anomalies in the inverted model. We thus used a Gaussian-filter width parameter of 1, corresponding to a low-pass filtering at ∼0.5 Hz. Each individual receiver function is visually inspected, and traces with no obvious P-wave to S-wave (Ps) conversion at the Moho or with strong long-period oscillations are omitted. We retrieved 26 859 RF for 670 stations in eastern Tibet (triangles in Fig. 2). Fig. S2 provides examples of the P-wave receiver functions for 4 seismic stations. RF waveforms within a delay time window of –2 to 50 s are considered in the subsequent inversion. Such a time-window length is set to include crustal converted phases (e.g. Ps) and most of their multiples (e.g. PpPs and PsPs/PpSs) even for very thick crust such as in the Tibetan plateau and exclude fake or weak mantle phases that may lead to biased mantle structures.

3 METHODOLOGY

The joint RF–SWD inversion to produce a 3-D seismic velocity model is designed on the basis of independent linearized inversions of the SWD and RF data for a 3-D model. We therefore first introduce the key steps of the inversions of SWD and of RF, and then present our joint tomographic inversion methodology.

3.1 Relationship between surface wave dispersions and 3-D Vs

SWD are primarily sensitive to Vs distributions, making them an ideal candidate to invert for Vs structure. Surface wave tomographic inversions are generally partitioned into two main steps: period-by-period 2-D tomographic inversions for regionalized SWD and cell-by-cell inversions of regionalized SWD for 1-D Vs profiles (e.g. Pasyanos et al. 2001; Villaseñor et al. 2001; Ritzwoller et al. 2002; Huang et al. 2003; Feng et al. 2004). Unlike traditional inversions of surface wave tomography, Feng & An (2010) developed a single-step surface wave tomographic method that combines the traditional two-step formulations into one formula to directly invert for a 3-D Vs model. The method enables simultaneously the use of regional data (with rays completely inside the study region) and teleseismic data (with rays partly inside the study region), thus making it applicable to studies at continental scale (Feng & An 2010; An et al. 2015), regional scale (An et al. 2009; Feng et al. 2011) and local scale (Ramos et al. 2016). Furthermore, any a priori constraints and different kinds of observations can be readily added to the 3-D model inversion system. We thus adopted the methodology of Feng & An (2010) to invert for 3-D Vs structure here. The basic concepts of this approach are summarized below, and a detailed description can be found in Feng & An (2010).

In eq. (3), all the variables are known for a given reference model, with the exception of Δβ_ij. This leads to the direct relationship between Δt and Δβ_ij. Eq. (3) only represents one constraint for one period in one surface wave dispersion curve. By using more dispersion curves over more periods, a large number of linear equations similar to eq. (3) can be created.

Surface wave measurements in eqs (1)–(3) are given as travel times for regional data with rays completely inside the study region (‘regional’ in Fig. 3a). For teleseismic data with rays partly inside the study region (‘teleseismic’ in Fig. 3a), the measurements are given as traveltime differences between each pair of closely located paths (Feng & An 2010). Since part of a teleseismic path is outside the study region, teleseismic data contain larger uncertainties than regional data. However, teleseismic data give an important contribution to long-period measurements (>100 s), which are helpful to constrain lithospheric-mantle structures. A detailed description on how to choose closely located paths for teleseismic data and on how to include them in the inversion for 3-D Vs is given in Feng & An (2010). For both regional and teleseismic cases, the data to be directly fitted in the inversion are the differences between observations (travel times or traveltime differences) and predictions from a reference model.

3.2 Relationship between receiver functions and 3-D Vs

Figure 4.

Sketch of ray paths for two receiver functions beneath one station (RF1 and RF2). P = P wave, Ps = Ps wave, PxXs = one of the multiples (PsPs, PpPs or PpSs). Piercing points of Ps for RF1 and RF2 are, respectively, labelled as P₁ and P₂, which may locate in different cells (ith and (i+1)th). Piercing point of multiples for RF1 is labelled as P_1x. The shaded crustal area is sampled by P, Ps and multiples of RF1. RF1 is used to constrain structures beneath the ith cell where P₁ is located.

Open in new tabDownload slide

A RF of a seismic station is primarily related to the seismic-velocity structure along the ray path to the station. If the rays from various backazimuths and/or incident angles to one seismic station pass through the same cell, the 3-D approach of eq. (6) and the 1-D approach of eq. (5) will give the same results. However, if RF traces for one seismic station pass through different horizontal cells, the two approaches may give different results, as illustrated in Fig. 4. The piercing points at the Moho for RF1 and RF2 (P1, P2 in Fig. 4) recorded by the same seismic station fall into two different cells, the ith and the (i+1)th cell, respectively. Then, RF1 constrains Vs beneath the ith cell, while RF2 constrains the Vs of the (i+1)th cell. In this case, the 3-D approach will be superior to the 1-D approach since it is able to reveal lateral variation of the Moho discontinuity. Note that this superiority gradually emerges as the data distribution improves and the cell size becomes smaller.

In the case that many RF rays cross through the same cell, the partial derivatives at each delay time t of different receiver functions in eq. (6) will be similar since RF traces within a small area (in one cell) are similar and the reference model for one cell is the same. To decrease the redundancy of similar equations from similar rays that cross through the same cell, we superpose such equations into one equation. The superposition is different from directly stacking RFs, since receiver functions of the same station with different incident angles have different equations.

The idea to directly invert for the 3-D velocity model from RF using eq. (7) is actually based on partial derivatives of the 1-D velocity model beneath each cell, so the matrix G_rf in eq. (7) is sparse and ill-posed. In any case, the 3-D approach adopted here is more scalable than the 1-D approach not only because it enables the identification of lateral variation in the Moho beneath one station if different RF rays of one station pass through different cells, but also because 3-D a priori constraints (e.g. 3-D model smoothness/flatness) can be easily considered in the inversion for a 3-D model. These constraints intend not only to produce reasonably distributed structures for neighbouring cells, but also suppress the effects of noisy phases in the RF.

3.3 Joint inversion of RF and SWD directly for 3-D Vs

At this point, we can mathematically resolve Vs perturbations Δβ from the linearized eq. (9) by the program LSQR (Paige & Saunders 1982), where both λ_rf and λ_c are determined by trial and error. Then the 3-D Vs solution β is β₀ +Δβ. The inverted model β can offer a better data fitting than β₀. The reference model at the first iteration, β₀, is also called the starting model or initial reference model. By taking β as the new reference model and then repeating the inverse processes for another iteration, a new solution can be obtained. The solution at subsequent iterations should be better than that at former iterations. In a linearized inversion, multiple iterations are normally required for solution convergence.

For a practical study, several key points need to be emphasized:

1. The 3-D inversion requires an initial reference (or starting) model as real as possible to decrease the non-linearity and trade-off between Vp/Vs and Vs in the RF inversion, and to give well-customized a priori flatness constraints (weighted by the velocity contrast between neighbouring cells/layers).

2. Eq. (7) for the 3-D inversion of RF is not designed to replace 1-D RF inversion, but as an adaptation for a 3-D joint inversion of SWD and RF. Independent 3-D inversion of RF is recommended only for a study region covered by very densely distributed seismic stations.

3. SWD and RF have different spatial resolving power. SWD of long wavelength character produces longer spatial-resolution length in the model than RF of short wavelength character. Therefore, the resolution of a jointly inverted model from SWD and RF strongly varies in space.

4. The effects of radial and azimuthal anisotropies are not considered in the 3-D joint inversion of SWD and RF presented here. The inverted 3-D Vs model is thus radially and azimuthally isotropic.

5. Both reference and resulting models only contain Vs but not crustal thicknesses. Crustal thicknesses are estimated from the inverted Vs structures by the approach of An et al. (2015) to be introduced in the next section.

4 RESULTS

The study region is discretized into a 3-D grid, with a 0.5° spacing in both latitude and longitude, and a variable layer thickness from the surface to 350 km depth. The layer thickness is 2 km for depths of ≤80 km, 5 km for depths of 80–120 km, 10 km for 120–200 km, 20 km for 200–300 km and 50 km for 300–350 km. As mentioned above, an initial reference (or starting) model with confident crustal structure (Vp/Vs, Vs and thickness) is required. The crustal structure is firstly interpolated from CRUST1.0 (Laske et al. 2013), and mantle structure is from IASP91 (Kennett & Engdahl 1991). The crustal thickness at 729 points from previous RF studies (Li et al. 2014b; He et al. 2014; Feng et al. 2017) is used to obtain the Moho depth beneath the related cells in the starting model. Fig. S3(a) shows crustal thicknesses from CRUST1.0 and from previous RF studies that are used to construct the starting model. As most previous studies did not provide Vp/Vs ratios, we keep them the same as in CRUST1.0. Therefore, the starting model already contains information about Moho depth constrained by receiver functions. As Q strongly influences the amplitudes of RF, by synthesizing previous experiences (Kennett et al. 1995; Dziewonski & Anderson 1981; Durek & Ekström 1996; Dalton et al. 2008), we set 300 and 374 for crustal Qs and Qp, and 200 and 352 for mantle Qs and Qp, respectively, in the starting model. We denote the starting model as START. Given the starting model and the prepared RF and SWD data, we can obtain a 3-D Vs model by iteratively inverting the linearized eq. (9).

A non-linear inverse problem solved by linearized approximation, like the joint inversion here, requires multiple iterations. However, tests of 1-D linearized joint inversion of SWD and RF (Feng et al. 2017) show that a two-step inversion, using the output of the SWD inversion as the reference model for an inversion using RF, with a few iterations is much better than a direct SWD-RF joint inversion using many iterations. Thus, we set a very small λ_rf in eq. (9) to give a very small weight to the RF data compared to the SWD data in the 1st inverse iteration, to produce a Vs model mainly constrained by SWD data. The inverted 3-D model after the 1st inverse iteration is called SURF. Synthetic tests of SWD inversions (Feng & An 2010) showed that the solution of the 1st inverse iteration is similar enough to the true synthetic model and can fit the observations reasonably well. Thus, SURF can represent the model inverted only from SWD. The 2nd inverse iteration is run by taking SURF as the reference model and by giving a balanced weight between SWD and RF data. More iterations are run, but produce trivial model variations relative to the solution after the 2nd iteration, as the starting model contains the information on Moho depths from previous RF studies. We therefore take the solution after the 2nd iteration as the final model and call it JOINT.

Before the formal inversion for each iteration, a pilot inversion is run to discard bad observations with a misfit greater than 2σ (two times the standard deviation) of the total misfit distribution, respectively, for the RF data and for the four sources of SWD data (Fig. 3a). The formal inversion is then run with the remaining good data after about 3.5–14 per cent of the observations are discarded.

4.1 Model appraisal

Fitting and misfit (standard deviation, SD) between observed and predicted data are goodness indicators of the 3-D models. Figs S4 and S5 show examples of SWD and RF fits between observed and predicted data, respectively. Generally, surface waves at short periods (<60 s) have a better fitting than at longer periods (Fig. S4), since the quantity of observations at periods of <60 s (Fig. 3) is larger, and then they have a higher contribution in the inversion than the long-period data. When the START model contains large uncertainties in the Moho depth or in the velocities around the Moho (e.g. station HWS in Fig. S5), the structures around the Moho are obviously improved in the JOINT inversion with inclusion of RF. Otherwise, if the START model is close to being realistic, the JOINT inversion does not significantly improve the structures (e.g. station DXI in Fig. S5).

We separately calculated misfit for the four sources of SWD data (Fig. 3) and found that misfit for the teleseismic event data is largest, and misfit for the regional noise data is smallest, as expected. Here we only show the mean misfit for all the four sources of SWD data to evaluate the models. Fig. 5 gives misfits of group velocities for the starting model (START), the inverted model SURF after the 1st inverse iteration, and the inverted model JOINT after the 2nd inverse iteration. The JOINT and SURF models give similar misfits, but about 50 PER CENT reduction of misfit relative to the START model (Fig. 5), confirming that surface wave observations can be significantly fitted after only one inverse iteration (Feng & An 2010). As a lot of a priori crustal thicknesses from previous RF studies were included in our START model, inclusion of RF in the JOINT inversion did not strongly improve structures down to and around the Moho for most of the stations. Thus, there is no remarkable effect on SWD fitting at periods around 40–80 s which are primarily sensitive to structures around the Moho (Fig. 5).

Figure 5.

Misfit between observed and predicted group velocities for the START (blue line), SURF (green line) and JOINT models (red line).

Open in new tabDownload slide

However, the JOINT model gives a certain misfit variation at periods of >110 s, such as a bigger misfit at periods of 110–145 s but a smaller misfit at periods of >145 s, compared to the SURF model. Such strong misfit variation occurs at longer periods because most of the long-period data are contributed to by the ‘teleseismic event’ data containing the largest uncertainty compared to the other three sources of data (Fig. 3a). Besides, long-period observations (<∼1000) are dramatically fewer than short-period data (>∼10 000, Fig. 3a), and thus have a much smaller relative weight in the joint inversion. Consequently, the fit or misfit for long-period data is easily affected by any adjustment of the relative weights among the SWD, RF and flattening. Although converted waves and multiples from the Moho dominate the RF waveforms and mainly constrain structures above and around the Moho, side effects of wave trains within the time-window of –2 to 50 s on deep structures emerge where the contribution of SWD data is weak. In any case, the JOINT model gives about 50 PER CENT reduction of misfit relative to the starting model.

Fig. 6 shows the misfit of RF waveforms for the START, SURF and JOINT models. The START model gives the widest misfit distribution with a mean misfit of 0.027 and a maximum value of 0.049. The SURF model slightly worsens the mean misfit to 0.028 (>0.027 for START) and slightly narrows the distribution to a maximum value of 0.045. It is not surprising that the SURF model is slightly worse than START for RF, since the Moho depths in START were obtained using results of previous RF studies and surface waves are not sensitive to the sharp velocity contrast around the Moho. However, the JOINT model, with a mean misfit of 0.024, produces 11 per cent of misfit reduction compared to the mean misfit (0.027) of START, and narrows the distribution to a much smaller maximum value of 0.04 (<0.049, Fig. 6).

Figure 6.

Misfit between observed and predicted RF waveforms for the START, SURF and JOINT models. The value in the top right of each panel is the mean value of the misfit.

Open in new tabDownload slide

In total, compared to the START model, the JOINT model produces about 50 PER CENT of misfit reduction for SWD data, and 11 PER CENT of misfit reduction for RF data (Figs 5 and 6). The JOINT model fits the group velocities at most periods as well as the SURF model and fits the RF waveforms best. Therefore, a 3-D model simultaneously constrained by SWD and RF is superior to a model purely constrained by one type of observations, as expected.

Checkerboard tests are performed to assess the ability of the joint inversion to resolve structural details of the central–east Tibetan plateau, using the inversion kernels of the last joint inverse iteration. Alternating checkers with ±0.07 km s^–1 of Vs perturbations were set as structures in the input 3-D model space. Since surface waves at different periods have different sensitivity and lateral resolution at different depths, we performed synthetic tests for four different 3-D checker sizes (Fig. 7): 2° × 2° × 16−50 km in x, y and z directions, respectively; 3° × 3° × 20–60 km; 4° × 4° × 20−70 km and 5° × 5° × 30−80 km. Random noise with a standard normal distribution was added to the synthetic data, and the joint inversion was then performed to determine how well the checkers were recovered. Fig. 7 shows a series of horizontal and vertical slices from the recovered models. For most of the study region, the model recovered 2° × 2° × 16−30 km checkers down to ∼50 km, 3° × 3° × 20−40 km checkers down to ∼80 km, 4° × 4° × 20−50 km checkers down to ∼150 km, and 5° × 5° × 30−80 km checkers down to ∼200 km depth. In summary, the resolution length varies from ∼200 to ∼500 km laterally and from ∼15 to ∼80 km vertically. The resolution decreases with depth and is small at the margins where ray path coverage is low.

Figure 7.

Horizontal and vertical Vs slices of the checkerboard tests. The sizes of the 3-D checkers are as follows from top to bottom: 2° × 2° × 16−50 km (in x, y and z direction, respectively), 3° × 3° × 20−60 km, 4° × 4° × 20−70 km and 5° × 5° × 30−80 km. Horizontal slices at four depths are shown in the left-hand column, from top to bottom: 10 km; 50 km; 120 and 150 km. Two vertical transects through each retrieved model are shown in the right-hand column. The locations of the vertical transects are marked on the horizontal slices.

Open in new tabDownload slide

4.2 Crustal thickness

Figure 8.

Crustal thickness extracted from our jointly inverted 3-D Vs model. Abbreviations are as follows: AS = Ailaoshan Suture, ATF = Altyn Tagh Fault, BM = Burma Microplate, BNS = Bangong–Nujiang suture, CAOB = Central Asian Orogenic Belt, IT = Indochina Terrane, JRS = Jinsha River Suture, KF = Kunlun Fault, LMSF = Longmenshan Fault, LS = Lanchangjiang Suture, MFT = Main Frontal Thrust, NCC = North China Craton, NQT = North Qaidam Thrust, NT = Naga Thrust, QSFT = Qilian Shan Frontal Thrust, RRF = Red River Fault, SCC = South China Craton, SF = Sagaing Fault, SGT = Songpan–Ganzi Terrane, YZS = Yarlung Zangbo Suture.

Open in new tabDownload slide

The crustal thickness map (Fig. 8) highlights a series of notable features that correspond to tectonic domains, with strong lateral variations across the major boundaries shaping the Tibetan plateau. These features are consistent with previous studies, including deep seismic soundings (Zeng et al. 1995; Teng et al. 2003) and passive seismic studies in the south (Acton et al. 2011), northeast (Vergne et al. 2002; Karplus et al. 2011; Ye et al. 2015), east (Zhang et al. 2009; Qian et al. 2018; Liu et al. 2014) and southeast (Wang et al. 2018a) of the plateau. For example, the thickest crust of 70–75 km is observed between the northern Lhasa Terrane and the Kunlun Fault (KF). To the south of the plateau, there is an abrupt crustal thinning from ∼60 km in the central Himalaya Terrane to ∼35–40 km across the Main Frontal Thrust (MFT). To the east of the eastern Himalayan syntaxis, the crust thins gradually from ∼65 km in the Hengduan Mountains to ∼30–35 km in the Indochina (or Indosinian) Terrane. To the east of the Tibetan plateau, there is an abrupt crustal thinning from the plateau (∼60 km) to the Sichuan Basin (∼40 km) across the eastern basin-mountain transition, the Longmen Shan Fault (LMSF). However, crustal thinning becomes more gradual to the north, with ∼10 km of crustal thinning across the Kunlun Fault to the Kunlun–Qilian Orogen (50–55 km) and a further ∼10 km of crustal thinning to the North China Craton (40–45 km). A further description of crustal thicknesses (or depth of the Moho discontinuity) will be given in the next section when showing Vs along vertical profiles.

4.3 3-D Vs

Fig. 9 shows 2-D horizontal and vertical slices of Vs from the final 3-D tomographic model JOINT. The Vs of the 3-D model are provided in Table S2. Most of the Tibetan plateau possesses a crustal thickness >60 km (Fig. 8), whereas the regions surrounding the plateau possess an average crustal thickness of ∼40 km. The velocities at 45 km depth (Fig. 9a) are thus representative either of the lower crustal structure inside the Tibetan plateau or of the upper mantle structure outside the plateau. The distribution of Vs at 45 km depth correlates well with the intensity of crustal lateral deformation (GPS movements, Fig. 9a). Strong low velocities occupy the whole plateau area south of the KF, where the largest GPS movements have been consistently measured [arrows in Fig. 9a, (Zhao et al. 2015)]. Weaker low velocities are distributed in most of the Kunlun–Qilian Orogen north of the KF, except the Qaidam Basin, where intermediate GPS movements have been found. High velocities are distributed in the southern, eastern and northern surrounding regions, where the smallest GPS movements have been found. The above Vs distribution pattern at 45 km depth which is obtained from RF and Rayleigh wave group velocities is quite similar to the Rayleigh phase velocity map of Yang et al. (2010) at 35 s period, which is most sensitive to Vs at 45 km depth.

Figure 9.

Horizontal slices of Vs at 45, 100 and 150 km depths (a–c) and vertical transects along profiles a–a’, b–b’, c–c’, d–d’ and e–e’ (d–h). In (a), arrows are GPS measurements from Zhao et al. (2015). In (b), magenta and purple lines are iso-velocity contours at 4.37 km s^–1 marking the limits of the Z-shaped low-velocity belt at 100 km depth. Black dashed line around 90°E, 32°N denotes a NE–SW trending fault from recent geological maps (Ren et al. 2013; Wang et al. 2013). In (c), black lines labelled with letters at the ends are locations of the vertical transects shown in (d–h). Other symbols in (a–c) are the same as in Fig. 1. In the vertical transects (d–h), the grey shading on the top is exaggerated topography, together with red bars marking major sutures or faults. Small circles are historic earthquakes from 1960–2016 from the EHB catalogue (Engdahl et al. 1998). Magenta and purple crosses are the intersecting points with the magenta and purple iso-velocity contours of 4.37 km s^–1 in (b), respectively. Crosses and short black vertical bars with labels at the depths of 100–150 km mark possible plate-boundary positions discussed in the text. Dashes mark the Moho from Fig. 8. Strong Moho variations are highlighted by ‘offset’ labels. Dotted lines are sketched to facilitate correlating surface sutures/faults with Moho variations, and crustal and mantle structures. The 5 transects share the same colour bar below transect a–a’. Abbreviations are as follows: ALP = Alashan Platform, AP = Asian Plate, AS = Ailaoshan Suture, ATF = Altyn Tagh Fault, BM = Burma Microplate, BNS = Bangong–Nujiang suture, CAOB = Central Asian Orogenic Belt, HT = Himalaya Terrane, IP = Indian Plate, IT = Indochina Terrane, JRS = Jinsha River Suture, KF = Kunlun Fault, KQO = Kunlun – Qilian Orogen, LMSF = Longmenshan Fault, LS = Lanchangjiang Suture, LT = Lhasa Terrane, MFT = Main Frontal Thrust, NCC = North China Craton, NQT = North Qaidam Thrust, NT = Naga Thrust, OB = Ordos Basin, QB = Qaidam Basin, QSFT = Qilian Shan Frontal Thrust, QT = Qiangtang Terrane, RRF = Red River Fault, SB = Sichuan Basin, SCC = South China Craton, SF = Sagaing Fault, ST = Sibumasu Terrane, SGT = Songpan–Ganzi Terrane, TB = Tarim Basin, TP = Tibetan Plate, YZS = Yarlung Zangbo Suture.

Open in new tabDownload slide

Velocities at 100 km depth (Fig. 9b) reflect structures of the lithospheric upper mantle, since the lithosphere beneath most of the study region is thicker than ∼100 km (An & Shi 2006). A notable feature at 100 km depth in the whole study region is a Z-shaped low-velocity belt defined by a Vs value of 4.37 km s^–1 (magenta and purple lines in Fig. 9b) separating two high-velocity areas. The eastern boundary of our Z-shaped low-velocity belt (purple lines) is nearly the same boundary which separates lithospheric low-to-high Vp anomalies around the Eastern Tibetan plateau in a recent body wave tomography (Zhang et al. 2018). The low-velocity belt starts from the Kunlun–Qilian Orogen in the north, passes through the Songpan–Ganzi terrane (SGT) and Hengduan mountains in the east, and ends in the Sibumasu and Indochina Terranes in the south. The southwestern high-velocity area includes the east Qiangtang, Lhasa and Himalaya Terranes around the eastern Himalayan syntaxis, which are directly affected by continental collision/subduction of the Indian Plate. The northern and eastern high-velocity areas cover the stable blocks of the Asian Plate, including the Tarim Craton, North China Craton (NCC) and South China Craton (SCC).

At 150 km depth (Fig. 9c), high velocities generally correspond to the lithosphere, and low velocities to the asthenosphere or weak lithosphere. At this depth, the Z-shaped low-velocity belt disappears, and most of the Tibetan plateau is dominated by high velocity, including the northern Qiangtang terrane along the JRS, the eastern Songpan–Ganzi Terrane, regions around the eastern Himalayan syntaxis and the Hengduan mountains, implying that the lithosphere may be thicker than ∼150 km in these areas. However, most of the areas surrounding the plateau, including the Central Asian Orogenic belt (CAOB), Qaidam basin, NCC, Qinling Orogen, northern SCC, and Indochina terrane (IT), are all characterized by low velocities, implying that the lithosphere in these regions may be weaker or thinner than ∼150 km. Weak low velocities roughly around the longitude of 90°E in the southern Tibetan plateau can be found at 150 km depth. Similar low-velocity anomalies were interpreted as tearing of Indian mantle lithosphere (e.g. Li & Song 2018).

Vertical slices of Vs are given for 5 profiles through different regions of the central–east Tibetan plateau (profiles from a–a’ to e–e’ in Figs 9d–h). Profile a–a’ is a south–north striking transect from the Indian Plate in the south to the Altyn Tagh fault in the north, in the same direction as the current Indian–Eurasian relative movement (Kreemer et al. 2014, marked as a big white arrow in Fig. 9c). The most prominent feature on this profile are the three segments of high upper-mantle velocities separated by two dotted lines: a northward deepening one from India to the BNS in the south, a flat one between the BNS and the KF in the middle, and a northward deepening one beneath the Tarim and Qaidam basins in the north. A notable feature in the crust is that the Qaidam basin (QB) has higher middle and lower crustal Vs than the SGT south of it and the Tarim basin (TB) north of it. The Moho quickly deepens from ∼40 km beneath northern India to ∼75 km beneath the central Himalaya in a distance as short as ∼220 km and remains at ∼70–75 km beneath the Himalaya, Lhasa, Qiangtang and Songpan–Ganzi terranes. It suddenly shallows to ∼50 km beneath the Qaidam and Tarim basins north of the KF. An obvious Moho offset occurs beneath the KF.

Profile b–b’ (Fig. 9e) is a SW-NE striking transect from the eastern Himalayan syntaxis in the southwest to the Alashan platform (ALP) in the northeast (westernmost part of the NCC). Again, there are three segments of high upper-mantle velocities: the northward deepening one south of the Jinsha River Suture (JRS); the flat one between the JRS and the North Qaidam Thrust (NQT); a weakly southward deepening one beneath the Kunlun–Qilian orogen and the Alashan platform. Strong crustal-velocity variations are found beneath the Kunlun–Qilian orogen (KQO) between the KF and NQT. Moho depths gradually deepen from India to the BNS in the south, shallow in a stepwise-like manner from the JRS to the KF in the middle, and then gradually shallow beneath the Kunlun–Qilian orogen and the Alashan platform in the north. The Moho step is ∼15–20 km and sharp beneath the KF, and ∼10 km and relatively gentle beneath the JRS.

Profile c–c’, located east of profile b–b’, is also a SW-NE striking transect from the eastern Himalaya syntaxis in the southwest to the Ordos basin (NCC) in the northeast. It shows quite a similar appearance of upper-mantle velocities to profile b–b’: one segment of northeastward deepening high velocities south of the JRS; another flat one between the JRS and the KF or NQT; a weakly southward deepening one beneath the KQO and the NCC. The strongest crustal velocity variation occurs beneath the NQT. The sharpest Moho depth variation occurs beneath the KF.

Profile d–d’ is a NW-SE trending transect crossing through the eastern SGT and the Sichuan basin. It is roughly perpendicular to profile c–c’ and to the LMSF between the SGT and the Sichuan basin. Three segments of high upper-mantle velocities extend beneath the whole profile. The one beneath the Sichuan basin (SCC) has the strongest intensity. Crustal velocities are laterally homogeneous in general. Moho depth varies abruptly beneath the LMSF, beneath which the region between two of the high upper-mantle velocity segments also occurs.

Profile e–e’ is a W–E trending transect that mainly crosses through regions east of the eastern Himalayan syntaxis. Along the profile, there are only two high-velocity segments. An eastward deepening high upper-mantle velocity block is observed beneath the Burma microplate (BM) and the Sibumasu Terrane (ST, Fig. 9h). The other eastward thickening high-velocity block is observed beneath the SCC. Major crustal velocity anomalies locate around the Sagaing Fault. Though elevation strongly varies, and several large faults are exposed at the surface along the profile, Moho depths keep a stable value of ∼35–40 km with no sudden variations.

Besides the prominent feature of three high-velocity segments in the mantle, another common feature on the vertical transects is that the middle high-velocity segment along profiles a–a’, b–b’, c–c’ and d–d’ and the transition zone between the two high-velocity segments along profile e–e’ are almost always overlain by about 50 km of uppermost mantle with slower velocities (labelled as ‘soft’ in Figs 9d–h), correspondent to the Z-shaped low-velocity belt at 100 km depth in Fig. 9(b). Especially in the profiles b–b’ and c–c’ (Figs 9e and f), the base of the Z-shaped low-velocity belt looks like a discontinuity at ∼110 km. A weak discontinuity was detected in these areas by Feng et al. (2011). The Z-shaped low-velocity belt looks like an upper-mantle soft channel lying between the Moho and the discontinuity. Previous studies reported widespread low Vs in the middle-to-lower crust along the periphery of the Tibetan plateau (e.g. Yang et al. 2012). In our 3-D Vs model, these low velocities likely extend down into the uppermost mantle.

5 DISCUSSION

The Vs and crustal thicknesses (or depths of the Moho discontinuity) obtained from our joint tomography of SWD and RF, as described above, provide a wealth of structural information for the central–east Tibetan plateau in a 3-D perspective view, from which new understanding on regional dynamics can be extracted.

5.1 Deformation zonation of the Tibetan plateau

The Tibetan plateau has been elevated by compression after the Indian–Eurasian collision occurred at 55–60 Ma (Yin & Harrison 2000). Thus, there is consensus that the large area bounded by the Main Frontal Thrust in the south, the LMSF in the east, and the Altyn Tagh Fault (ATF) and Qilian Shan Frontal Thrust (QSFT) in the north, with an average elevation of >∼3000 m, must have been impacted by the Indian–Eurasian collision. As to what extent the different parts of the plateau have been involved in the post-collisional deformation is still under debate.

The plateau regions bounded by the KF in the north, the LMSF in the east, and the Yunnan-Guizhou plateau in the southeast, except for the triangular part of the Songpan–Ganzi Terrane (SGT) east of 102°E, are well shaped by thick crust of >60 km (Fig. 8). The Kunlun–Qilian Orogen (KQO), except for the central Qaidam and the part roughly east of 102°E, is dominated by a smaller crustal thickness of ∼55 km. The regions surrounding the plateau have a much thinner crust of <45 km. The Vs at 45 km depth (Fig. 9a) show the same pattern. Regions bounded by the KF, LMSF and the Yunnan-Guizhou plateau are covered by very low velocities of <3.6 km s^–1, while the Kunlun–Qilian Orogen has higher velocities of ∼3.7 km s^–1. Exceptions again appear in the central Qaidam and regions east of 102°E where velocities are as high as 3.8–3.9 km s^–1. The regions surrounding the plateau have much higher average velocity of >4.1 km s^–1. This zonal distribution of the crustal structures suggests that the region of the plateau south of the KF, the Kunlun–Qilian Orogen north of the KF, and the regions surrounding the plateau have experienced deformations of decreasing intensity.

Within the Kunlun–Qilian Orogen, crustal thickness and Vs at 45 km depth in the central Qaidam and regions east of 102°E are different from the other regions of the orogen. Deep structures across the orogen (profiles a–a’, b–b’ and c–c’ in Figs 9d–f) show a certain correlation to the shallow structures. On profile a–a’ that crosses through the western Kunlun–Qilian orogen (including Qaidam), a high upper-mantle velocity block beneath the Tarim and Qaidam basins seems overthrusting southwards to the northern Tibetan plateau and forms a dome-shaped Moho offset beneath southern Qaidam. Numerical modelling results (Sun & Liu 2018) showed that different rheological contrasts between the plateau and surrounding blocks caused different ways of lateral growth of the plateau. Therefore, it is not surprising that the Qaidam basin, characterized by higher crustal and upper-mantle velocities, and thus by higher viscosities, behaves like a rigid block when the plateau tries to expand across the KF. The ‘dome-shaped’ Moho offset and the ‘overthrusting’ high upper-mantle velocities both imply that the western Kunlun–Qilian Orogen, except for the rigid Qaidam basin itself, has experienced strong deformation. For the middle Kunlun–Qilian Orogen east of the Qaidam and west of 102°E (profile b–b’), a Moho offset beneath the KF and gently southward dipping high upper-mantle velocities below the Kunlun–Qilian Orogen and Alashan platform are observed, as previously imaged by Kind et al. (2002) and Yue et al. (2012). The Moho offset is not as large as the ‘dome-shaped’ one on profile a–a’ and the lateral variation of crustal velocity across the KF is not as sharp as that on profile a–a’, implying that the middle Kunlun–Qilian Orogen has experienced weaker deformation than the western orogen. The eastern orogen roughly east of 102°E imaged no high upper-mantle velocities but a gap of normal velocities between the SGT and the Ordos (profile c–c’). The Moho offset beneath the east KF on profile c–c’ becomes even smaller and the lateral variation of crustal velocity across the KF becomes more gradual, implying that the eastern Kunlun–Qilian Orogen roughly east of 102°E has experienced the weakest deformation. Therefore, our model shows that the deformation intensity of the Kunlun–Qilian Orogen is gradually reduced from west to east. Other controlled and passive source seismic studies have also imaged a decrease in the Moho offset across the KF from west to east (e.g. Zhang et al. 2011).

A similar structural variation is observed across the eastern and southeastern boundaries of the Tibetan plateau. For the east boundary (profile d–d’), the image of high upper-mantle velocity beneath the Sichuan basin (SCC) and the large Moho offset around the LMSF is quite similar to the case of the Qaidam basin and the KF in northern Tibet shown on profile a–a’, implying that eastward expansion of the plateau has been resisted by the rigid Sichuan basin and very strong deformation has occurred in regions around the LMSF. For the southeast boundary (profile e–e’), there exists a gap between the rigid SCC with high upper-mantle velocity and the eastward subducting Indian slab beneath the Burma microplate, and no Moho offset is observed around the southeast boundary of the plateau (the Lanchangjiang Suture, LS), implying weaker deformation has occurred here.

In summary, the highly elevated region of the plateau with very thick crust has experienced the strongest deformation since the Indian–Eurasian collision. The boundary regions of the plateau have undergone deformation to different degrees. The eastern boundary around the LMSF is maturely and strongly deformed, while the southeastern boundary around the LS is newly and weakly deformed. The Kunlun–Qilian Orogen is intermediately deformed and has decreasing deformation intensity from west to east. The northeast and southeast, where the rigid NCC or SCC are distant from the boundaries of the plateau and bounded by wider and softer mantle materials, might be two potential expansion directions of the Tibetan plateau in the future.

5.2 Indian, Tibetan and Asian lithospheres beneath the Tibetan plateau

The geometry of the Indian slab and the location of its leading edge beneath the Tibetan plateau are important to understand the tectonic evolution of the plateau. Receiver function analysis (Kosarev et al. 1999) along the INDEPTH profile observed a northward dipping structure from the Moho 50 km north of the YZS to 200 km beneath the BNS, which was interpreted as the northward subducting Indian lithosphere. Later, body wave tomography using INDEPTH data (Tilmann et al. 2003) showed that the Indian Plate underthrusted the plateau to as far as the BNS in central Tibet, and that the Indian Plate might then sink subvertically to at least 400 km depth. However, regional body wave tomography (Li et al. 2008) imaged no high velocities down to about 200 km depth between the YZS and the BNS and suggested that Indian lithosphere underlay only the southwestern part of the plateau and that the central and northeastern part was underlain by Asian lithosphere. A S-wave receiver function study (Zhao et al. 2011) observed Sp waves converted from the lithosphere–asthenosphere boundary (LAB) at various depths beneath the Tibetan plateau and correlated them, respectively, with the Indian, Tibetan and Asian plates. They followed the interpretation of Tilmann et al. (2003) that the Indian slab stopped or sank beneath the BNS, but they further proposed that the Indian slab stopped/sank there because it confronted the southward subducting Asian Plate from the northeast. Several surface wave models also show the Indian Plate underthrusting the plateau as far north as to the southern half of the Qiangtang terrane (e.g. Friederich 2003; Feng et al. 2011). In contrast, the surface wave model of Priestley et al. (2006) and P-wave model of Zhou & Murphy (2005) show that all of Tibet is underlain by high velocity Indian lithosphere. The different results and interpretations have attracted widespread concerns on the geometry and location of the northern front of the subducting Indian slab, as well as the possible dynamics.

Our tomographic model imaged segmented high velocities from the Moho down to ∼200 km depth, possibly correlating with different lithospheric plates. There are three segments of high lithospheric velocities from south/southwest to north/northeast beneath the central–east Tibetan plateau. The boundaries of the high-velocity bodies may represent the leading or trailing edges of different lithospheric plates. Based on the aforementioned knowledge from previous studies, it is straightforward to correlate the northward deepening high-velocity segments in the south of profiles a–a’, b–b’ and c–c’ (Fig. 9) to the subducting Indian slab, and the southward shallowing or subducting high velocity segments in the north of profiles a–a’, b–b’ and c–c’ (Fig. 9) to the Asian Plate. We cannot tell whether a vertically sinking Indian slab exists because of our limited resolution for depths >200 km.

5.2.1 Indian slab

The northern edge of the Indian slab has been detected previously. Global (e.g. Lebedev & Van Der Hilst 2008) and Asian surface wave tomographic studies (e.g. Priestley et al. 2006) showed that the leading edge of the Indian slab occurred as far north as the KF. Our regional tomographic study using a dense coverage of surface waves imaged similar high upper-mantle velocities extending from the MFT northward to the KF (black bar labelled with IP/TP on profile c–c’ in Fig. 9f). However, an obvious drop occurs within the high velocities around the JRS (black bar labelled with IP in Fig. 9f) which seems to divide the high velocities into two segments. Furthermore, on profiles a–a’ and b–b’ in Figs 9(d) and (e), the southern high-velocity segment (south of the magenta cross labelled with IP) is more obviously disconnected from the middle high-velocity segment. Therefore, the north edge of the Indian slab, corresponding to the southern high-velocity segment, should not be far north of the JRS. Profiles a–a’ and b–b’ show that the edge of the Indian slab locates, respectively, around the BNS and the JRS (indicated by the magenta crosses labelled IP in Figs 9d and e, which coincide with the magenta boundary of the Z-shaped low-velocity belt at 100 km depth in Fig. 9b). To the west of our study region the edge of the Indian Plate may again locate north of the BNS (Hung et al. 2011).

In the southern part of our study area, a recent receiver function image down to ∼100 km depth (Zheng et al. 2020) showed that the eastward subducting Indian slab ends around the Sagaing Fault (SF), consistent with the position of the boundary of our Z-shaped belt at 100 km depth (magenta cross on profile e–e’ in Fig. 9h). However, our model with a deeper and wider coverage shows that the subducting Indian slab can be traced farther east of the SF (black bar labelled with IP in Fig. 9h) to a depth greater than ∼150 km.

Based on the boundaries of the Z-shaped low-velocity belt (magenta lines in Fig. 9b or magenta crosses, sometimes labelled with IP, in Figs 9d–h) and of the southern high-velocity segment on the vertical transects (black bars labelled with IP in Figs 9d–h), we give a rough delineation of the leading edge of the Indian slab (blue lines in Fig. 10a). As the magenta crosses locate at the end of the Indian slab on profiles a–a’ and b–b’ and in the middle of the Indian slab on profiles c–c’ and e–e’, the proposed leading edge of the Indian slab coincides with the magenta boundary of the Z-shaped low-velocity belt in the west but deviates in the east (blue lines in Fig. 10a).

Figure 10.

Illustration of lithospheric-mantle dynamics of the central–east Tibetan plateau. The base map in (a) is the relief slope calculated from the relief data of SRTM30 (Farr et al. 2007). Magenta and purple lines are the same boundaries of the Z-shaped low-velocity belt at the depth of 100 km in Fig. 9(b). Blue (labelled a–f) and dark green lines (labelled A–F) in (a) mark the boundaries of the lithospheric mantle of the Indian and Asian plates, according to mantle Vs in the horizontal and vertical slices in Fig. 9. Red triangle labelled as ‘NBS’ marks the Namche−Barwa syntaxis. Black arrows in (a) and circle with a point in (b) represent crustal escape directions from GPS in Fig. 9(a). CAOB = Central Asian Orogenic Belt, MFT = Main Frontal Thrust, YZS = Yarlung Zangbo Suture. Other symbols are the same as in Fig. 1.

Open in new tabDownload slide

The proposed leading edge of the Indian slab (blue lines in Fig. 10a) does not seem to correlate with topography or geology (e.g. elevation and tectonics in Fig. 1). However, it does show a good correlation with surface relief slopes as shown in Fig. 10(a). Surface relief slope is an indicator of uplift rate of the plateau: a steep slope corresponds to a fast uplift. The region above the proposed subducting Indian slab in the eastern plateau (bounded by the blue lines) has much larger surface relief slopes (i.e. faster surface uplift) than the western and northern plateau. A possible reason for the faster uplift is due to the Indian Plate sinking and thus causing a loss of denser mass in the uppermost parts of the mantle, which then provokes faster surface uplift, as an isostatic response. Meanwhile, light and buoyant crust will be strongly shortened by the Indian–Eurasian convergence, which will accelerate the surface uplift of the region, as illustrated in Fig. 10(b).

Different from most of the E–W striking faults or sutures exposed at the surface in the central Tibetan plateau, part of the Indian slab edge derived from our upper-mantle Vs (black dashed line in Fig. 9b, or the b–c segment of the blue lines in Fig. 10a, close to Bange and Ando) trends NE–SW. Geological maps (Ren et al. 2013; Wang et al. 2013) do show some NE-trending faults around the b–c edge segment. Besides, the central Himalayan and Lhasa terranes developed several N–S or NE–SW striking rift systems, while the eastern Plateau shows no rift development (Yin & Harrison 2000), implying that the proposed NE-trending segment of the Indian slab edge may mark a tectonic boundary separating the central and eastern Tibetan plateau. The high-velocity block east of the b–c edge segment (bounded by the magenta lines in Fig. 9b) accompanied by large surface relief slopes (bounded by the blue lines in Fig. 10a) looks like a giant syntaxis in the lithospheric mantle. According to Ren et al. (2013), the NE–SW faults around the b–c edge segment are characterized by sinistral strike-slip, implying that the Indian slab east of the edge segment tends to move farther north, consistent with the appearance of our proposed giant lithospheric-mantle syntaxis. The Eastern Himalayan or Namche-Barwa syntaxis (triangle labelled ‘NBS’ in Fig. 10a) locates around the centre of the giant syntaxis. Similar strongly curved segments appear both at the northeastern corner of the giant syntaxis (labelled ‘d’) and at the NBS on the Yarlung-Zangbo suture (Fig. 10a). This gives an impression that the former resulted by moving the latter northeastward in the approximate direction of the Indian–Eurasian convergence which, in turn, confirms that the giant high-velocity block delineated by the blue lines in Fig. 10(a) reflects the Indian slab.

The length of the Indian Plate along the three profiles a–a’, b–b’ and e–e’ is about 600 km. These three profiles have their southern or western ends relatively close together at about 26°N and about 70 km south of the Main Frontal Thrust in the case of profiles a–a’ and b–b’ and about 70 km west of the Naga Thrust in the case of profile e–e’. Along profile c–c’ whose southern end is at about 27.5°N, also about 70 km southwest of the Main Frontal Thrust, the length of the Indian Plate is 400–500 km. However, if one would largely ignore the protrusion defined by the Main Frontal Thrust and the Naga Thrust in this area and extend profile c–c’ south to about 26°N, then presumably the length of the Indian Plate would also be about 600 km. Thus, on all four profiles, a–a, b–b, c–c and e–e, the length of the Indian Plate measured from about 70 km south or west of the plate boundary at the surface, and ignoring the protrusion in NE India, is approximately 600 km.

5.2.2 Asian mantle lithosphere

According to Zhao et al. (2011), the middle segment of high velocity on profiles a–a’, b–b’ and c–c’ in Fig. 9 is a pre-existing Tibetan Plate underlain by subducting Asian lithosphere. Along profile b–b’, the middle segment (labelled as TP/AP in white in Fig. 9e) seems connected with the northern segment (the Asian Plate). Thus, the Asian Plate might at least extend as far south as the KF. However, the region between the middle and northern high-velocity segments shows a small drop of Vs, implying that the middle high-velocity segment is also possibly not associated with the Asian Plate. Furthermore, the middle high-velocity segment on profile a–a’ is obviously disconnected with the northern segment (the Asian Plate). Thus, the southern edge of the Asian Plate beneath the northern Tibetan plateau is possibly located around the NQT (purple crosses or black bar labelled with AP in Figs 9d and e, respectively). All profiles except for profile b–b’ show that the southern edge of the Asian mantle lithosphere occurs around the northern/eastern boundary of the Z-shaped low-velocity belt at 100 km depth (purple lines in Fig. 9b or purple crosses labelled with AP in Figs 9d–h).

Based on the northern/eastern boundary of the Z-shaped low-velocity belt at 100 km depth (purple lines in Fig. 9b or purple crosses labelled with AP in Figs 9d–h) and of the northern high-velocity segment on the vertical transects (black bar labelled with AP in Fig. 9e), we provide a rough delineation of the southern edge of the Asian Plate (dark green lines in Fig. 10a). As the purple crosses on all profiles except for profile b–b’ locate around the end of the northern high-velocity segment, the proposed edge of the Asian Plate (dark green lines in Fig. 10a) coincides with the purple boundary of the Z-shaped low-velocity belt, except for the C–D edge segment around profile b–b’ (dark green dashes).

5.2.3 Tibetan mantle lithospheres

Previous studies imaged two segments of high velocity beneath the Tibetan plateau, corresponding to the Indian slab and the Asian Plate, respectively. The mantle overlying the Indian slab and the Asian Plate is often taken as pre-existing Tibetan lithosphere (e.g. Priestley et al. 2006; Zhao et al. 2011; Agius & Lebedev 2013; Feng et al. 2011). From this point of view, the Z-shaped upper-mantle low velocity channel in our model should be the Tibetan mantle lithosphere. As the thickness of the Z-shaped channel with the overlying Tibetan crust is close to that of a general non-cratonic lithosphere (∼110 km), it is reasonable to take the mantle part as original Tibetan lithosphere (Feng et al. 2011).

According to Zhao et al. (2011), the middle segment of high velocity on profiles a–a’, b–b’ and c–c’ in Fig. 9 is a pre-existing Tibetan Plate underlain by subducting Asian lithosphere. According to Priestley et al. (2006), the middle high-velocity segment on profiles a–a’, b–b’ and c–c’ (Fig. 9) would belong to the Indian Plate. The cross-sections of Priestley et al. (2006) and Zhao et al. (2011) both lie between profiles a–a’ and b–b’ in the northern half of the plateau. However, besides the two major high-velocity segments from the south (corresponding to the Indian slab) and from the north (corresponding to the Asian Plate), our model imaged another smaller high-velocity segment in between which is overlain by a Z-shaped uppermost mantle low-velocity channel. Then, what could the smaller middle segment be?

When taking a closer look at the middle high-velocity segment on different vertical profiles, we find that the middle segment is disconnected with either of the other two segments on profile a–a’, almost connected with the Asian Plate segment from the north on profile b–b’, and almost connected with the Indian-slab segment from the south on profile c–c’. Thus, it is possible that the middle high-velocity segments in different parts of the Tibetan plateau do not belong to a unique plate but are just some isolated pieces. Vs maps at the depths of 100 and 150 km (Figs 9b and c) do show that the middle high-velocity segments on profiles a–a’, b–b’ and c–c’ are not connected with each other.

The Tibetan region has been built by sequential accretion through subduction and collision of several microcontinents and island arcs from north to south since the early Palaeozoic (Allégre et al. 1984; Dewey et al. 1988; Yin & Harrison 2000; Pan et al. 2012), including the Songpan-Ganzi terrane accreted to the Kunlun-Qilian Orogen, the Qiangtang terrane to the Songpan-Ganzi terrane, the Lhasa terrane to the Qiangtang terrane, etc. The Indian-slab segment in the mantle was easily recognized due to the known Indian–Eurasian collision. The Asian Plate segment was recognized due to the convergence between the Kunlun-Qilian Orogen and the Asian Plate. Thus, high-velocity blocks might also exist at depths below 100 km caused by collisions/subductions prior to the Indian–Eurasian collision. After the Indian–Eurasian collision, parts of such high-velocity blocks might be kept. Thus, the middle high-velocity segments in the Tibetan upper mantle may be remnants of some microplates associated with collisions or subductions prior to the Indian–Eurasian collision, during the accretion of the Tibetan plateau. The high-velocity remnants are not related to tectonic events associated with the Indian–Eurasian collision but are related to the formation of the Tibetan terranes, and can thus be termed as ‘pre-existing Tibetan lithosphere’ to discriminate with the ‘present Tibetan lithosphere’ (the Z-shaped channel with overlying crust).

In summary, the Tibetan mantle lithosphere might be composed of two parts (Fig. 10b): the Z-shaped low-velocity channel in the uppermost mantle and the high-velocity lithospheric remnants of the microplates or slabs during the accretion of the Tibetan terranes prior to the Indian–Eurasian collision, which correspond to soft/hot lithosphere and hard/cold lithospheric remnants, respectively.

There is another possibility for the middle high-velocity segment on profile b–b’. Since a less significant decrease in S-velocity occurs between the middle and northern high-velocity segments (black bar labelled AP in Fig. 9e), the middle segment along profile b–b’ may also belong to the Asian Plate (labelled as TP/AP in white in Fig. 9e). This would then mean that the northwestern high-velocity segment along profile d–d’, which ends close to where the middle high-velocity segment occurs along profile b–b’ south of the KF, may also be Asian lithosphere (labelled as TP/AP in white in Fig. 9g).

5.3 Upper mantle Z-shaped channel

As mentioned above, a Z-shaped upper-mantle low-velocity channel, which corresponds to the Tibetan lithospheric mantle, occurs roughly along the northern and eastern periphery of the Tibetan plateau (Fig. 9b). One possible reason for the low velocities in the Tibetan lithospheric mantle is that the uppermost lithospheric mantle could be heated up by increased radioactivity from the thickened Tibetan crust (McKenzie & Priestley 2008). This would then cause the mantle just below the Moho to be hotter and softer, and thus have slower velocities, than the lower colder and harder parts. Hotter and softer materials are more deformable than colder and harder materials. The Z-shaped low-velocity Tibetan lithospheric mantle may promote or have promoted the deformation of the Tibetan crust above it during the interplay of the different tectonic blocks (Fig. 10b).

Several oceanic or continental subductions/collisions occurred between the Tibetan terranes during their accretion (Allégre et al. 1984; Dewey et al. 1988; Yin & Harrison 2000; Pan et al. 2012). The Z-shaped channel mainly overlies the middle high-velocity segments. As the high-velocity segments are attributed to the remnants of the cold lithospheres in collisions/subductions during the Tibetan accretion, the overlying low velocity channel is also possibly related with old subduction. In a subduction system, water is released from the subducted slab due to break-down of hydrous minerals at high pressure and temperature. The derived water flux will go up and consequently influence the overlying lithosphere. The addition of water into the overlying mantle not only causes a marked lowering of seismic velocities (Karato & Jung 1998; Hirth & Kohlstedt 1996), but also has the potential to influence lithospheric strength dramatically (Jackson 2002), which will consequently promote the deformation of the overlying Tibetan crust.

6 CONCLUSIONS

A well-constrained 3-D lithospheric Vs model for the central–east Tibetan plateau was constructed by simultaneously inverting body wave receiver function waveforms and surface wave dispersion curves in a joint tomographic system. Crustal thicknesses of the region were then extracted from the Vs model. The crust is 55–75 km thick within the Tibetan plateau, and ∼40 km thick on average outside the plateau. The biggest Moho offset is found across the LMSF and the eastern KF just south of the Qaidam basin. A Moho offset widely exists along the eastern KF but becomes smaller from west to east.

In the uppermost mantle beneath the study region, a Z-shaped low-velocity channel, which was taken as the present Tibetan lithospheric mantle, and two large areas of high velocity corresponding to the Indian and Asian lithospheres are the major features. The Z-shaped low-velocity channel, reflecting easily deforming mantle material, may promote or have promoted the crustal deformation during the interplay among the Indian, Tibetan and Asian lithospheres. Segmented high upper-mantle velocities dominate most of the study region with different occurrences and different depth extents. We correlate the slightly northward deepening segment south of the BNS in central Tibet and south of the JRS in eastern Tibet with the subducting Indian slab, the southward shallowing segment beneath the Qaidam basin and the southward deepening segment beneath the eastern Kunlun–Qilian Orogen with the overthrusting or underthrusting Asian Plate, and the smaller flat segments in between mainly with the remnants of cold microplates/slabs associated with collisions/subductions prior to the Indian–Eurasian collision during the accretion of the Tibetan region.

The eastern boundary of the plateau (the LMSF), characterized by the largest Moho offset and the greatest upper-mantle velocity contrast, is proposed to be the most strongly deformed boundary of the plateau. The northern boundary (the KF), characterized by a weaker Moho offset, crustal-velocity variation and mantle-velocity contrast from west to east, has undergone an intermediate and variable degree of crustal deformation. The southeastern boundary of the plateau (LS), characterized by non-equilibrated crust and deep focus earthquakes, might be the most recently and weakly formed boundary. The northeast and southeast margins of the plateau are two potential regions for Tibetan plateau expansion in the future because they have a wide and thus more deformable upper-mantle low-velocity channel and both the NCC and SCC are distant from the boundaries of the plateau and thus there may be less resistance in these regions.

SUPPORTING INFORMATION

Figure S1. Epicentral distances and path coverage of Rayleigh wave measurements. (a) Average epicentral distance for different periods of Rayleigh wave group velocities for regional and teleseismic measurements, respectively. (b) Ray path density (colour shades) of surface wave measurements at 50 s period. The path density is defined as the number of group-velocity measurements in each 0.5° × 0.5° cell. Other symbols are the same as in Fig. 1.

Figure S2. Examples of P-wave receiver functions for seismic stations (a) DXI, (b) HWS, (c) DXX and (d) A07. Station locations are shown in the inset map of Fig. S5.

Figure S3. A priori crustal thicknesses and their differences with the result here. In (a), the coloured image shows crustal thicknesses from CRUST1.0 (Laske et al. 2013) and colours in circles represent thicknesses from previous RF studies (He et al. 2014; Li et al. 2014b; Feng et al. 2017). All those thicknesses are used to construct the START model. In (b), the differences between the results in Fig. 8 and the previous RF studies in (a) are shown.

Figure S4. Examples of data fits between observed (triangles) and predicted group-velocity dispersion curves from the START, SURF and JOINT models, respectively. The inset map shows the paths with stations (triangles) and epicentres (circles) for the four dispersion curves. The number in the upper left of each frame corresponds to the ray path indicated by the same number in the inset map. For path 2, the station is LSA (Lhasa), the M_s 6.3 earthquake occurred on 2004/08/10 in Yunnan, China. For path 4, the station is ENH, and the M_w 5.1 earthquake occurred on 2004/04/06 in Xizang (Tibet), China. The two stations belong to the network IC, and the data are available in IRIS. The data for paths 1 and 3 are from portable stations.

Figure S5. Examples of data fits between observed and predicted receiver functions from the START, SURF and JOINT models, respectively, for four cells. The cell ID and station code are given in the upper left of each frame and the locations are indicated in the inset map. Only receiver functions of the given station with Ps phases penetrating the given cell are shown. We can see that although the receiver functions mainly constrain structures down to the Moho in the JOINT inversion (red lines) using the model SURF as the reference model (green lines), change is still observable for Vs in layers far below the Moho because of the co-effects of crustal and mantle structures on the SWD.

Table S1. Crustal thicknesses of the study region.

Table S2. Vs of the study region.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

This research has been supported by the Natural Science Foundation of China (grants 41974051, 41574049 and 41174039). James Mechie received funding from the Deutsche Forschungsgemeinschaft and the Helmholtz-Zentrum Potsdam - Deutsches GeoForschungsZentrum GFZ, and the provision of instruments from the Geophysical Instrument Pool of the Deutsches GeoForschungsZentrum GFZ. We thank the team members of the INDEPTH-IV and Qilian Shan seismic arrays for making seismic data available. Waveform data from Chinese permanent seismic stations for this study were provided by the Data Management Centre of China National Seismic Network at Institute of Geophysics, China Earthquake Administration (SEISDMC, doi:10.11998/SeisDmc/SN). The facilities of the IRIS Data Management System, and specifically the IRIS Data Management Center were used for access to the waveform and metadata required in this study. We thank two anonymous reviewers for their constructive comments and suggestions. Most of the figures were generated using Generic Mapping Tools (GMT, Wessel & Smith 1991).

REFERENCES