Lithospheric Structures beneath the Tibetan Plateau Constrained by Sn Propagation Efficiency Open Access

After the Neo-Tethyan Ocean closed, the Indian continent collided with the Asian continent (50-60 Ma) [1], resulting in shortening deformation of the continent [2–6] and uplift of the Tibetan Plateau on a scale of 1500-3000 km, as shown in Figure 1 [7–9]. There is still no consensus on the tectonic mechanisms for plateau uplift

and deformation. Previous studies [10–12] suggested that the Indian continental lithosphere subducted under the Tibetan Plateau from the south, while the Asian continental lithosphere subducted beneath the Tibetan Plateau from the north, causing crustal thickening. It was proposed that shortening leads to thickening of the crust and uppermost mantle of the Tibetan Plateau [13–16], and a large amount of shortening and thickening may make the lithospheric mantle unstable, resulting in partial or total delamination. Later, the uplift of the Tibetan Plateau was attributed to the upwelling of hot material from the base of the upper mantle [17, 18], but the uplift of the surrounding terranes was caused by multidirectional continental subduction. Previous results from tomographic images [19–22] reveal that the uppermost mantle in the northern and southeastern Tibetan Plateau shows obvious low-velocity characteristics, implying that there is a high-temperature structure with strong seismic attenuation in the uppermost mantle. Based on this, some researchers consider that asthenospheric hot material beneath the Tibetan Plateau might upwell and consequently partially melt in the uppermost mantle. However, the spatial distribution of seismic attenuation structures, resulting from partial melts in the uppermost mantle, is still unclear and requires additional seismic evidence.

Pn or Sn waves propagate through the uppermost mantle [23], and their propagation is closely related to the seismic attenuation structure, composition, and velocity gradient of the uppermost mantle [11, 24]. Seismic attenuation is usually compelling evidence for high temperature or partial melting structures or their combination [25, 26]. Since a weak seismic waveform amplitude is a direct reflection of strong seismic attenuation, many previous works have discussed the characteristics of seismic attenuation through Sn wave amplitude information in the Tibetan Plateau [27–29]. The Sn wave propagation efficiency was qualitatively evaluated by comparing the Sn amplitude to the Pn amplitude [27]. They first reported an inefficient Sn wave zone in the Qiangtang terrane, indicating that asthenospheric material might intrude upward to a position near the base of the crust.

Figure 2 shows an example of Sn waves propagating in the Tibetan Plateau. The source is located in the southwestern Tibetan Plateau; at station XC.SANG, waves do not pass through the Qiangtang terrane, and at XC.WNDO, waves pass through the Qiangtang terrane. Please see Figure 1 for the locations of events (red stars and black dots) and stations (blue and cyan triangles). The Sn signal of station XC.SANG does not weaken, but that of station XC. The WNDO weakens significantly with increasing frequency, indicating the obvious frequency dependence of Sn wave propagation in the Tibetan Plateau.

The inefficient Sn wave zone is consistent with the location of Quaternary basaltic volcanism in the Qiangtang terrane, which suggests high temperature and attenuation structures for this area [27]. With the ratio of Sn wave amplitude to Pg coda wave amplitude from seismic data for 11 stations and 173 events on the Tibetan Plateau, McNamara et al. [28] quantitatively measured Sn wave propagation efficiency. They found a wide range of inefficient Sn wave zones and pointed out the frequency dependence of Sn wave propagation in the Tibetan Plateau. With developments in seismic observation systems, more seismic waveform data in the Tibetan Plateau were used to obtain the distribution of Sn wave propagation efficiency [29]. They improved the spatial resolution of Sn wave propagation efficiency and studied the frequency dependence of Sn wave propagation efficiency in more detail.

