Complex structure of upper mantle beneath the Yadong-Gulu rift in Tibet revealed by S-to-P converted waves

JID:EPSL

AID:115954 /SCO

[m5G; v1.279; Prn:15/11/2019; 12:04] P.1 (1-9)

Earth and Planetary Science Letters ••• (••••) ••••••

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

Complex structure of upper mantle beneath the Yadong-Gulu rift in Tibet revealed by S-to-P converted waves Zhen Liu a,b , Xiaobo Tian a,c,∗ , Xiaohui Yuan d , Xiaofeng Liang a,b , Yun Chen a,b , Gaohua Zhu e , Hongshuang Zhang f , Wei Li a,b , Ping Tan g , Sicheng Zuo a,b,h , Chenglong Wu a,b , Shitan Nie a,b,h , Gaochun Wang a,b,h , Guiping Yu a,b,h , Beibei Zhou a,b,h a

State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China Innovation Academy for Earth Science, CAS, Beijing 100029, China c CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing 100101, China d Deutsches GeoForschungsZentrum-GFZ, Telegrafenberg, 14473 Potsdam, Germany e Earth System Science Programme, Faculty of Science, Chinese University of Hong Kong, Hong Kong, China f Institute of Geology, Chinese Academy of Geological Sciences, Beijing 100037, China g Institute of Disaster Prevention, Hebei, Sanhe 065201, China h University of Chinese Academy of Sciences, Beijing 100049, China b

a r t i c l e

i n f o

Article history: Received 1 June 2019 Received in revised form 24 September 2019 Accepted 5 November 2019 Available online xxxx Editor: A. Yin Keywords: Tibetan plateau Yadong-Gulu rift structure of upper mantle S-to-P converted waves

a b s t r a c t The convergence between the Indian and Eurasian plates has produced the thick crust and uplifted the Tibetan plateau since about 50 Ma. However, the deformation of the mantle lithosphere is still a puzzle. The geometry of the subducting Indian mantle lithosphere beneath the plateau and the thickening or/and delamination of the Tibetan mantle lithosphere are the keys for understanding the continental collision process and the evolution of the plateau. However, knowledge has been restricted due to sparse data coverage in Tibet. In this study, S-wave receiver functions are calculated using tele-seismic waveforms recorded by two broadband arrays in central Tibet to image the lithospheric structure, mainly the depth variation of the lithosphere-asthenosphere boundary (LAB). Our results show that the depth of the Tibetan LAB decreases from ∼150 km in the west to ∼120 km in the east across the north-south trending Yadong-Gulu rift. Similarly, the LAB depth of the subducting Indian slab decreases from ∼270 km in the west to ∼200 km in the east, and the northernmost subducting Indian slab can be observed beneath the Bangong-Nujiang suture. These observations suggest that the thickness of the Tibetan lithosphere and the depth of the underlying Indian slab are segmented across the Yadong-Gulu rift in different degrees. The abrupt changes imply that the subducting Indian slab has been torn, which provided an upwelling channel for the asthenosphere contributing to the development of the rift. © 2019 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction The Tibetan Plateau has been formed by the collision between the Indian and Asian plates and their subsequent convergence since approximately 50 Ma (e.g., Molnar and Tapponnier, 1978; Mulch and Chamberlain, 2006; Tapponnier et al., 2001). The geometrical morphology of the subducting Indian mantle lithosphere (SIML) is one of the most important pieces of information to un-

*

Corresponding author at: State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China. E-mail address: [email protected] (X. Tian).

derstand the collision and the consequent mountain-building process. This geometry has been studied extensively as a major objective of geophysical observations, such as broadband seismic arrays. The SIML is usually imaged in two ways: imaging the lithosphere itself by tomography or detecting the lithosphere-asthenosphere boundary (LAB) by receiver functions. Several north-trending receiver function profiles across the Tibetan Plateau have imaged the LAB (e.g., Kumar et al., 2006; Zhao et al., 2010, 2011), and many seismic tomography results have been obtained for the upper mantle structures beneath the Tibetan Plateau (e.g., Agius and Lebedev, 2013; Chen et al., 2017; Li et al., 2008; Li and Song, 2018; Liang and Song, 2006; Liang et al., 2016a). However, these results show significant differences in the geometry of the SIML due to various station coverages, diverse seismic phases used and varying sensi-

https://doi.org/10.1016/j.epsl.2019.115954 0012-821X/© 2019 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

JID:EPSL

2

AID:115954 /SCO

[m5G; v1.279; Prn:15/11/2019; 12:04] P.2 (1-9)

Z. Liu et al. / Earth and Planetary Science Letters ••• (••••) ••••••

Fig. 1. Tectonic map around the Tibetan Plateau. Geological structures are based on Styron et al. (2010). Pink filled circles mark distribution of magmatic episodes at ∼150 Ma (Chung et al., 2005). The brown filled circles mark high 3 He/4 He ratio anomaly sample points (Hoke et al., 2000; Hou et al., 2004). The yellow area indicates the low seismic velocity zone in the upper mantle beneath the Yadong-Gulu rift (YGR) (Liang et al., 2016a). The violet dashed line delineates a 2.5% anomaly contour at a depth of 200 km beneath Tibet from surface wave tomography (Agius and Lebedev, 2013). The green dashed rectangles delineate the Moho offset (Zhang and Klemperer, 2005; Tian et al., 2015). The short green bar marks the highly electrically conductive body in the lower crust (Wang et al., 2017). The violet lines delineate the LAB of the subducting Indian lithosphere (Zhao et al., 2010, 2011), and the dark green line denotes the locations of their observation profiles. Black arrows indicate motion of the Indian Plate relative to stable Eurasia. Blue triangles and red inverted triangles denote stations of the SANDWICH and Tibet-31N seismic arrays, respectively. The locations of these stations are listed in Table S1. The bottom-right inset is a map showing the location of the events used in this study. The red rectangle is the study area. The blue pluses and red crosses are events recorded by the SANDWICH and Tibet-31N arrays, respectively. The blue (the event occurred on 30th April 2015 at UTC 10:45:05) and red (the event occurred on 26th June 2010 at UTC 05:30:19) stars denote the location of the events shown in Fig. 2 as examples. Abbreviations: BNS: Bangong-Nujiang suture; IYS: Indus-Yarlung suture; JRS: Jinsha River suture; PXR: Pumqu-Xianza rift; TYR: Tangra Yum Co rift; YGR: Yadong-Gulu rift. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

