Upper-Mantle Velocity Heterogeneity of Eastern Tibetan Plateau from Teleseismic P-Wave Tomography and Its Tectonic Implications

0. INTRODUCTION

The Tibetan Plateau, known as the roof of the Earth and the highest plateau or 'third pole' is not only a plateau with high elevation, but is also famous for its complex history and evolutionary theories. It is supposed to be formed more or less 50 million years ago (Yin and Harrison, 2000; Molnar and Tapponnier, 1975). Some famous statements about Tibetan Plateau state that the plateau was caused by the northward subduction of the Indian Plate lithosphere beneath the southern Tibetan Plateau (Replumaz et al., 2013; Zhao et al., 2011; Li et al., 2008; Tilmann et al., 2003) and southward subduction of Eurasian Plate lithosphere beneath northern Tibetan Plateau (Zhao et al., 2011; Kind et al., 2002; Willett and Beaumont, 1994). Also, the continental subduction of both plates contributes greatly to the plateau growth in the southern and northern Tibetan Plateau (Replumaz et al., 2013). Although good number of studies have been done and going on in this area, still many unrevealed facts and questions arise about this plateau. This means undoubtedly it is essential to keep continuing research around this region for a better understanding. There are many geodynamic models based on the evolution and structure of the Tibetan Plateau. Most of the models have tried to explain the reason for high topography but the mantle structure of the eastern Tibetan Plateau and its effects on topography is not revealed clearly (Lei et al., 2019). Maximum large earthquake epicenters of the eastern Tibetan Plateau are located at boundaries of the plateau which consists of heterogeneous low and high velocity anomalies (Lei and Zhao, 2016a). More and more precise research is needed to clear up the idea about upper mantle dynamics. Many studies have been done along the eastern margin of Tibetan Plateau to reveal the crust and upper mantle structure although differences among tomographic maps can be seen, which means all results are not giving same interpretations (Zheng et al., 2019; Guo and Chen, 2017; Bao et al., 2015).

A large number of studies have been carried out in the eastern Tibetan Plateau, but only a few studies focus on the lateral variations of underground structures or the velocity variation beneath the eastern Tibetan Plateau. To keep consistency with the tradition of finding unknown facts we attempt to focus on some unclear facts to bring them to light. For that reason, an attempt has been made to visualize the upper mantle heterogeneity of the eastern Tibetan Plateau.

1. TECTONIC SETTING AND RESEARCH BACKGROUND

In the last few decades many research works have been done to reveal the detailed tectonics of the Tibetan Plateau using seismic tomography, magneto telluric study, shear wave splitting analysis, receiver function, GPS etc. (Zhang H et al., 2021, 2018, 2017, 2012; Lei et al., 2019, 2014; Pei et al., 2019; Xu et al., 2019; Zhang F X et al., 2018; Ye et al., 2017; Zhang X Y et al., 2017; Peng et al., 2016; Lei and Zhao, 2016a, b, 2009; Li M K et al., 2016; Zhao J et al., 2014; Hu et al., 2013, 2011; Wei W et al., 2013; Sun et al., 2011; Bai et al., 2010; Zhang Z J et al., 2010; Gan et al., 2007; Li C et al., 2008; Royden et al., 2008; Lev et al., 2006; Zhang P Z et al., 2004; Wei W B et al., 2001).