There remain some shortcomings in previous research work. (1) Ray density is low in the northeastern Tibetan Plateau, leading to unclear northern and eastern boundaries of the inefficient Sn wave zone; (2) only the lateral variation in Sn wave propagation efficiency in the Tibetan Plateau is obtained, and there are few constraints on the depth variations in Sn wave propagation efficiency. Here, we utilize the characteristics of Sn propagation to obtain more specific spatial variations in the seismic attenuation structure. To show the outline of the inefficient Sn wave zone more clearly, we obtain the Sn wave propagation efficiency from a total of 31,538 earthquake-station pairs across the entire Tibetan Plateau after data selection. We apply the fast marching method [30–33] to track ray paths in the constructed 3D model, in which the uppermost mantle structures have been extensively investigated (e.g., [34–36]). This investigation would help constrain the depth variations in Sn wave propagation efficiency.

Seismic data for 3272 events and 247 stations on and around the Tibetan Plateau (Figure 1) in this study are from the Incorporated Research Institutions for Seismology (IRIS). To ensure the quality and reliability of seismic waveforms, earthquake magnitudes are greater than 4 (Richter scale, earthquake catalog from the U.S. Geological Survey (USGS)). This study uses events located in the crust to pick Sn phases. According to the CRUST 1.0 model [37], focal depths shallower than 65 km and 35 km inside and outside the Tibetan Plateau are selected so that all paths would pass through the uppermost mantle. We finally collect waveforms for 31,538 event-station pairs during 1991-2010.

Sn wave propagation efficiency is calculated from the amplitude ratio of the Sn wave to Pg coda in specific windows as $r=Sn/Pgcoda$ from the transverse component of seismic data. As the averaged Pg-wave velocity is ~6.4 km/s beneath Tibet [28, 29], the Pg coda window is estimated by the velocity ranging from 5.0 to 6.0 km/s. We define the Sn window with a velocity between 4.2 and 5.0 km/s [28, 29]. As illustrated in Figure 3, the time window of the Sn wave is marked by a red segment (⁠$TS$⁠), the time window of the Pg coda wave is marked by a blue segment (⁠$TP$⁠), and $T0$ is the intercept time of the Sn phase.

As the amplitude of Pg coda is generally steady, we take the root-mean-squared (RMS) amplitude of the envelope of Pg coda. The Sn window is split into five equal time segments, and the Sn amplitude is taken as the maximum RMS amplitude of the envelope function in each segment. In addition, we calculate the signal-to-noise ratio, i.e., (Pg coda)/noise, to remove noisy waveforms. The noise is defined in the time window of [-35, -20] s relative to the predicted Pn arrival from the ak135 model. Those waveforms with SNRs smaller than 2.0 are discarded. After selection, the raypath number in each frequency band is listed in Table S1.

Pn/Sn and Pg/Sg must be strictly separated to calculate the Sn wave propagation efficiency. To constrain the minimum epicentral distance for the seismic data, we apply the fast marching method (FMM) [30, 31] to simulate the travel time curves of Pn/Pg and Sn/Sg (Figure 4). The models of the modified ak135 (Figures 4(a) and 4(c)) [38] and ESUPP (Figures 4(b) and 4(d)) [39] are both involved in simulations to constrain the crossover distance.

As shown in Figure 4, the epicentral distances where Pn precedes Pg according to the two models are both approximately 3.65° (Figures 4(a) and 4(b)), and those where Sn exceeds Sg are approximately 3.70° for model ak135 (Figure 4(c)) and 3.90° for model ESUPP (Figure 4(d)). In previous work on Pn tomography, the separation of Pn and Pg was considered to be 3.0° [40]. To ensure that the Sn wave is the first S arrival and that the Pg wave is not affected by the Pn wave, we therefore select 4° as the minimum epicentral distance for seismic data to discern Sn phases as the first S arrival.

In Figure 5, we first show a scatter diagram of the distribution of Sn wave propagation efficiency values at 1.0-5.0 Hz in the Tibetan Plateau. Figure 5(a) shows that the efficiency values of small epicentral distances (4°–7°) are mainly less than 3, those of middle epicentral distances (7°– 12°) are mainly less than 2, and those of large epicentral distances (12°–16°) are mainly less than 1, indicating that most Sn wave propagation efficiency values are no greater than 3.