tivities to deep structures. Some studies show that the dip angle of the SIML increases from subhorizontal in the west to a steeper angle in the east (Li et al., 2008; Zhao et al., 2010), whereas others claim an overall subhorizontal Indian subduction beneath Tibet (e.g., Chen et al., 2017). A series of Cenozoic north-south-trending rifts and pairs of conjugate strike-slip faults have developed in the Tibetan Plateau (e.g., Armijo et al., 1986, 1989; Kapp and Guynn, 2004; Yin and Taylor, 2011), indicating east-west extension in the crust. The mechanism of these extensional structures is still unclear. Possible candidates are gravity collapse after a maximum uplift and crustal thickening of the plateau or delamination of the lithospheric mantle (England and Houseman, 1989; Molnar and Tapponnier, 1978). Recently, NS-trending rifts have been suggested to reflect SIML slab tearing (e.g., Chen et al., 2015; Liang et al., 2012; Tian et al., 2015; Yin, 2000). Variations in the SKS phase arrival time and shear wave birefringence across southern Tibet (Chen et al., 2015), alternately distributed low seismic velocity anomalies beneath the plateau (Li and Song, 2018; Liang et al., 2016a) and high electrical conductivity beneath southern Tibet (e.g., Wang et al., 2017; Wei et al., 2010) indicate that the SIML may be fragmented into several pieces, which would induce mantle upwelling and introduce east-west extension from the bottom of the lithosphere (Wu et al., 2019). The SIML slab tearing is consistent with the observations of mantle-derived helium isotopes in hot springs around southern Tibet and ultrapotassic magmatic rocks along the rifts (Guo et al., 2007; Hoke et al., 2000; Hou et al., 2004). However, the geodynamic process relating these geological extension structures to the proposed slab fragmentation is unclear. Currently, distinguishing why, where and how the SIML has been torn still needs more observations, especially LAB observations as the most straightforward evidence. Because the S-wave

receiver functions (SRFs) isolate S-to-P converted energy from Swave wavefields and avoid the contamination of multiple phases, SRFs are more suitable for detecting the LAB than P-wave receiver functions. Since the previous station coverage in south-central Tibet is insufficient to conduct a detailed SRF study, the previous stations (Tibet-31N array) and our newly deployed stations (SANDWICH array) in central Tibet are combined in this study. We obtain the SRFs using the teleseismic waveform data recorded by the two seismic arrays and image the morphology of the Indian and Tibetan LABs. 2. Data and methods The Seismic Array iNtegrated Detection for a Window of Indian Continental Head (SANDWICH, Liang et al., 2016b) is a 2-D broadband seismic array in central Tibet with a total of 61 stations distributed on both sides of the Bangong-Nujiang Suture (BNS). The interstation distance is ∼40 km (blue triangles in Fig. 1). More than 2000 teleseismic events with magnitudes greater than 5.3 were recorded by the SANDWICH from November 2013 - October 2016. We selected events with epicentral distances of 50◦ -120◦ to calculate SRFs. S phases were used in the epicentral distance range of 50◦ -90◦ , while SKS phases were used in the distance range of 90◦ -120◦ . To improve the resolution in the southern part of the research region, we added data recorded by the Tibet-31N array (55 stations), which was operated from September 2009 to November 2010 (red inverted triangles in Fig. 1) (Chen et al., 2015; Tian et al., 2015). Receiver function provide an effective mean of detecting underground discontinuities. Direct waves produce converted waves when they encounter velocity discontinuities. The converted waves can be considered as the convolution of the direct wave with the

JID:EPSL

AID:115954 /SCO

[m5G; v1.279; Prn:15/11/2019; 12:04] P.3 (1-9)

Z. Liu et al. / Earth and Planetary Science Letters ••• (••••) ••••••

3

Fig. 2. S-wave receiver functions of the two events (red and blue stars marked in the bottom-right inset of Fig. 1). The station names are listed on the left and the information about the events is shown at the top. The SRFs filled with gray and black were deconvolved with a Gaussian parameter of 1.0, and the blue or red unfilled SRFs were deconvolved with a Gaussian parameter of 0.6. The dashed rectangles mark the approximate arrival times of the obvious LAB phases.

effects of structure beneath the station. Using deconvolution of the direct wave component from the converted wave component, we can extract the receiver function, which represents the structure beneath the station (Ammon, 1991; Vinnik, 1977). In the SRF, the S-to-P converted phase (Sp) arrives earlier, while the multiple phases (SsSp, SpSp+SsPp) arrive later than the incoming S phase. Thus, using SRFs to detect the LAB can avoid contamination by multiple phases. We used a window of 90 s (60 s before and 30 s after the S/SKS wave arrival) to trim waveforms. A zero-phase 0.03-0.4 Hz bandpass filter was used to enhance the signal-to-noise ratio. The East and North component waveforms were rotated to radial and transverse components according to the theoretical back azimuths. The SRFs were calculated in the time domain (Ammon, 1991) with an iterative deconvolution (100 times), by removing the radial components from the vertical components. The amplitude of the converted phase in the receiver function depends on the velocity contrast across the discontinuity as well as the sharpness of the velocity discontinuity. A downward increase in velocity, such as that occurs at the Moho, generates a converted phase with negative amplitude in the SRF, while a decrease in velocity, such as that occurs at the LAB, produces a converted phase with positive amplitude. The polarities of the SRF phase are inverted to make them consistent with those of P-wave receiver function. We can