The eastern part of the Tibetan Plateau belongs to Phanerozoic orogeny (Zheng et al., 2013; Xiao et al., 2010; Windley et al., 2007; Jahn et al., 2000; Şengör et al., 1993). Consequently, eastern Tibetan Plateau started lifting with time and reached its present-day elevation (more or less 3 000 m) by the Late Miocene or 11 Ma (Xu et al., 2016). Eastern Tibetan Plateau is topographically higher than the gently dipping south-eastern Tibetan Plateau thus their margin is mantled by the elevated low-relief relict landscape which was formed before the uplift of the eastern Tibetan Plateau (Clark et al., 2005). Previous observation and geodetic study shows that since the last 4 million years there has been a little crustal shortening along the eastern margin of the Tibetan Plateau (Royden et al., 1997). Also some results suggest that the uplift of the eastern Tibetan Plateau during the Late Cenozoic might have been produced by eastward over-thrusting and crustal shortening (Liu et al., 2017). The horizontal shortening and vertical stretching, in short 'weak lithosphere deformation' caused the plateau growth in the eastern Tibetan Plateau (Zhang F X et al., 2018). The topography, crustal thickness and sutures are assumed to be formed during the continental-continental collision of the Indian Plate and Asian Plate (Royden et al., 2008; Tapponnier et al., 2001; Yin and Harrison, 2000). Analysis from velocity variation and anisotropy reveals that the Tibetan Plateau might have extruded and its mantle materials might have flowed around the East Himalayan syntaxis (Pei et al., 2007). Another research shows lithosphere mantle thinning under the eastern Tibetan Plateau where the depth of the lithosphere-asthenosphere boundary is 140 km (Tunini et al., 2016). Recent ambient-noise tomographic result (Bao et al., 2020) suggests that, the eastern Tibetan plateau may be uplifted mainly by crustal shortening with later disordered crustal flow modulating the landscape. The topography of the eastern Tibetan Plateau is visibly different from northern and southern boundaries of the Tibetan Plateau thus eastern Tibetan Plateau has two long continuous slopes without distinct topographic boundaries (Burchfiel, 2004). Research reveals that, eastern Tibetan Plateau growth is induced by the push-squeeze mechanism based on a lithospheric scale (Zhang F X et al., 2018). The eastern margin of the Tibetan Plateau is surrounded by the Qaidam Basin and Qilian Orogen in the north, Alashan Block and Ordos Block in the northeast, Sichuan Basin and Yangtze Craton in the east. Few distinct sutures called the Jinsha River suture, Bangong-Nujiang suture and Indus-Tsangpo suture, Anyimaqen-Kunlun-Muztagh suture can be found across the eastern Tibetan Plateau (Fig. 1).

Figure  1.  Basic tectonic settings of the Tibetan Plateau with the station location indicated by red triangles. Sutures and thrust faults are created by using an open database available from previous research work (Styron et al., 2010; Taylor and Yin, 2009). The upper left map is an Azimuthal equidistant map showing 628 teleseismic events denoted by red circles, where the blue triangle indicates the study area. The epicentral distance between events ranges from 30° to 90°.

DownLoad: Full-Size Img PowerPoint

From the results of shear wave splitting measurements, it is evident that the northeastern part of the Tibetan Plateau is an undoubtedly complex region to portray all splitting measurements in a unique geodynamic model (Li et al., 2011). Also, seismic anisotropy results indicate that in Tibetan Plateau, mountain building processes caused significant deformation of the lithosphere beneath Tibetan Plateau but the asthenosphere seems to be less affected by the mountain building process (Huang et al., 2011). Many results show the diversity of this area. That is why a more detailed study can reveal unsorted and undefined facts.

2. DATA

Seismic data, as well as travel time data from the Namche Barwa Broadband Seismic Network (Meltzer, 2003) and ASCENT project, has been adopted in this study. Most of the data has been downloaded from IRIS (Incorporated Research Institutions for Seismology) and sorted according to the research area. Although datasets from these projects have been used in several studies (Yang et al., 2012; Zhang H et al., 2012, 2011; Fu et al., 2010; Singh and Kumar, 2009; Li et al., 2008; Sol et al., 2007) but a different approach with fast marching tomography method has been applied with the combined datasets, to find out the upper mantle heterogeneity. By using a MATLAB software package (Yu et al., 2017) with a waveform cross-correlation technique, at least 6 P-wave arrivals for each event were manually picked. As this research is based on teleseismic analysis, 628 distant earthquake data which has a magnitude of more than Mw 5.5 was considered (Fig. 1). These earthquakes have been recorded from 2003 to 2009 at 95 seismic stations. In total 8 532 P-wave arrivals were inverted for the upper mantle structure. The ray coverage in the study area is considerable (Fig. 2).

Figure  2.  A depth slice (at 400 km) and a vertical slice (93°E longitude) along with the ray coverage detected by stations.

DownLoad: Full-Size Img PowerPoint

3. METHODOLOGY AND INITIAL MODEL

For 3-D travel-time tomography, FMTOMO (fast marching tomography), a FORTRAN-90 software has been used (Rawlinson and Urvoy, 2006; Rawlinson and Sambridge, 2004a, b). FMTOMO software package is based on the fast marching method (Sethian and Popovici, 2012; Sethian, 1996). FMTOMO is based on a fast-marching method (FMM) for the forward step of travel-time prediction and also has a subspace inversion scheme to adjust model parameters for satisfying data observations. The eikonal equation of fast marching, which represents the propagation of seismic waves in the high frequency looks like this