To further investigate the distribution of Sn wave propagation efficiency values, we divide the study region into 16 local areas (Figure 5(b)) and then count the number of raypaths for each Sn wave propagation efficiency value in every local area (Figure 6). As shown in Figure 6, the number of raypaths with Sn wave propagation efficiency values $r<2$ in regions 7, 8, 11, 12, and 16 accounts for approximately 93.75%, 90.82%, 95.64%, 91.51%, and 89.28%, respectively. The number of raypaths with $r<2$ in other regions decreases significantly, especially in regions 1, 5, and 10, where the proportions of raypaths with $r<2$ are 59.61%, 41.80%, and 59.70%, respectively.

McNamara et al. [28] and Barron and Priestley [29] took $r=2$ as the boundary between high and low efficiency. As the proportions of raypaths change obviously (Figure 6), we also consider the efficiency value $r=2$ to be the critical value for inefficient/efficient propagation of the Sn wave. To validate our work, we perform the same workflow with the same parameters and data as Barron and Priestley [29]. We obtain a pattern of low Sn propagation efficiency similar to that of Barron and Priestley [29], which is provided in Supplementary Text S1.

As illustrated in the previous section, we consider $r=2$ to be the critical value for inefficient/efficient propagation of the Sn wave. After filtering seismic data with a five-pole Butterworth bandpass filter, we extract the amplitude information and obtain all the Sn wave propagation efficiency values at 0.3-0.4 Hz, 0.6-0.7 Hz, 0.8-1.0 Hz and 1.5-2.0 Hz. We then calculate the SNR and select raypaths with $SNR>2$⁠, which leads to 14,662, 22,103, 24,168, and 20,087 raypaths in each frequency band (Table S1). We finally derive the lateral variations in the mean Sn wave propagation efficiency at these frequency bands in the Tibetan Plateau (Figures 7(a), 7(d), 7(g), and 7(j)). These maps are gridded at 0.25°×0.25°, and the mean value of Sn wave propagation efficiency is then calculated in each cell. In the end, a simple smoothing of 0.5°×0.5° is applied, which is the same as Figure S1b.

As shown in Figure 7, at low frequencies (0.3-0.4 Hz), a region of weak inefficient Sn wave propagation appears in the northern-central parts of the Qiangtang terrane and the western Songpan-Ganze terrane. With increasing frequency, the inefficient Sn wave propagation in the northern-central parts of the Qiangtang terrane and the western Songpan-Ganze terrane gradually becomes stronger. The strong inefficient Sn wave propagation extends to the east of the Songpan-Ganze terrane (0.6-0.7 Hz) and extends to the Sanjiang area and the southern margin of the Qilian orogenic belt (0.8-1.0 Hz, 1.5-2.0 Hz).

Compared with previous research results [28, 29], this study outlines the inefficient Sn wave zone more clearly, especially at the northern and eastern boundaries. Since we use more raypaths to the south of the Qilian orogenic belt, a small inefficient Sn wave zone on the southern margin of the Qilian orogenic belt (37° N/102° E) is newly identified (Figures 7(g) and 7(i)). Due to the low raypath density in the southwestern and northeastern corners of the Tibetan Plateau (Figures 7(b), 7(e), 7(h), and 7(k)), the lateral variations in the mean Sn wave propagation efficiency in these two regions are not discussed in this study.

The crust is as thick as 65-70 km beneath the central Tibetan Plateau [37]. Earthquakes in the upper crust and the lower crust might have different influences on the Sn wave propagation efficiency. We examine the focal depth and find that most earthquakes (2,813 out of 3,272) are shallower than 35 km. We use those earthquakes shallower than 35 km to map the mean Sn wave propagation efficiency (Figure S2) with the same mapping process as Figure 7. Furthermore, we also map the mean Sn wave propagation efficiency (Figure S3) from two groups of events (focal depths in the ranges of 0-10 km and 30-60 km) with comparable earthquake numbers. All the results show a similar pattern of Sn wave propagation efficiency. It is concluded that the focal depth has little influence on the calculation of the Sn wave propagation efficiency.