obtain images of such velocity discontinuities by stacking the amplitudes to the position where the conversion occurred after the time-depth conversion using a background velocity model (common conversion point (CCP) stacking) (Zhu, 2000). To assure the data quality, there are two stages of visual inspection process applied to both the original waveform data and the SRFs. First, the waveforms of same events were aligned according to direct S waves after the bandpass filter (0.03-0.4 Hz) to remove low signal-noise-ratio waveforms before calculating the SRFs. Then, the extracted SRFs of same stations were aligned to compare and remove poor ones. We selected a total of 316 events to calculate SRFs, including 252 events recorded by the SANDWICH array and 64 events recorded by the Tibet-31N array, and obtained a total number of 4595 SRFs, including 3664 from SANDWICH and 931 from Tibet-31N, respectively. The locations of events recorded by SANDWICH and Tibet-31N are displayed in the bottom-right inset of Fig. 1 with light blue pluses and light red crosses respectively. The SRFs of two example events for SANDWICH and Tibet-31N arrays are displayed in Fig. 2. The three-component waveforms and corresponding SRFs from stations GACO of SANDWICH and C12 of Tibet-31N are displayed in Fig. S1. More SRFs are shown in Figs. S2 and S3. Due to the filter response, side lobes may be generated accompanying the real signal during bandpass filtering. The arrival time of the peak amplitude of a side lobe varies with different fre-

JID:EPSL

4

AID:115954 /SCO

[m5G; v1.279; Prn:15/11/2019; 12:04] P.4 (1-9)

Z. Liu et al. / Earth and Planetary Science Letters ••• (••••) ••••••

Fig. 3. CCP stacking images for SRFs. (a) The green lines indicate the locations of the CCP stack profiles. The red pluses and blue crosses denote the locations of the piercing points at depths of 150 km and 200 km, respectively. (b) to (i) The CCP stack images along the profiles shown in (a). The negative signals are colored blue, while the positive signals are colored red. The dark green dashed lines indicate the Tibetan LAB, the blue dashed lines indicate the Indian LAB, and the thin orange dashed lines indicate an uncertain discontinuity. The shaded area has a stack number of less than 100. Abbreviations: BNS: Bangong-Nujiang suture; YGR: Yadong-Gulu rift.

quency bands. In contrast, the arrival time of a real signal remains stable in different frequency bands and depends on the velocity structure only if a zero-phase filter has been used. Since the amplitudes of SRFs are usually weak, distinguishing a true converted signal from an artificial side lobe is important. We calculated SRFs with different Gaussian parameters of 0.6 and 1.0 to test whether negative amplitudes are consistent in different bands. The aligned troughs of two SRF series with different Gaussian parameters help to differentiate the reliable converted phases as LAB signals in SRFs. The results show that the troughs of the SRFs calculated with two Gaussian parameters match well; thus, we identify the troughs as signals from negative velocity discontinuities rather than side lobes. On the other hand, a negative signal shouldn’t be a side lobe of a positive phase if they are separated from each other. An evaluation criterion could be the number of inflection points on the slope between the neighboring trough (negative phase) and peak (positive phase). If a peak (or trough) is a side lobe of its neighboring trough (or peak), no more than one inflection point exists on the slope between them (Fig. S4a), because two or more inflection points would separate the adjacent peak and trough. As a low-frequency SRF may merge the neighboring phases and smooth

the waveform, we also surveyed SRFs with higher-frequency content deconvolved with a Gaussian parameter of 1.0 (Fig. S4b). The trends of the slopes between the adjacent peaks and troughs of the SRFs show that the troughs are separated from their neighboring peaks, indicating that the negative LAB phase should be a real signal rather than a side lobe of the positive phase. 3. Results We obtained a total of 4595 SRFs, including 3664 from SANDWICH and 931 from Tibet-31N. The LAB structures along five approximately north-south-trending and three east-west-trending profiles were imaged with the CCP stacking technique (Zhu, 2000) using the IASP91 model as a reference (Kennett and Engdahl, 1991) (Fig. 3). The CCP images were smoothed with a bin of 20 km along the profiles, 1 km in vertical dimension and 70 km in width. In the CCP images, regions with stacked numbers of SRFs larger than 100 were highlighted. The hit-count maps of each profile are shown in figure S5. The Moho as the most prominent positive velocity discontinuity lies at a depth of ∼70 km (red signals, Fig. 3b-3i). Below the Moho, our results reveal several negative velocity discontinuities

JID:EPSL

AID:115954 /SCO

[m5G; v1.279; Prn:15/11/2019; 12:04] P.5 (1-9)

Z. Liu et al. / Earth and Planetary Science Letters ••• (••••) ••••••

