The splitting 660 km discontinuity associated with lithospheric delamination in the northern part of the North-South seismic zone, China

The north-south seismic zone (NSSZ) is a destructive zone of large-scale earthquakes in China, and the earthquake mechanism associated with deep structures remains unclear. Previous studies have indicated that lithospheric delamination or absence of lithospheres in the western part of the NSSZ may facilitate the eastern extrusion of the Tibetan Plateau and lead to stress accumulation and release. However, the deep process of lithospheric delamination needs to be further clariﬁed. In this study, I collect abundant high-quality teleseismic data recorded by permanent seismic stations and perform common conversion point (CCP) stacking of receiver functions in the north part of the NSSZ. The results show that lithospheric delamination might result in the splitting 660 km discontinuity and a thickening region of the mantle transition zone (MTZ).

. Portions of the model are not shown where the recovery from the input velocity model is below 20%.(He et al., 2014).

Introduction
The north-south seismic zone (NSSZ) on the Chinese continent is an important earthquake-prone region.From historically documented records, approximately 1/3 of the larger earthquakes (Ms>7) in China occurred in this area; moreover, the zone is a boundary between the highlands in the western part and lowlands in the eastern part (Li et al., 2006) and is also a zone of strong lateral gradient in gravity anomaly and crustal thickness (He et al., 2014;Zhang et al., 2013).
The study region located in the northern part of the NSSZ includes the Alashan terrane, Qilian orogenic belt, Ordos terrane, Songpan-Ganzi terrane and Western Qinling orogenic belt (Fig. 1) (Zhang et al., 2013;Zhang and Wang, 2009), which have undergone a complex tectonic evolution process (Zhao et al., 2010;Santosh, 2010).The Ordos terrane is a large intracratonic compressional basin (Wang et al., 2005;Li, 1996) that is affected by the collision between the South and North China Blocks in the Triassic (Ratschbacher et al., 2000;Ames et al., 1996).In the western part of this region, the Alxa terrane, Qilian terrane and Qaidam terrane assembled and collided during the Late Ordovician to Devonian (Xu et al., 2006).Specifically, the eastward extrusion of the Tibetan Plateau due to the collision between the Indian and Eurasian Plates in the Cenozoic might lead to the formation of an earthquake-prone zone in this area (Deng et al., 2003).
In this area, a number of geophysical studies have been carried out, such as tomography (e.g., Liu et al., 1989;He and Santosh, 2017), shear-wave splitting (Wang et al., 2008), receiver functions (e.g., He et al., 2014) and deep seismic sounding (e.g., Li et al., 2002;Gao et al., 2006).A receiver function study revealed that delamination might result in variations in the crustal structure and topography of the mantle transition zone (He et al. 2014).A tomographic study revealed high-velocity anomalies in the mantle transition zone, which might be linked with the delamination of lithospheric delamination (He and Santosh, 2020).However, the deep processes that are associated with the delamination in this area need to be further elucidated, which may be played in a key role for the formation of the earthquake-prone region.
To understand the delamination model in this area, I collected abundant highquality teleseismic data recorded by permanent seismic stations (Fig. 1

Data and method for the CCP stacking of receiver functions
A total of 982 teleseismic events were collected in this area, which were recorded by 214 permanent seismic stations of the China Seismic Network from 2007 to 2014 (Fig. 1, insert figure of left panel).The events were limited to Ms >5.8, and the earthquake epicentral distances ranged from 30°to 90°for individual event-station pairs.The raw waveforms with a 50 Hz or 100 Hz sampling rate were cut from 15 s before to 200 s after the P-wave arrival time and filtered using a Butterworth bandpass filter ranging from 0.1 to 0.2 Hz; waveform sample rates were decimated to 0.1 s.The waveform cross-correlation technique (VanDecar and Crosson, 1990) was employed to select consistent waveform data (an example can be found in Fig. S1).A modified frequency domain deconvolution using a Gaussian factor equal to 1 and a water level equal to 0.01 was employed to extract the receiver functions (Zhu and Kanamori, 2000).Finally, 13335 high-quality receiver functions (for example, please see Fig. S2), whose Ps phase and PPms phase can be clearly seen, were extracted and used for CCP stacking of receiver functions.This dataset is far larger than any dataset used for previous receiver function studies in this area.
A CCP stacking of receiver functions (e.g., Dueker and Sheehan, 1997;Eagar et al., 2010;Zhu, 2000) is used to image the MTZ of the northern part of the NSSZ.The piercing points at depths of 410 and 660 km are calculated by using the IASP 91 1-D velocity model (Kennett et al., 1991) (Fig. S3).Spherical coordinates are used to calculate the Ps-P differential time T Ps (Ps: PS converted wave, P: P-wave, T Ps : travel-time from Ps phase to P phase) (Eagar et al., 2010).The effect of velocity heterogeneities in the upper mantle is removed by a global 3-D S-and P-wave velocity model (Lu et al., 2019) and a local 3-D P-wave velocity model (He and Santosh, 2017).Lateral grid intervals of 0.5°and depth intervals of 1 km are designed for the CCP stacking of receiver functions, and the migrated receiver functions are searched within a radius of 75 km (Xu et al., 2018).

Results
Based on previous tomography (He and Santosh, 2017), large high-velocity anomalies and large low-velocity anomalies (Hv1 and Lv1) are present at depths of 70-200 km (Fig. S4).Hv1 might be the lithospheric root of the Ordos Basin; Lv1 may be a manifestation of absence of the lithosphere root in the western part of the study region.At depths of 300-700 km, many high-velocity anomalies that may be associated with lithospheric delamination during different periods (He and Santosh, 2017) are designated Hv2.Three profiles of P-wave velocity perturbation and CCP stacking of receiver functions were created simultaneously along the same location (location of profiles, see left panel of Fig. 1).Fig. 2. Latitudinal profiles of the P-wave velocity perturbation (a, b and c) (He and Santosh, 2017) and Latitudinal profiles of CCP stacking of receiver functions (a, b and c) (Location of profiles, see left panel of Fig. 1), The bootstrapping method is used to calculate the stacked amplitudes (resampling 2000 times in the dataset), and the 95% confidence level is used to calculate the final mean receiver functions (middle lines).Blue and black horizontal lines: depths of 410 km and 660 km.
The results show that a high-velocity anomaly appears at the MTZ (Fig. 2, upper panel of profiles a, b and c).CCP stacking profiles of receiver functions indicate that the topographic variation of the 410 km discontinuity is not obvious, which is basically consistent with the global average depth (Houser et al., 2008;Flanagan and Shearer, 1999).In contrast, the topography of the 660 km discontinuity shows amazing variation at profiles a and b (Fig. 2, lower panel of profiles a and b, rectangular region); the 660 km discontinuity splits into multiple interfaces, and its amplitude obviously becomes weak.The splitting interface of the 660 km discontinuity is distributed at depths from 630 to 800 km.Fig. 3. Profiles of CCP stacking of receiver functions overlap with those of the P-wave velocity perturbation (He and Santosh, 2017).
I overlaped two kinds of profiles (Fig. 3).The results indicate that the highvelocity anomaly (Hv2) almost overlaps the 660 km discontinuity in profiles a and b, whereas the high-velocity anomaly (Hv2) is basically above the 660 km discontinuity in profile c, whose variations are relatively small.Fig. 4. Topography of the 410 km and 660 km discontinuities and the thickness of the MTZ.The bootstrapping method is used to calculate the stacked amplitudes (resampling 2000 times in the dataset), and the 95% confidence level is used to calculate the final mean receiver functions; the number of stacking amplitude points is greater than 10 (a, b, c and d).Depths of the 410 km discontinuity (a) and 660 km discontinuity (c) and the thickness of the MTZ (e), which have been corrected on the basis of a global 3-D local velocity model (Lu et al., 2019).Depths of the 410 km discontinuity (b) and 660 km discontinuity (d) and the thickness of the MTZ (f), which have been corrected on the basis of a local 3-D local velocity model (He and Santosh, 2017), and S-wave velocity is inferred from the AK135 1-D velocity model (Kennett et al., 1991).(Kemp et al., 2019).Note that the depth stacks use a 3D model to convert from time to depth, while the predictions for the slowness stack use AK135 1-D velocity model (Kennett et al., 1991).
When conversions and surface multiples arrive at similar times, two kinds of wave interference may occur.Due to surface multiples with positive slowness and conversions with negative slowness relative to the direct P-wave, I performed slowness stacks to distinguish two kinds of waves.These stacks are created by shifting all the receiver functions in time to a common epicentral distance of 60°using relative slowness values between −0.8 and +0.8 s/deg compared with the direct P-wave slowness and then stacking the shifted receiver functions for each of these slowness values (Kemp et al., 2019).In slowness-time space, a bullseye shape shows negative and positive coherent amplitude arrivals.I used the AK135 1-D velocity model (Kennett et al., 1991) to calculate predicted lines and positions for multiples and Pds converted phases.The results show that the X-discontinuity, 410-km discontinuity and 660-km discontinuity are located at depths of 324, 411 and 667 km, respectively.There are no conversions or surface multiples of similar arrival times (Fig. 6).

Splitting 660 km discontinuity
The MTZ is characterized by abrupt changes in seismic velocities, which divide the Earth's upper and lower mantle.Based on seismic investigations, major global seismic velocity discontinuities appear at depths of approximately 410 and 660 km, which is the boundary of the top and bottom of the MTZ (van Stiphout et al., 2019).Generally, these discontinuities are considered polymorphic olivine phase change systems (e.g., Ringwood et al. 1975;Ito and Takahashi 1989;Katsura & Ito 1989).410 is linked with a phase transition from olivine to wadsleyite with a positive temperature-pressure gradient (or Clapeyron slope) (e.g., Katsura & Ito 1989).In contrast, 660 km is associated with a phase transition from ringwoodite into perovskite and magnesiowüstite with a negative Clapeyron slope (e.g., Ringwood et al. 1975;Ito and Takahashi 1989); when the temperature decreases, the topography of the 410 km discontinuity will become shallow.The non-olivine or garnet-dominant system occurs in a relatively hightemperature environment that may be connected with a positive Clapeyron slope (e.g., Hirose 2002;Yu et al. 2011); when the temperature decreases, the topography of the 660 km discontinuity will deepen.
A previous seismic study reported two discontinuities around the 660 km discontinuity or splitting 660 km discontinuity, suggesting that both phase transitions are occurring (Andrews and Deuss, 2008).The splitting 660 km discontinuity may be associated with both phases occurring in this area.However, the splitting 660 km discontinuity identified in this study obviously is not two discontinuities (Fig. 2, profiles a and b), and it may be three or four interfaces.Moreover, the amplitude of the 660 km discontinuity obviously becomes weak, which was first discovered in global seismic observations.

Lithospheric delamination
The lower continental crust, which might be a weak decoupling interface (Göğüş and Ueda, 2018), is suggested to provide a bridge between the felsic upper crust and the mafic/ultramafic upper mantle (Rudnick and Fountain, 1995;Christensen and Mooney, 1995).Under compression, the high-density lithosphere/lower crust might experience viscous dripping or convective removal due to Rayleigh-Taylor/gravitational instability (discrete high-velocity body delamination) or delamination along the weak decoupling interface (plate-like delam-ination) (Kay and kay, 1993;Göğüş et al., 2017) and contribute to velocity inhomogeneity in the mantle and variations in crustal composition.
A previous receiver function study indicated that the northern part of the NSSZ is not only a crustal thickness gradient zone (Fig. S5, a, white arrow: crustal thinning direction) but also a predominantly felsic lower crust induced by mafic/ultramafic lower crustal/lithospheric delamination (Fig. S4, b, rectangle region) (He et al., 2014).Tomographic images also indicate lithospheric absence in the western part of the study region induced by lithospheric delamination (Fig. S4, 70-200 km depth sections, Lv1) and several high-velocity anomalies (Hv2) at depths of 300-700 km, which might be associated with lithospheric delamination in this area (He and Santosh, 2017).Based on these images, it is suggested that a large amount of the lithosphere has been delaminated into the MTZ.Generally, it is widely accepted that the high-velocity structure or delaminated lithosphere in the MTZ is a cold domain (Foulger, 2012), which might result in a shallowing topography of the 410 km discontinuity and a deepening topography of the 660 km discontinuity as well as the thickening of the MTZ in the center part of the this study (Fig. 5).

Conclusions
1.The high-velocity anomalies in the mantle transition may lead to the 660 km discontinuity splitting into multiple discontinuities, which is worth further study.
2. Delaminated lithosphere in the MTZ resulted in thickening of the MTZ.
3. The results identified in this study support the notion that the lithosphere delaminated into the MTZ.

Fig
Fig. S1 P-wave velocity perturbation at different depths (He and Santosh, 2017).Portions of the model are not shown where the recovery from the input velocity model is below 20%.

Fig. S2 .
Fig. S2.Piercing points at depth of 410 km (a) and at depth of 660 km (b).Red points: piercing points at depths of 410 km and 660 km.Blue lines: profiles of P-wave velocity perturbation and CCP stacking profiles of receiver functions.

Fig. S4 .
Fig. S4.Crustal thickness (a) and Vp/Vs ratio (b) in the northern part of the NSSZ.Yellow triangles: seismic stations.Lower Vp/Vs ratios in the northern part of the NSSZ are generated by a deep process of lithospheric delamination(He et al., 2014).

Fig
Fig. S5 Vertical component at each seismic station for the 20120202133440 event.

Fig
Fig. S6 R and T receiver functions for the QH.LED station.
, insert figure of left panel) and carried out CCP stacking of receiver functions.The results indicate a splitting 660 km discontinuity and thickening of the MTZ beneath the northern part of the NSSZ, which might have been generated by the deep process of lithospheric delamination.

Fig. 1 .
Fig. 1.Left panel: location of the study region.Insert figure of the left panel: distribution of seismic events used in this study.Right panel: tectonic framework and distribution of seismic stations, blue lines: profiles of the Pwave velocity perturbation and CCP stacking of receiver functions.White lines of both panels: geological boundaries.
Fig. 5. Depth (left panel) and slowness (right panel) stacks of all 13335 receiver functions.The average amplitude at depth (solid black line) is plotted along with the lines reflecting 2 standard errors (dashed black line).Arrows indicate the significant positive peaks, with the depth in kilometres and individual symbols.Predicted lines for the conversion and multiple phases in slowness/time space using AK135 1-D velocity model are shown as Pds, PPvds (multiple) and PPvdp (multiple).The symbols indicate predicted arrivals for the corresponding peaks in the depth stacks: 324 km -orange Square, 411 km -green circle, 667 km -purple triangle for the low/high-frequency stacks(Kemp et al., 2019).Note that the depth stacks use a 3D model to convert from time to depth, while the predictions for the slowness stack use AK135 1-D velocity model(Kennett et al., 1991).