Previously, the Sn propagation efficiency was the mean value of all the propagation efficiency values in each cell. We further examine the robustness of our results by using the median instead of the mean Sn propagation efficiency with the same mapping process as above, as shown in Figure S4. Comparing Figure S4 with Figure 7 (left), while there is little difference on the south boundary of the inefficient Sn wave zone at higher frequencies (0.8-1.0 Hz, 1.5-2.0 Hz), the distributions of Sn propagation efficiency are similar. It is concluded that the use of the mean and median values of the Sn wave propagation efficiency does not change our interpretation.

Sn phases propagate deeper into the uppermost mantle as the epicentral distance of events increases [34–36]. As shown in Figure S5, the 3D Sn phases can propagate downward to 130 km at an epicentral distance of 16° with the flat Moho at 65 km, i.e., the averaged Moho of the CRUST 1.0 model [37]. These results illustrate that the propagation of Sn waves is much more complicated than the simplified assumption of running along paths immediately below a flat Moho.

Based on the ESUPP velocity model [36], we utilize the FMM [30, 31] to carry out a 3D forward simulation of the raypaths and travel time. The input velocity model is interpolated into a uniform lateral interval of 0.5°◊0.5° and depth interval of 5 km from the original model. The flat Moho does not change the pattern of Sn propagation efficiency, as raypaths are primarily dependent on the 3D velocity model. We finally obtain 31,538 valid raypaths (Table S1). Figure 8 shows three vertical profiles of 3D raypaths projected into a west-east direction within a distance range of 2° perpendicular to the profile. According to this 3D simulation, we demonstrate the mean Sn wave propagation efficiency variations at different depths below the Moho interface in Figure 9.

In the northern-central parts of the Qiangtang terrane and the western Songpan-Ganze terrane (approximately 32°N-36°N and 84°E-92°E), the strong inefficient Sn wave structure can extend down to 130 km (lower lithosphere) from a position near the base of the crust (Figures 8(c) and 9(a); Figure S9a). In the eastern Songpan-Ganze terrane and the Sanjiang area (approximately 29°N-34°N and 95°E-103°E), the raypaths with inefficient Sn waves are mainly concentrated from 90 to 130 km (Figure 8(b)). Accordingly, the region of inefficient Sn waves here is mainly concentrated at 90-130 km (Figures 9(b) and 9(c); Figure S6b, c) and becomes weaker at shallower depths (Figure 9(a); Figure S6a). Since there are few raypaths involved in the 3D travel time simulation on the southern margin of the Qilian orogenic belt (37° N/102° E), the depth range of Sn wave propagation efficiency is not discussed.

Efficient Sn wave propagation corresponds to relatively cold seismic structures with low attenuation in the uppermost mantle. In contrast, inefficient Sn wave propagation corresponds to relatively high temperatures or partial melting in the uppermost mantle [24, 28, 29].

In this study, we obtain a much clearer outline of the inefficient Sn wave zone in the Tibetan Plateau and reveal larger distributions of the seismic attenuation structure in the uppermost mantle of the Tibetan Plateau. As illustrated in Figure 7, the northern-central parts of the Qiangtang terrane, the Songpan-Ganze terrane, the Sanjiang area, and the southern margin of the Qilian orogenic belt present distinctly inefficient Sn features, which provide excellent evidence for a strong seismic attenuation structure or partial melting in the uppermost mantle.

Previous seismic imaging results, such as Pn [40] or S wave velocity models [19–21, 41], show that the uppermost mantle beneath the northern and southeastern Tibetan Plateau exhibits obvious low-velocity characteristics, indicating that there may be high-temperature structures with strong seismic attenuation. These features in tomographic images are consistent with our results from Sn propagation efficiency. Moreover, the inefficient Sn wave zone from our results in the Tibetan Plateau matches well with the geological boundaries around the Tibetan Plateau, including the southern boundary of the Alashan block, the northern edge of the Bangong-Nujiang suture, and the western boundary of the Ordos block and Yangtze block, which are also similar to the tomographic images [19–21, 40, 41]. This indicates that there may be strong deformation coupling between the lithospheric mantle and the crust in the Tibetan Plateau [26, 40, 41].