(blue signals in Fig. 3b-3i), which may represent the LAB, generally consistent with previous results (Kumar et al., 2006; Zhao et al., 2010, 2011) (Fig. S6). Three east-west cross-sections (Fig. 3b-3d) show obvious variations in the negative velocity discontinuities from west to east, with an apparent step across the YGR. The strong negative signals are marked by dashed lines with different colors and projected to the approximately north-south cross-sections (Fig. 3e-3i). We interpret the continuous discontinuities as the Indian and Tibetan LABs. The Tibetan LAB is at a depth shallower than 150 km. An offset of the Tibetan LAB is observed at ∼90◦ E beneath the central part our research region (Fig. 3c, ∼32◦ N). The thickness of the Tibetan lithosphere decreases from ∼150 km in the west to ∼120 km in the east (Fig. 3c). The deeper signal beneath the Tibetan LAB could represent the Indian LAB. To the east of the YGR, the Indian LAB is at a depth of ∼200 km. The structure to the west of the YGR becomes more complex. In Figs. 3c, 3d, 3f, we can see that the discontinuity at ∼270 km extends southward much farther than the ones at ∼210 km. Previous studies, such as tomography (Li et al., 2008; Liang et al., 2016a) and anisotropy (Chen et al., 2015), indicated a high dip angle of the Indian slab. Thus, we identify the Indian LAB as a negative velocity discontinuity reaching ∼270 km depth (Figs. 3c and 3d). More E-W profiles (∼0.5◦ spacing in latitude) indicate that the discontinuity at ∼270 km west of YGR extends further southward than the one at ∼210 km (Fig. S7). The discontinuity at ∼210 km locates slightly north of BNS (thin orange dashed line in Fig. 3c). This discontinuity might be a local tectonic structure, such as the bottom boundary of the high velocity body beneath the northern Lhasa detected by Liang et al. (2016a) or less likely the Asian LAB from the north, because the discontinuity does not extend farther south and we do not know whether it extends farther north due to the low resolution. 4. Discussion 4.1. Influence of different velocity models used for the CCP stacking The lateral variation in the lithospheric structure may affect the CCP stacking because it could change the ray paths dramatically given the relatively large ray parameters. Therefore, the IASP91, as a horizontal layered model, might introduce some artifacts into the CCP images. We used two different seismic velocity models to obtain the CCP images of all SRFs to evaluate the velocity model dependence of the final LAB images. First, we used a 2-D model with different velocity structures for the areas east and west of the YGR (Figs. 4a and 4b) to evaluate the potential effect of the velocity variation across the YGR. The EW stepped velocity model was established using the depths of the velocity discontinuities from this study. The lithospheric thickness of the Tibetan Plateau was set to 115 km to the east and 150 km to the west of the YGR, and the LAB depth of the SIML was set to 200 km to the east and 275 km to the west of the YGR. The thickness of the SIML was fixed at 65 km. The S-wave velocities are from Zhao et al. (2011), and the V P / V S ratios are fixed at 1.80 and 1.75 for mantle and crust, respectively. The result shows that the discontinuities are similar to those imaged using the IASP91 model (Fig. S8). The hit-count maps are shown in Fig. S9. We then tested a 3-D hybrid velocity model combining a regional seismic velocity model from double-difference seismic travel-time tomography for depths shallower than 150 km (Xin et al., 2018) and body wave tomography for depths greater than 150 km (Liang et al., 2016a). Because the results of Liang et al. (2016a) provide only velocity perturbations, we smoothed the mean values of the seismic velocity at 150 km depth from Xin et al. (2018) and the IASP91 model at 400 km as local references to convert the velocity perturbations to absolute velocities. This model was

5

used to evaluate the influence from 3-D P- and S-wave velocity models, since the SRF is sensitive to the V P / V S ratio (Tian et al., 2005). The results show that there is not significant difference between the discontinuities obtained using the 1-D layered model and those using the 3-D model (Fig. S10). The hit-count maps are shown in Fig. S11. 4.2. Synthetic SRFs from different velocity models The prominent feature of the LABs is the depth offset across the YGR. Two synthetic SRFs were generated using the two 1-D velocity models for the areas east and west of the YGR (Fig. 4). The synthetic SRFs show structures similar to those of the CCP images. The direct and converted ray paths might pass through different velocity structures due to the large ray parameters. To test the influence of the lateral variable velocity model, synthetic SRFs were calculated based on two 2-D step models (Fig. 5) (Zhao et al., 2008) with or without a slab window. The ray paths of the synthetic SRFs were taken from the left (west) and right (east), as the models are homogeneous in the north-south direction. Then, the CCP stacking method was used to image the velocity discontinuities. The CCP images show that the major discontinuities can be recovered, but that the slab window cannot be recovered, which might be caused by diffraction of seismic waves at the slab margins or smearing effect during CCP stacking. Two models with and without a northward-dipping subducting Indian slab were also tested (Fig. 6) using the RF code from Frederiksen and Bostock (2000). The dip angle of the subducting slab was set to 20◦ . The velocity values were the same as in previous models. Synthetic SRFs from different back azimuths recorded at a certain station (the inverted triangle in Fig. 6) were calculated. The SRFs change dramatically with the back azimuth for the subducting slab model (Fig. 6c), especially for the events parallel to the slab dip direction. However, the events used in this study mostly came from the ESE direction nearly vertical to the dip direction. The waveforms plotted in red have back azimuths similar to those in the real data (Fig. 6). The synthetic test shows that the polarity of the negative amplitude does not inverse in the back-azimuth we used. These waveforms were used to conduct CCP stacking, and the input discontinuities can be recovered reasonably well with limited back azimuthal coverage as the real dataset. Thus, a dipping structure can be recovered using our observing system. The nearly horizontal Indian LABs in our results indicate horizontal SIML beneath the BNS. The positive phase in the synthetic test representing the interface at the top of the SIML is absent in the real data. The reason could be that the velocity model used to synthesize seismograms has a sharp discontinuity at the top of the slab while the real velocity boundary may be a gradual transition at the top of the SIML. The 1-D velocity model with a gradual interface verifies this hypothesis (Fig. S12). 4.3. Deep structure and dynamics inferred from the CCP images We observed a LAB offset of the SIML beneath the YGR. The depth of the bottom of the Indian slab decreases from ∼270 km in the west to ∼200 km in the east. This offset has a good correlation with the segmented gravity data (Hetenyi et al., 2016). The low S-wave velocity anomaly at depths of ∼170-250 km (Liang et al., 2016a) is spatially consistent with the LAB offset in our SRF images. Li and Song (2018) also proposed the tearing of the SIML from their Pn tomography result, and their proposed dip angle change across the YGR is consistent with our result. Our result implied that the tearing occurred at least in the central Lhasa terrane, consistent with the results of Chen et al. (2015), Bao et al. (2015) and Shi et al. (2015). Other tomography results