Here, Δ_x is the gradient operator, T is travel time and s(x) is the slowness function of position x (Rawlinson et al., 2006). In short, the purpose of FMTOMO software is to let any velocity or interface structure be inverted by using most classes of travel time datasets (especially teleseismic which is used in this research work) (Rawlinson and Kennett, 2008). For the inversion step the iterative non-linear approach is the following equation (Kennett et al., 1988).

where m = model parameter, m₀ = starting model or reference model parameter g(m) = predicted residuals, d_obs = observed residuals, C_d = priori data covariance matrix, C_m = priori model covariance matrix, D = second derivative smoothing operator, ε = damping parameter, η = smoothing parameter.

These last two parameters ε, η controls the level of imposed regularization and also control the trade-off between the best way that satisfies data (Rawlinson and Kennett, 2008).

This multi-stage FMM (fast marching method) was given an extension form by introducing 3-D coordinates that included huge flexibility by allowing teleseismic phases to it (de Kool et al., 2006). This method is accomplished by computing travel time from distant source through a global reference model to the boundary nodes of a given or predefined model volume (Rawlinson and Kennett, 2008), in this case the global reference model is, ak135 (Kennett et al., 1995) which has been used in our research as well. More detailed information is available in these published papers (Rawlinson and Urvoy, 2006; Rawlinson and Sambridge, 2004a, b). For the reference model, the 1D reference or starting model ak135 (Kennett et al., 1995) has been used (Fig. 3).

Figure  3.  The starting or reference 1-D model ak135 (Kennett et al., 1995).

DownLoad: Full-Size Img PowerPoint

4. TOMOGRAPHIC RESULTS

4.1 Resolution Analysis

Synthetic resolution tests can evaluate the constraints of the dataset on the model parameters to justify the robustness of datasets. This means resolution tests give us a chance to find out where we have good data coverage and where we do not. Based on adequate ray-path coverage and arrival time uncertainty (Rawlinson and Kennett, 2008) we have conducted several resolution tests. The checkerboard test has been used in this research where the starting or initial model's checkerboard pattern is a combination of alternating fast and slow anomalies in the 3-D model (Rawlinson and Kennett, 2008; Rawlinson and Sambridge, 2003). Synthetic resolution tests showed good resolution in most places around the station location. Also, in some places where fewer stations are available but the good resolution is achieved due to the ray path coverage around those stations. The tomographic inversion scheme beneath the 95-station array is represented with a model volume comprising 46 464 velocity nodes with 33 km spacing in depth, 17.21 km in both latitude and longitude. Also, this model spans 700 km in depth, 12° in latitude and 17° in longitude. Before establishing tomographic results as well as a solution model, synthetic checkerboard tests (Rawlinson and Sambridge, 2003) have been done to find out the robustness of the model. The synthetic test was done using the identical sources, receivers and phase types the same as the observational dataset to determine the arrival-time residuals. These residual values are synthetic. After getting these values, they were inverted with the same procedure of solution model inversion (Fig. 4).

Figure  4.  Input model image and recovered model image in checkerboard resolution test.

DownLoad: Full-Size Img PowerPoint

4.2 Solution Model

Six iterations of the tomographic inversion scheme had been applied to the 8 532 relative arrival-time residuals to produce a solution model that satisfies the data to an acceptable level. The subspace dimension value was set at 20 for inversion where SVD (singular value decomposition) was used. The solution model has reduced data variance from 0.655 55 to 0.354 74 s², which represents a 46% variance reduction. This corresponds to an RMS reduction from an initial value of 809.61 ms to a final value of 595.57 ms. After running a big number of tests, the damping and smoothing parameters were selected by considering data variance. For the final solution model damping, ε = 5 and smoothing, η = 150 were imposed during inversion. Several depth slices have been created at 200, 300, 400 and 500 km (Fig. 5) from the tomographic solution of the eastern Tibetan Plateau.

Figure  5.  The tomographic solution model at 200 km (a), 300 km (b), 400 km (c) and 500 km (d) depth. In addition, because we use relative travel time residuals for the inversion, we usually adopt velocity variations in percentage (%) because we cannot derive absolute velocity variation.

DownLoad: Full-Size Img PowerPoint

5. RESULTS AND DISCUSSION