The depth of the seismic attenuation structure or partial melt in the uppermost mantle has varied effects on Sn waves at different frequencies [28, 29]. According to the frequency dependence of Sn wave propagation [29] in the Tibetan Plateau (Figure 7), the strong seismic attenuation zones in the uppermost mantle of the northern-central parts of the Qiangtang terrane and the western Songpan-Ganze terrane are relatively thicker, while those in the Sanjiang area, the eastern Songpan-Ganze terrane and the southern margin of the Qilian orogenic belt are thinner.

Furthermore, we specifically constrain the depth range of the strong seismic attenuation in the uppermost mantle (Figures 8 and 9). The strong seismic attenuation structures in the uppermost mantle of the northern-central parts of the Qiangtang terrane and western Songpan-Ganze terrane can extend down to 130 km from almost the base of the crust, indicating increased temperature or partial melting. Similarly, those of the eastern Songpan-Ganze terrane and the Sanjiang area are mainly concentrated at 90-130 km, indicating an increased temperature or a lower degree of partial melting [27, 28, 40].

Due to the low ray coverage involving 3D ray forward simulation on the southern margin of the Qilian orogenic belt (Figure 8, A3, 102°E), this study poorly constrains the Sn wave propagation depth there. However, according to the lateral variations in Sn wave propagation efficiency (Figure 7), the inefficient Sn wave zone at high frequencies in this region, with sufficiently high raypath density, can prove that partial melting of the uppermost mantle might occur in a thinner layer.

According to the lateral and vertical variations in Sn wave propagation efficiency, this study reveals the spatial variations in seismic attenuation structures of the lithospheric mantle in the Tibetan Plateau. Seismic attenuation structures appear from the Moho downward to a depth of 130 km beneath the northern Qiangtang terrane and western Songpan-Ganze terrane. The attenuation structure tends to a thinner layer (90-130 km) along the southern boundary of the Alashan block and the northern edge of the Bangong-Nujiang suture. After encountering the western boundary of the Yangtze block, the seismic attenuation structures turn to the southeast of the Tibetan Plateau.

From the northwestern to southeastern Tibetan Plateau along the Bangong-Nujiang suture, the potassic basaltic magma assemblage (45-30 Ma) presents a cluster-type intermittent distribution (areas A, B, and C in Figure 10(a)) [42–46]. These basaltic magmas were derived from the fertile lithospheric mantle metasomatized and the upwelling of the asthenosphere. Moreover, potassium alkaline-carbonate rock assemblages (27-7 Ma) are present in the Ailaoshan area and southeast of the Qilian Orogen (areas D and E, Figure 10(a)) [42–46]. This type of magma assemblage results from the carbonation of metasomatized lithospheric mantle at the edge of a rigid block, reflecting the lateral flow of upwelling asthenospheric material [1, 46]. This indicates that the seismic attenuation structure of the lithospheric mantle in the Tibetan Plateau is quite consistent with the series of magmatic assemblages. The thicker seismic attenuation structure (from the Moho to a depth of 130 km) of the lithospheric mantle in the northern-central parts of the Qiangtang terrane and the western Songpan-Ganze terrane corresponds to the older potassic magmatic association (45-35 Ma) [1, 45, 46]. The seismic attenuation structure in a thinner layer (90-130 km) of the lithospheric mantle in the Sanjiang area and the western Songpan-Ganze terrane corresponds to the relatively new potassic magmatic assemblage (40-32 Ma) and the potassium alkaline carbonate assemblage (27-7 Ma) in the Ailaoshan area and the southeastern Qilian orogen [1, 45, 46].

Previous studies based on seismic tomography and receiver functions [47, 48] showed that the subduction front of the Indian Plate has reached the Bangong-Nujiang suture, which agrees well with the south boundary of the inefficient Sn wave zone in this study (Figure 7). The uppermost mantle appears as a high-velocity feature beneath the Lasha terrane in the south, while it shows a low-velocity feature beneath the Qiangtang terrane in the north [19–21, 41]. The observations are consistent with the distribution of the inefficient Sn propagation zone in our results (Figure 7). These results indicate that the upwelling of asthenospheric material beneath the Qiangtang terrane may have occurred due to the subduction of the Indian Plate [48], leading to a thinner lithosphere beneath the central Qiangtang terrane [49].