JID:EPSL

6

AID:115954 /SCO

[m5G; v1.279; Prn:15/11/2019; 12:04] P.6 (1-9)

Z. Liu et al. / Earth and Planetary Science Letters ••• (••••) ••••••

Fig. 4. Comparisons of the synthetic S-wave receiver functions with the real CCP images. Two 1-D S-wave velocity models were used for the west and east of the YGR (a and b). The S-wave velocities are taken from Zhao et al. (2011), and the V P / V S ratios are fixed at 1.75 in the crust and 1.80 in the mantle. The locations of CCP profiles (c), (e), (d) and (f) are the same as those in Fig. 3e, 3f, 3h and 3i. The synthetic waveforms in (c) and (e) are calculated using the S-wave velocity model in (a), and the waveforms in (d) and (f) are calculated using the velocity model in (b). The negative amplitude of the synthetic SRFs is filled with gray color.

(Liang et al., 2012; Agius and Lebedev, 2013; Nunn et al., 2014) also suggested that the Indian slab subducts at a shallower angle east of YGR than that in the west. We indicated that the LAB of the SIML is ∼200 km east of the YGR near the BNS, which is different from the interpretations of Kumar et al. (2006) and Zhao et al. (2011), who argued that the LAB belongs to Asia. In fact, the poor resolution beneath the Qiangtang terrane in Kumar et al. (2006) leads to a gap between the imaged LAB and the Asian lithosphere. In Zhao et al. (2011), the image of the 180-km-deep LAB is also connected to the Indian slab beneath the IYS. In contrast, the results of body wave tomography (Li et al., 2008; Liang et al., 2012; Tilmann et al., 2003), receiver functions (Mechie and Kind, 2013), and surface wave tomography (Ceylan et al., 2012; Chen et al., 2010) suggest that the

Indian lithosphere has reached the BNS or even farther north. Furthermore, P-wave receiver function images in Zhao et al. (2011) show that the top of the Asian mantle lithosphere (Moho) reaches a depth of ∼160 km with a high dip angle (Kind et al., 2002) beneath the Jinsha River suture (JRS) to the east along the GolmodLhasa highway, which is inconsistent with the speculation that the Asian slab is shallower than 160 km beneath the BNS and the Qiangtang terrane. The widespread Cenozoic volcanic rocks south of the JRS (Chung et al., 2005) indicate a broad and deep channel for asthenospheric upwelling beneath the northern Qiangtang terrane. Body wave tomography results also suggest no or limited subduction of the Asian lithosphere beneath the Songpan-Ganzi fold belt for eastern Tibet (Liang et al., 2012).

JID:EPSL

AID:115954 /SCO

[m5G; v1.279; Prn:15/11/2019; 12:04] P.7 (1-9)

Z. Liu et al. / Earth and Planetary Science Letters ••• (••••) ••••••

7

Fig. 5. CCP stacking images of synthetic S-wave receiver functions (b and d) based on two 2-D step models. (a) and (c) show two different 2-D step models. The seismic velocities in the models are the same as in Fig. 4. The Tibetan lithosphere is in yellow, while the Indian lithosphere is in blue. The receivers are located at the dark green triangles in (a) and (c). The positive and negative velocity discontinuities are marked as gray lines and dashed gray lines, respectively. The ray paths are from the left (west) and right (west) of the models.

Based on our SRF images and the discussion above, we propose a model to interpret our observations (Fig. 7). The LABs of the Tibetan Plateau and the SIML are offset beneath the YGR. The sudden change in the LAB depth for both the Indian lithospheric slab and the Tibetan lithosphere is direct and reliable evidence for the tearing of the subducting Indian slab and different degree of thickening for the plateau lithosphere on both sides of YGR. The tearing of the Indian lithosphere may provide a channel for asthenospheric upwelling (Chen et al., 2015; Liang et al., 2016a; Wu et al., 2019), which itself is beyond the resolvability of our SRF results. The upwelling of the asthenosphere may induce local mantle convection, which provides driving forces for the formation of the rift (Yin, 2000; Tian et al., 2015) and velocity contrast decrease. The rift penetrates the whole Tibetan Plateau and provides a channel for deep magma to migrate to the surface, which is consistent with the widespread occurrence of ultra-potassic volcanic rocks along the rifts (Guo et al., 2007). The tearing of the SIML may be caused by the variation in its inherited structures. Moho offsets are found among different blocks of the Indian plate (Hetenyi et al., 2016; Mandal et al., 2013), which is an assembly of multiple Precambrian cratonic blocks (Gahalaut and Kundu, 2012; Rao, 1973). Furthermore, several crustal ridges have been observed in northern India, which are subparallel to the rifts in Tibet (Gahalaut and Kundu, 2012; Rao, 1973). Thus, the variations in the Indian lithosphere may result in sufficient strain localization to cause the tearing of the SIML. Alternatively, the tearing may arise from the lateral variations in the geometry of the SIML. The arched nature of the Himalayas and the IYS might lead to different onset locations for the Indian subduction, which may also facilitate tearing of the SIML in the weak areas with the same mechanism as the slab tearing in the Calabrian and Apennine subduction systems (Rosenbaum et al., 2008). The negative amplitudes in northern Tibet is wider than that in the south (Figs. 3b-3d) indicating gradual changes of velocity rather than a shape discontinuity across the LAB. The gradual changes of velocity could be caused by a large region of low velocity due to upwelling asthenosphere. This might indicate that