Figure 5 shows depth slices at 200, 300, 400 and 500 km depth respectively where high velocity anomalies can be seen on the southern and northern parts of the eastern Tibetan Plateau. In our results down to 200 km low velocity zone can be seen longitudinally from 92°E to 100°E and the corresponding latitude from 28°N to 38°N but at 30°N to 32°N, it seems to be disconnected by the high velocity zone. In between 30°N to 32°N and at 94°E longitude low velocity anomalies were observed significantly. At 300 km depth, the distribution of low velocity anomaly seems more dispersed. But an intense low velocity zone is observed under Songpan-Ganzi terrain in all depth slices. At 200, 300, and 400 km depth under the Bangong-Nujiang suture, low velocity anomalies can be seen to be disconnected at 30°N to 32°N by high velocity anomalies. This high velocity anomaly is supposed to be the Indian lithosphere (Replumaz et al., 2013; Zhao et al., 2011; Li et al., 2008; Tilmann et al., 2003). Also a very high velocity anomaly can be seen at 200, 300, and 400 km depth slices spread from 37°N to 40°N, which is beneath Qilian Orogen indicating possible Asian lithosphere (Zhao et al., 2011; Willett and Beaumont, 1994).

There have been a good number of studies conducted in the eastern Tibetan Plateau. For example, the upper mantle S velocity model created from Rayleigh wave tomography also shows two subducting continental plates facing opposite to each other indicating that the Indian Plate is subducting under the southern Tibetan Plateau and the Asian Plate is subducting under the northern Tibetan Plateau (Li Y H et al., 2013). P-wave speed is observed low beneath the central and eastern Tibetan Plateau (Li et al., 2008). GPS measurements also suggest that mantle flow promotes block extrusion under Tibetan Plateau although hot asthenospheric upwelling under Tibetan Plateau can lead to post-collisional magmatism under the southern Tibetan Plateau widely (Zou et al., 2019). Another study suggests that low velocity anomalies of the upper mantle under the Songpan-Ganzi terrain portray hot and wet mantle upwelling which possibly occurred and was induced via the Moho interface to add fluids to the crustal flow of the eastern Tibetan Plateau (Lei et al., 2019). Also it is suggested that the upwelling counter flow might be the reason for the low velocity zone beneath the Songpan-Ganzi terrain (Yue et al., 2012). The Indian lithospheric mantle is underthrusting sub-horizontally beneath the eastern Tibetan Plateau also on the northern side Asian lithospheric mantle is detected close to Qaidam Basin (Zhang et al., 2012; Nábělek et al., 2009; Tilmann et al., 2003; Kind et al., 2002). This Asian lithospheric mantle is also clearly visible in our results. From a recent study by (Zheng et al., 2019) it is assumed that seismic velocities of the uppermost mantle are slow in the northeastern Songpan-Ganzi terrain in comparison to the south although the outward expansion of the Tibetan Plateau might have evolved by a combination of crustal flow along with the southeastern margin. Zheng et al. (2019) also suggested that mantle upwelling could be triggered by the removal of the thickened lithosphere in its eastern marginal places.

Velocity heterogeneity of the eastern Tibetan Plateau can be portrayed by observing our tomographic results. Beneath the Bangong-Nujiang suture, where the Indian Plate is subducting, a sudden change in velocity can be observed (Fig. 6). Another sudden change in velocity can be seen under Qaidam Basin and Qilian Orogen, where possible Asian lithosphere is present. Also, between 33°N to 35°N low velocity zone seems to be disconnected by a relatively high velocity zone. This portion could be a part of subducting lithosphere or some detached block of pre-existed lithosphere probably causing large scale thermal convection, which might be a reason for plateau growth and uplift. A recent study shows strong evidence of the convective removal of the lithosphere (what we called detached patches in our results) beneath the Bangong-Nujiang suture in the central Tibetan Plateau (Xiao et al., 2020). In between the subducting Indian Plate and Asian Plate, low velocity anomalies can be seen from 200 to 500 km mostly (in some places more than 500 km) (Fig. 5 and Fig. 6). From results of other parts of the eastern Tibetan Plateau it is assumed that the low velocity zone ranges from 200 to 400 km on average (Fig. 6). Although in some places this depth increases down to 660 km (Fig. 6). This low velocity zone ranging from 200 to 400 km is interpreted as a possible mantle wedge in a prominent recent study (Xiao et al., 2020; Lei et al., 2019). In our results we interpret this zone as hot asthenosphere disconnected by subducted continental lithosphere.