Combining previous information with our results, we suggest that upwelling and lateral flow of asthenospheric material have occurred beneath the Tibetan Plateau (Figures 10(b) and 10(c)), providing a thermal driving force for the uplift of the Tibetan Plateau and the deformation of the crust-mantle structure [1, 46]. After the northward subduction of the lithospheric mantle of the Indian continent (Figure 10(b)) [47, 50, 51], the asthenospheric material inevitably upwells beneath the central Qiangtang terrane according to the principle of material conservation [48]. Due to the obstruction by surrounding rigid rock layers (Figure 10(b)), such as the Tarim block, Alashan block, Ordos block, and Yangtze block, this upwelling asthenospheric material flows laterally to the east of the Songpan-Ganze terrane, the Sanjiang area, and the southern margin of the Qilian orogenic belt (Figure 10(c)) [1, 46], consistent with the SKS splitting measurements [47, 52]. The upwelling and lateral flow of asthenospheric material beneath the Tibetan Plateau caused partial melting of the lithospheric mantle (Figure 10(c)) [46]. Given these processes, the lithospheric mantle injects magma derived from the mantle into the plateau crust through partial melting, which contributes to the growth of the plateau crust [1, 44]. Moreover, the thermal driving force from asthenospheric material upwelling and lateral flow and the extrusion of surrounding rigid continental blocks may provide an explanation for the eastward and southward extrusion of material in the Tibetan Plateau [46].

Based on 31,538 Sn arrivals from 3,272 events and 247 stations, we obtain the lateral variations in Sn wave propagation efficiency across the entire Tibetan Plateau and surroundings, including the northern-central parts of the Qiangtang terrane, the Songpan-Ganze terrane, the Sangjia area, and the southern margin of the Qilian orogenic belt, revealing the widespread distribution of the seismic attenuation structure in the uppermost mantle. Moreover, we obtain the depth variations in Sn wave propagation efficiency by simulating the 3D raypath of 31,538 arrivals, which restricts the specific depth range of the seismic attenuation structure in the uppermost mantle of the Tibetan Plateau; the attenuation zone in the uppermost mantle in the northern-central parts of the Qiangtang terrane is thicker, extending from the base of the crust to 130 km, while that in the eastern Songpan-Ganze terrane and Sanjiang area is thinner, mainly concentrated at 90-130 km, and does not extend upward to the base of the crust. Combining our results with previous studies, we suggest that upwelling and lateral flow of asthenospheric material might occur in the northern-central parts of the Qiangtang terrane and the western Songpan-Ganze terrane. Due to obstruction by the surrounding rigid rock layers, the upwelling asthenospheric materials then flow laterally to the eastern Songpan-Ganze terrane, the Sanjiang area, and the southern margin of the Qilian orogenic belt. The upwelling and lateral flow of asthenospheric material provide a thermal driving force for the uplift of the Tibetan Plateau and the deformation of the crust-mantle structure, maintaining crustal growth and driving the extrusion of plateau material toward the southeast.

The data is available at IRIS Data Management Center (http://ds.iris.edu/ds/nodes/dmc) for retrieving waveform data recorded by permanent and portable stations.

The authors declare that they have no conflicts of interest.

The waveforms are retrieved from IRIS Data Management Center (http://ds.iris.edu/ds/nodes/dmc), including stations from networks XC, XD, XF, XG, XR, XW, X4, YA, YL, and ZV. The work is supported by the National Natural Science Foundation of China (grant nos. 41774060, 42022026, and 41720104006). The support from the Youth Innovation Promotion Association CAS is also acknowledged. The Supercomputing Laboratory of Institute of Geology and Geophysics, Chinese Academy of Sciences, is thanked for the computational resources.