the mantle lithosphere in northern Tibet has been totally removed (McNamara et al., 1997). 5. Conclusions Based on the CCP images of 4595 SRFs, offsets of the LAB in both the Tibetan and Indian lithosphere are observed beneath the YGR (Fig. 3c). The thickness of the Tibetan lithosphere is ∼30 km greater to the west of the YGR than to the east. Similarly, the Indian lithosphere has an offset of ∼70 km across the YGR and is deeper in the west. The proposed offset of the Indian lithosphere might indicate a variation in subduction angle from west to east beneath the Tibetan Plateau. Our results provide valuable evidence for lateral variations in the lithospheric structures and suggest tearing of the SIML along the YGR, which may provide a pathway for asthenospheric materials to the surface of the Earth. Acknowledgements Most figures were plotted using the Generic Mapping Tools (Wessel and Smith, 1998). The 1-D synthetic tests were conducted using the Computer Programs in Seismology (Herrmann, http:// www.eas.slu.edu/People/RBHerrmann/ComputerPrograms.html). We thank Prof. Liang Zhao, Prof. Zengxi Ge, Dr. Zhenbo Wu, Dr. Chao Lv, Dr. Lin Chen, Tingzi Li and Haiyan Yang for their kindly and helpful discussions. This research was supported by the National Key R&D Program of China (2016YFC0600301), the Strategic Priority Research Program (B) of the Chinese Academy of Sciences (Grant No. XDB03010700) and the National Natural Science Foundation of China (Grant Nos. 41604051, 41574048, 41574056, and 41804056). Seismic instruments were provided and maintained by the Seismic Array Laboratory, Institute of Geology and Geophysics, Chinese Academy of Sciences (IGGCAS). The receiver function SAC files are available in the supporting information of this paper. All raw data can be downloaded from the data center of the Seismic Array Laboratory, IGGCAS (http://www.seislab.cn/).

JID:EPSL

8

AID:115954 /SCO

[m5G; v1.279; Prn:15/11/2019; 12:04] P.8 (1-9)

Z. Liu et al. / Earth and Planetary Science Letters ••• (••••) ••••••

Fig. 6. Synthetic SRFs and their CCP stacking images from two 3-D velocity models. (a) and (b) show two different velocity models with a dipping structure or only horizontal layers. The values of the seismic wave velocities are the same as in Fig. 4. (c) and (d) show the SRFs from different back azimuths recorded by the stations located at the inverted triangles in (a) and (b), respectively. The red waveforms have back azimuths similar to those of real data. (e) and (f) show the CCP stacking images used the synthetic waveforms for two different models. The ray paths used to calculate the synthetic SRFs are the same as the real data. The gray lines and dashed gray lines in Fig. 6e and 6f show the locations of the discontinuities in Fig. 6a and 6b, respectively.

Appendix A. Supplementary material Supplementary material related to this article can be found online at https://doi.org/10.1016/j.epsl.2019.115954. References

Fig. 7. 3-D interpretative diagram showing the slab tearing model for the SIML with different subduction angles across the YGR. Upwelling of the asthenosphere through the tearing window may have fed the low seismic velocity in the uppermost mantle and contributed to the development of the rift. The southern part (without data coverage in this study) has been dimmed to a light color scheme and shaded by a semitransparent plate. The upward gradient color of the SIML indicates that the seismic velocity in the uppermost mantle varies gradually and the top of the subducted slab is not well detected.

Agius, M.R., Lebedev, S., 2013. Tibetan and Indian lithospheres in the upper mantle beneath Tibet: evidence from broadband surface-wave dispersion. Geochem. Geophys. Geosyst. 14, 4260–4281. Ammon, C.J., 1991. The isolation of receiver effects from teleseismic P-wave-forms. Bull. Seismol. Soc. Am. 81, 2504–2510. Armijo, R., Tapponnier, P., Mercier, J.L., Han, T.L., 1986. Quaternary extension in southern Tibet - field observations and tectonic implications. J. Geophys. Res. B, Solid Earth Planets 91, 13803–13872. Armijo, R., Tapponnier, P., Tonglin, H., 1989. Late Cenozoic right-lateral strike-slip faulting in southern Tibet. J. Geophys. Res. 94, 2787–2838. Bao, X.W., Song, X.D., Li, J.T., 2015. High-resolution lithospheric structure beneath Mainland China from ambient noise and earthquake surface-wave tomography. Earth Planet. Sci. Lett. 417, 132–141. Ceylan, S., Ni, J., Chen, J.Y., Zhang, Q., Tilmann, F., Sandvol, E., 2012. Fragmented Indian plate and vertically coherent deformation beneath eastern Tibet. J. Geophys. Res. 117, 11. Chen, M., Niu, F., Tromp, J., Lenardic, A., Lee, C.-T.A., Cao, W., Ribeiro, J., 2017. Lithospheric foundering and underthrusting imaged beneath Tibet. Nat. Commun. 8, 15659. Chen, Y., Badal, J., Hu, J.F., 2010. Love and Rayleigh wave tomography of the QinghaiTibet Plateau and surrounding areas. Pure Appl. Geophys. 167, 1171–1203.

JID:EPSL

AID:115954 /SCO

[m5G; v1.279; Prn:15/11/2019; 12:04] P.9 (1-9)

Z. Liu et al. / Earth and Planetary Science Letters ••• (••••) ••••••

Chen, Y., Li, W., Yuan, X., Badal, J.H., Teng, J.W., 2015. Tearing of the Indian lithospheric slab beneath southern Tibet revealed by SKS-wave splitting measurements. Earth Planet. Sci. Lett. 413, 13–24. Chung, S.L., Chu, M.F., Zhang, Y.Q., Xie, Y.W., Lo, C.H., Lee, T.Y., Lan, C.Y., Li, X.H., Zhang, Q., Wang, Y.Z., 2005. Tibetan tectonic evolution inferred from spatial and temporal variations in post-collisional magmatism. Earth-Sci. Rev. 68, 173–196. England, P., Houseman, G., 1989. Extension during continental convergence, with application to the Tibetan Plateau. J. Geophys. Res. 94, 17561–17579. Frederiksen, A.W., Bostock, M.G., 2000. Modelling teleseismic waves in dipping anisotropic structures. Geophys. J. Int. 141, 401–412. Gahalaut, V.K., Kundu, B., 2012. Possible influence of subducting ridges on the Himalayan arc and on the ruptures of great and major Himalayan earthquakes. Gondwana Res. 21, 1080–1088. Guo, Z.F., Wilson, M., Liu, J.Q., 2007. Post-collisional adakites in South Tibet: products of partial melting of subduction-modified lower crust. Lithos 96, 205–224. Hetenyi, G., Cattin, R., Berthet, T., Le Moigne, N., Chophel, J., Lechmann, S., Hammer, P., Drukpa, D., Sapkota, S.N., Gautier, S., Thinley, K., 2016. Segmentation of the Himalayas as revealed by arc-parallel gravity anomalies. Sci. Rep. 6, 10. Hoke, L., Lamb, S., Hilton, D.R., Poreda, R.J., 2000. Southern limit of mantle-derived geothermal helium emissions in Tibet: implications for lithospheric structure. Earth Planet. Sci. Lett. 180, 297–308. Hou, Z.Q., Gao, Y.F., Qu, X.M., Rui, Z.Y., Mo, X.X., 2004. Origin of adakitic intrusives generated during mid-Miocene east-west extension in southern Tibet. Earth Planet. Sci. Lett. 220, 139–155. Kapp, P., Guynn, J.H., 2004. Indian punch rifts Tibet. Geology 32, 993–996. Kennett, B.L.N., Engdahl, E.R., 1991. Traveltimes for global earthquake location and phase identification. Geophys. J. Int. 105, 429–465. Kind, R., Yuan, X., Saul, J., Nelson, D., Sobolev, S.V., Mechie, J., Zhao, W., Kosarev, G., Ni, J., Achauer, U., Jiang, M., 2002. Seismic images of crust and upper mantle beneath Tibet: evidence for Eurasian plate subduction. Science 298, 1219–1221. Kumar, P., Yuan, X.H., Kind, R., Ni, J., 2006. Imaging the colliding Indian and Asian lithospheric plates beneath Tibet. J. Geophys. Res. 111, 11. Li, C., Van der Hilst, R.D., Meltzer, A.S., Engdahl, E.R., 2008. Subduction of the Indian lithosphere beneath the Tibetan Plateau and Burma. Earth Planet. Sci. Lett. 274, 157–168. Li, J.T., Song, X.D., 2018. Tearing of Indian mantle lithosphere from high-resolution seismic images and its implications for lithosphere coupling in southern Tibet. Proc. Natl. Acad. Sci. USA 115, 8296–8300. Liang, C.T., Song, X.D., 2006. A low velocity belt beneath northern and eastern Tibetan Plateau from Pn tomography. Geophys. Res. Lett. 33, 5. Liang, X.F., Sandvol, E., Chen, Y.J., Hearn, T., Ni, J., Klemperer, S., Shen, Y., Tilmann, F., 2012. A complex Tibetan upper mantle: a fragmented Indian slab and no southverging subduction of Eurasian lithosphere. Earth Planet. Sci. Lett. 333, 101–111. Liang, X.F., Chen, Y., Tian, X.B., Chen, Y.J., Ni, J., Gallegos, A., Klemperer, S.L., Wang, M.L., Xu, T., Sung, C.Q., Si, S.K., Lan, H.Q., Teng, J.W., 2016a. 3D imaging of subducting and fragmenting Indian continental lithosphere beneath southern and central Tibet using body-wave finite-frequency tomography. Earth Planet. Sci. Lett. 443, 162–175. Liang, X.F., Tian, X.B., Zhu, G.H., Wu, C.L., Duan, Y.H., Li, W., Zhou, B.B., Zhang, M.H., Yu, G.P., Nie, S.T., Wang, G.C., Wang, M.L., Wu, Z.B., Liu, Z., Guo, X., Zhou, X.P., Wei, Z., Xu, T., Zhang, X., Bai, Z.M., Chen, Y., Teng, J.W., 2016b. SANDWICH: a 2D broadband seismic array in central Tibet. Seismol. Res. Lett. 87, 864–873. Mandal, B., Sen, M.K., Rao, V.V., 2013. New seismic images of the Central Indian Suture Zone and their tectonic implications. Tectonics 32, 908–921. McNamara, D.E., Walter, W.R., Owens, T.J., Ammon, C.J., 1997. Upper mantle velocity structure beneath the Tibetan Plateau from Pn travel time tomography. J. Geophys. Res., Solid Earth 102, 493–505. Mechie, J., Kind, R., 2013. A model of the crust and mantle structure down to 700 km depth beneath the Lhasa to Golmud transect across the Tibetan plateau as derived from seismological data. Tectonophysics 606, 187–197. Molnar, P., Tapponnier, P., 1978. Active tectonics of Tibet. J. Geophys. Res. 83, 5361–5375. Mulch, A., Chamberlain, C.P., 2006. The rise and growth of Tibet. Nature 439, 670–671. Nunn, C., Roecker, S.W., Priestley, K.F., Liang, X.F., Gilligan, A., 2014. Joint inversion of surface waves and teleseismic body waves across the Tibetan collision zone: the fate of subducted Indian lithosphere. Geophys. J. Int. 198, 1526–1542.

9

Rao, M.B.R., 1973. Subsurface geology of Indo-Gangetic plains. J. Geol. Soc. India 14, 217–242. Rosenbaum, G., Gasparon, M., Lucente, F.P., Peccerillo, A., Miller, M.S., 2008. Kinematics of slab tear faults during subduction segmentation and implications for Italian magmatism. Tectonics 27, 16. Shi, D.N., Wu, Z.H., Klemperer, S.L., Zhao, W.J., Xue, G.Q., Su, H.P., 2015. Receiver function imaging of crustal suture, steep subduction, and mantle wedge in the eastern India-Tibet continental collision zone. Earth Planet. Sci. Lett. 414, 6–15. Styron, R., Taylor, M., Okoronkwo, K., 2010. Database of active structures from the Indo-Asian collision. Trans. Am. Geophys. Union 91, 181–182. Tapponnier, P., Xu, Z.Q., Roger, F., Meyer, B., Arnaud, N., Wittlinger, G., Yang, J.S., 2001. Oblique stepwise rise and growth of the Tibet plateau. Science 294, 1671–1677. Tian, X.B., Chen, Y., Tseng, T.L., Klemperer, S.L., Thybo, H., Liu, Z., Xu, T., Liang, X.F., Bai, Z.M., Zhang, X., Si, S.K., Sun, C.Q., Lan, H.Q., Wang, E.C., Teng, J.W., 2015. Weakly coupled lithospheric extension in southern Tibet. Earth Planet. Sci. Lett. 430, 171–177. Tian, X.B., Wu, Q.J., Zhang, Z.J., Teng, J.W., Zeng, R.S., 2005. Joint imaging by teleseismic converted and multiple waves and its application in the INDEPTH-III passive seismic array. Geophys. Res. Lett. 32, 4. Tilmann, F., Ni, J., Team, I.I.S., 2003. Seismic imaging of the downwelling Indian lithosphere beneath central Tibet. Science 300, 1424–1427. Vinnik, L.P., 1977. Detection of waves converted from P to SV in mantle. Phys. Earth Planet. Inter. 15, 39–45. Wang, G., Wei, W.B., Ye, G.F., Jin, S., Jing, J.N., Zhang, L.T., Dong, H., Xie, C.L., Omisore, B.O., Guo, Z.Q., 2017. 3-D electrical structure across the Yadong-Gulu rift revealed by magnetotelluric data: new insights on the extension of the upper crust and the geometry of the underthrusting Indian lithospheric slab in southern Tibet. Earth Planet. Sci. Lett. 474, 172–179. Wei, W.B., Jin, S., Ye, G.F., Deng, M., Jing, J.E., Unsworth, M., Jones, A.G., 2010. Conductivity structure and rheological property of lithosphere in southern Tibet inferred from super-broadband magnetotelluric sounding. Sci. China Earth Sci. 53, 189–202. Wessel, P., Smith, W.H.F., 1998. New, improved version of generic mapping tools released. Eos Trans. AGU 79, 1. Wu, C., Tian, X., Xu, T., Liang, X., Chen, Y., Taylor, M., Badal, J., Bai, Z., Duan, Y., Yu, G., Teng, J., 2019. Deformation of crust and upper mantle in central Tibet caused by the northward subduction and slab tearing of the Indian lithosphere: new evidence based on shear wave splitting measurements. Earth Planet. Sci. Lett. 514, 75–83. Xin, H., Zhang, H., Kang, M., He, R., Gao, L., Gao, J., 2018. High-resolution lithospheric velocity structure of continental China by double-difference seismic travel-time tomography. Seismol. Res. Lett. 90, 229–241. Yin, A., 2000. Mode of Cenozoic east-west extension in Tibet suggesting a common origin of rifts in Asia during the Indo-Asian collision. J. Geophys. Res. 105, 21745–21759. Yin, A., Taylor, M.H., 2011. Mechanics of V-shaped conjugate strike-slip faults and the corresponding continuum mode of continental deformation. Geol. Soc. Am. Bull. 123, 1798–1821. Zhang, Z.J., Klemperer, S.L., 2005. West-east variation in crustal thickness in northern Lhasa block, central Tibet, from deep seismic sounding data. J. Geophys. Res., Solid Earth 110, B09403. https://doi.org/10.1029/2004jb003139. Zhao, J.M., Yuan, X.H., Liu, H.B., Kumar, P., Pei, S.P., Kind, R., Zhang, Z.J., Teng, J.W., Ding, L., Gao, X., Xu, Q., Wang, W., 2010. The boundary between the Indian and Asian tectonic plates below Tibet. Proc. Natl. Acad. Sci. USA 107, 11229–11233. Zhao, L., Wen, L.X., Chen, L., Zheng, T.Y., 2008. A two-dimensional hybrid method for modeling seismic wave propagation in anisotropic media. J. Geophys. Res., Solid Earth 113, 21. Zhao, W.J., Kumar, P., Mechie, J., Kind, R., Meissner, R., Wu, Z.H., Shi, D.A., Su, H.P., Xue, G.Q., Karplus, M., Tilmann, F., 2011. Tibetan plate overriding the Asian plate in central and northern Tibet. Nat. Geosci. 4, 870–873. Zhu, L.P., 2000. Crustal structure across the San Andreas Fault, southern California from teleseismic converted waves. Earth Planet. Sci. Lett. 179, 183–190.