Figure  6.  Vertical slice along 97°E or AA' (profile location has been drawn in Fig. 1).

DownLoad: Full-Size Img PowerPoint

The main objective of this research work contains the determination of velocity variation of the upper mantle along the eastern part of the Tibetan Plateau. From these observations, the abrupt change in the low velocity zone can be seen in the middle of our study area (beneath the Bangong-Nujiang suture and Songpan-Ganzi terrain in the eastern part of the Tibetan Plateau) ranging from 200 to ~500 km depth. Following previous studies, the large-scale high velocity anomaly in the southern part and northern part of our image is explained as the subducting Indian Plate and Asian Plate accordingly. The subducted Indian Plate is visible with low velocity anomalies at longitude 94°E may be caused by slab tearing because it is located at EHS (eastern Himalayan syntaxis) (Peng et al., 2016; Li Z H et al., 2013).

Another feature is the mantle upwelling in and around the eastern Tibetan Plateau, which is sandwiched by two continental lithospheres (Indian lithosphere and Asian lithosphere). In between these two continental lithospheres, some high velocity anomalies of detached patches (remains of a plate that is supposed to be pre-existed or possibly detached from subducting plates) are visible which might be causing thermal convection and plateau uplift. Down to 500 km low velocity zones can also be observed with some lithospheric fragmentation. Hence the low velocity anomalies, which are detected in our results, may reflect the relatively hot asthenosphere. Moreover, either the asthenosphere upwelling or the mantle wedge can occur here due to the presence of the cold lithospheric mantle in the surrounding areas as we mentioned above.

6. CONCEPTUAL MODEL

From 92°E to 100°E good resolution has been achieved. To visualize the tomographic results in a 3-D way, a conceptual model has been established to portray our results distinctly. In consideration of depth slices and vertical slices (both along latitude and longitude) a 3-D model has been built (Fig. 7). In our model dark blue or navy-blue colored features are considered as two prominent continental lithospheres such as the Indian and Asian plates. The blue colored feature is considered as patches of the detachment of a plate that is supposed to be pre-existed or possibly detached from subducting plates, and the red colored features are low velocity anomalies. Possible mantle upwelling is indicated by black arrows. For the most part, a low velocity zone means weak material there, so it is easy to uplift under the strong pushing by the Indian Plate and form a plateau with a large range. Our results suggest the low velocity zone in the central Tibetan Plateau plays an important role to form a large and high plateau.

Figure  7.  Conceptual model created by considering tomographic results.

DownLoad: Full-Size Img PowerPoint

7. CONCLUSIONS

Teleseismic events recorded on various seismometers deployed by different projects and organizations have been used to visualize the 3-D P-wave velocity structure of the upper mantle of the eastern Tibetan Plateau. To reveal the image of the upper mantle structure of the eastern Tibetan Plateau, travel time residuals have been used for tomography by following the eikonal solution of the fast-marching method. Teleseismic data consists of 628 events from the period of 2003–2009 were focused mostly. From the tomographic results some points are summarized.

(1) A high velocity anomaly in the southern part indicates subducting Indian Plate or lithosphere, which shows low velocity at 94°E longitude, which might be caused by slab tearing as it is located on the eastern Himalayan syntaxis (EHS).

(2) High velocity anomaly in the northern part represents Asian Plate or lithosphere. This anomaly is observed down to 500 km depth at 37°N to 40°N beneath Qaidam Basin and Qilian Orogen. In between the subducting Indian Plate and Asian Plate, low velocity anomalies can be seen from a depth of 200 to 500 km mostly (in some places more than 500 km).

(3) At the central part of the study area high velocity detached patches surrounded by low velocity anomalies trigger the possibility of thermal convection and plateau uplift. This scenario is observed down to 500 km. The abrupt change in velocity can be seen between these patches and the Bangong-Nujiang suture. This indicates a possible asthenospheric upwelling. Also, another abrupt change in velocity is observed between detached patches and Songpan-Ganzi terrain on the eastern part of the Tibetan Plateau. This might represent the existence of a mantle wedge. In both cases these low velocity anomalies range from a depth of 200 to ~500 km.

We conclude that the low velocity zone in the central Tibetan Plateau play an important role to form large and high plateau and also the low velocity anomalies detected in our results may reflect either the hot asthenosphere upwelling or the mantle wedge due to the presence of the cold lithosphere.