Crust and upper mantle structure of East Asia from ambient noise and earthquake

Mengkui Li, Xiaodong Song, Jiangtao Li and Xuewei Bao

1 School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China

2 Hebei Hongshan National Observatory on Thick Sediments and Seismic Hazards, Beijing 100871, China

3 School of Earth and Space Sciences, Peking University, Beijing 100871, China

4 School of Earth Sciences, Zhejiang University, Hangzhou 310027, China

ABSTRACT The complex tectonic background of East Asia makes it an ideal region for investigating the evolution of the continental lithosphere, for which high-resolution lithospheric structural models are essential. In this study, we measured Rayleigh-wave phase-velocity dispersion curves at periods of 10-120 s and group velocity dispersion curves at periods of 10-140 s using event records from more than 1,000 seismic stations in and around China. By jointly inverting new and previously published dispersion data from ambient noise and earthquakes, we developed a high-resolution shear-wave velocity model down to a depth of ～300 km beneath East Asia. Our model revealed heterogeneous lithospheric structures beneath East Asia, and allowed us to investigate the velocity structure of the entire lithosphere. We also derived crustal and lithospheric thickness models from the three-dimensional (3D) shear-wave model, revealing strong spatial heterogeneity and a general thinning trend of lithospheric thickness from west to east across the study region. Overall, our models reveal important lithospheric features beneath East Asia and provide a valuable baseline dataset for understanding continental-scale dynamics and evolution.

Keywords: East Asia; Rayleigh wave tomography; shear-wave velocity model; lithospheric thickness

1. Introduction

Our study region consists of Chinese mainland and the eastern marginal seas in East Asia including the Japan Sea,Bohai Sea, Yellow Sea, East China Sea, Taiwan Strait, and South China Sea (Figure 1). This region has drawn much attention because of its various lithospheric styles and complex tectonic background (e.g., Zhu JS et al., 2002;Liang CT et al., 2004; Priestley et al., 2006; Li YH et al.,2013; Pandey et al., 2014; Bao XW et al., 2015; Chen M et al., 2017). For example, the ongoing deformation of this region has been mainly affected by the continental collision of the Indian and Eurasian Plates to the west (e.g.,Tapponnier et al., 1986; Yin A and Harrison, 2000) and the subduction of the oceanic Pacific and Philippine Plates to the east and southeast (e.g., Northrup et al., 1995;Lallemand et al., 2001; Schellart and Lister, 2005). The collision of the Indian and Eurasian Plates has created the highest, largest and most active area of the Tibetan Plateau(TP) in the west, and in eastern China, the subduction of the Pacific and Philippine Plates has caused the development and widespread extension of the western Pacific marginal basins and the Mesozoic basins (e.g., Ren JY et al., 2002). The subduction of the Pacific Plate has also induced widespread asthenosphere upwelling, largescale volcanic eruptions, and strong tectonic activity in eastern China. Such complex tectonic processes have caused strong lithospheric heterogeneities in East Asia,making it an ideal region for investigating the evolution of the continental lithosphere. Importantly, detailed knowledge of the thickness and velocity structures of the lithosphere in East Asia is crucial for understanding continentalscale tectonics and geodynamic processes.

Figure 1. Topographic map of the study area showing the main tectonic units in the East Asia region, and the distributions of seismic stations (blue triangles) and earthquakes (blue dots). The thick blue lines represent the boundaries of geological units (from Liang CT et al., 2004, which is modified from the original map by Zhang ZM et al., 1984). The abbreviations are:Jungger Basin (JGB), Qaidam Basin (QB), Sichuan Basin (SB), Ordos Block (OB), Songliao Basin (SLB), North China Basin(NCB), Qiang-Tang Terrane (QT), Lhasa Terrane (LT), Himalaya Terrane (HT), Songpan-Garzê Terrane (SGT), Qilian Orogen (QLO), Qinling Orogen (QO), Dabie-Sulu Orogen (DSO), Changbai Mountain Range (CBM), Yangtze Craton (YZC),North China Craton (NCC), South China Block (SCB), Cathaysia Fold Belt (CFB), Taiwan Island (TW), Taiwan Strait (TS),North South Gravity Line (NSGL), Bohai Sea (BHS), Yellow Sea (YS), East China Sea (ECS), South China Sea (SCS).

Seismic tomography (including body and surface wave tomography) has been widely used to explore the lithospheric velocity structure beneath East Asia at different scales (e.g., Shapiro and Ritzwoller, 2002; Liang CT et al., 2004; Priestley et al., 2006; Zheng SH et al.,2008; Sun XL et al., 2010; Li YH et al., 2013; Xu Z et al.,2013; Huang ZC et al., 2014; Pandey et al., 2014; Bao XW et al., 2015; Shen WS et al., 2016; Chen M et al., 2017;Yang ZG and Song XD, 2019; Peng J et al., 2020).Continental-scale teleseismic body wave tomography provides three-dimensional (3D) P and S velocity images,which have revealed large-scale features in East Asia,including relatively low velocities in eastern China and high velocities underneath some basins. However, because of its steep ray path requirements, teleseismic body wave tomography typically has a low vertical resolution, which limits its detection of small-scale features. Xin HL et al.(2019) constructed a lithospheric velocity model for China continent using double-difference seismic travel time tomography, which was updated by incorporating surface wave dispersion data (Han SC et al., 2021). These authors report bothvPandvSstructures with a high resolution,especially in eastern China; however, the maximum depth of their model was approximately 150 km, providing structural information above the uppermost mantle only.

Surface wave tomography is an important tool for imaging lithospheric structures because of the depth sensitivity of surface waves with different frequencies.Many earthquake-based surface wave tomography studies at both continental and regional scales have been conducted in East Asia (e.g., Zhu JS et al., 2002; Li YH et al., 2013; Pandey et al., 2014). These studies reveal significant lateral variations in lithospheric structures in different geological units, yet the resolutions provided by these studies may be significantly influenced by the distribution of earthquakes and stations. The revolutionary ambient noise tomography (ANT) technique, which is based on the idea that Green′s function between two receivers can be extracted from the cross-correlation of the diffuse field (noise interferometry) at the receivers, has been widely used to image the structure of the crust and uppermost mantle (e.g., Shapiro et al., 2005; Yao HJ et al.,2006; Zheng SH et al., 2008). Recent ANT studies at varying resolutions have also been conducted beneath East Asia, mainly in Chinese mainland and the eastern marginal seas (e.g., Zheng SH, 2008; Sun XL et al., 2010; Bao XW et al., 2015; Wang Q et al., 2017; Shen WS et al., 2016;Yang ZG and Song XD, 2019; Peng J et al., 2020). Most ANT studies have focused on Rayleigh waves; however,shear-wave models developed from the Love waves extracted from noise interferometry have also been proposed (e.g., Yang ZG and Song XD, 2019). Such studies have revealed the structure of the crust and upper mantle in greater detail. However, because ambient noise interferometry can only provide dispersion information over relatively short periods, ANT can only be applied to reveal the structure at relatively shallow depths (Bao XW et al., 2015; Peng J et al., 2020). One recent study (Zhang XZ et al., 2022) systematically compared several seismic tomography models of the Chinese continental lithosphere and revealed large discrepancies in the existing velocity models both in absolute values and perturbation patterns.This indicates that there is a significant opportunity for improving the seismic models of this region.

Earthquake-based two-station (TS) analysis, which measures the dispersion curves from the empirical Green’s functions extracted from the cross-correlation of the surface-wave chains generated by the same earthquake event, can provide dispersion information over longer periods (e.g., Yao HJ et al., 2006, 2008; Li YH et al.,2013; Chen HP et al., 2014; Wang WL et al., 2014).Compared with ANT, earthquake-based TS analysis can be applied to greater depths (e.g., Yao HJ et al., 2008; Li MK et al., 2018; Ojo et al., 2018). Thus, the joint analysis of dispersion data from both noise interferometry and earthquake-based cross-correlation (e.g., Xu Z et al., 2013;Bao XW et al., 2015) can increase path coverage and extend the period range, enabling the construction of lithospheric S-wave velocity models with higher lateral and vertical resolutions.

In this study, using event data from more than 1,000 seismic stations (Figure 1) from multiple national, global,permanent and temporary networks in and surrounding Chinese mainland (i.e., the China Earthquake Administration’s CEArray; the Incorporated Research Institutions for Seismology Data Management Center’s (IRIS DMS)Portable Array Seismic Studies of the Continental Lithosphere, PASSCAL; and the Global Seismographic Network, GSN), we systematically measured Rayleigh-wave phase-velocity dispersion curves at periods of 10-120 s and group velocity dispersion curves at periods of 10-140 s based on the earthquake-based TS method. We then combined the new dispersion measurements with three previously published datasets from Bao XW et al. (2015),Wang Q et al. (2017) and Ekström (2011) to construct a 3D lithospheric S-wave velocity model down to a depth of approximately 300 km in East Asia. The combined dataset provides a relatively uniform and perhaps the most complete surface wave data coverage of the study region.We show that our new model is generally consistent with previous models but with improved resolution, particularly at greater depths and in the marginal sea areas. Thus, our model can constraints the structure of the crust and upper mantle better, helping reveal the important lithospheric features of East Asia. Specifically, here we focus on variations in lithospheric thickness, the collision and interaction between the Indian and Eurasian Plates, and the delamination mechanism of Eastern China.

2. Data

We used waveform records from more than 1,000 seismic stations in East Asia to derive high-quality Rayleigh-wave phases and group velocities based on the earthquake-based TS analysis method. Most of the stations are included within the dense China Regional Seismic Networks (CRSNs) (Zheng XF et al., 2010); the other stations included temporary deployments in the TP, the Tianshan region, and some other permanent stations within the Full Range Seismographic Network (F-NET) of Japan,South Korea, and Taiwan Islands, which are available from the IRIS DMC (Figure 1). Earthquake records from the CRSNs stations and the permanent stations in Japan and Taiwan covered the period from 2008 to 2011, which provided the largest contribution to our waveform dataset.The other stations were mostly temporary stations covering 2001 to 2010. The network codes and operation periods of the temporary stations are presented in Table S1.

Table S1. List of network codes and corresponding operation period.

In addition to these newly measured dispersions, we also combined the following three datasets: (1) data from Bao XW et al. (2015) that across the Chinese continent,including Rayleigh-wave group and phase dispersion measurements derived from ambient noise correlations and earthquake group velocity measurements; (2) Rayleighwave group and phase dispersion measurements from Wang Q et al. (2017) for the marginal seas in East Asia derived from ambient noise correlations; and (3) Rayleighwave phase velocities between 130 and 250 s from Ekström (2011) during the inversions for the S-wave velocity structure. The addition of longer-period data helped constrain the S-wave velocity structure down to a depth of approximately 300 km, enabling us to investigate variations in lithospheric thickness beneath East Asia.

3. Methods

3.1. Earthquake-based TS analysis

We used an earthquake-based TS analysis method(Yao HJ et al., 2006; 2008; Chen HP et al., 2014; Ojo et al., 2018) to measure the high-quality phase and group velocity between selected station pairs. Earthquake-based TS analysis requires that the target event and the two selected stations should be on a great circle; however, this is often challenging but ensured that the great circle paths between the earthquake and the two stations only have small offsets. We used the following two parameters to measure the offset: (1) the azimuthal difference, α, of the earthquake to the two stations, and (2) the azimuthal difference, β, between the earthquake to the nearest station and the nearest station to the other station (Yao HJ et al.,2006; Figure S1). The maximum permitted values were set to 2° and 3° for α and β, respectively. For a given earthquake, a station pair that satisfied these criteria was selected, and the cross-correlation function (CCF) of the seismograms recorded by the two stations were calculated.The CCF can be considered as an approximation of the inter-station Green’s function (Landisman et al., 1969),which is used to extract the inter-station phase and group velocity. To ensure good excitation of the surface waves,only events with magnitudes ＞5.0 and depths ＜100 km were included. For each station pair, the distance between the event and either station was limited to 30°-120° to ensure good surface-wave waveforms. The distribution of selected events is shown in Figure 1.

3.2. Dispersion measurement and quality control

The dispersions were measured using the frequencytime analysis method (Ritzwoller and Levshin, 1998),which was integrated into the automatic frequency-time analysis (AFTAN) program. AFTAN is mainly used to analyze CCFs from ambient noise but can also be adapted to measure dispersion from event CCFs by setting the phase shift to pi/4 instead of -pi/4 (Chen HP et al., 2014).Here, we used AFTAN to measure Rayleigh-wave phasevelocity dispersions at periods of 10-120 s and group velocity dispersions at periods of 10-140 s. We then checked every dispersion curve by comparing it with the theoretical velocity derived from the GDM52 model(Ekström, 2011) to identify and reject poor dispersion measurements. Only the smooth part of the dispersion curve was retained. We also rejected dispersion curves that have less than 30 periods (Wang WL et al., 2014). Finally,for each selected station pair, all the dispersion curves derived from multiple earthquakes were averaged to provide the final dispersion measurement.

3.3. Strategy for combining dispersion measurements

Bao XW et al. (2015) measured phase dispersions from ambient noise correlations and group velocities from earthquakes across the Chinese continent. Wang Q et al.(2017) also measured the phase and group velocities from ambient noise correlations across marginal seas in East Asia. We combined our new dispersion measurements with these two dispersion datasets to construct Rayleighwave velocity maps. A comparison between the surfacewave velocities from the TS and noise interferometry at each overlapping period indicates a systematic bias, with the values from the TS generally greater than those from the noise interferometry (up to approximately 0.5%;Figures S2 and S3). This bias has also been observed in previous studies (e.g., Yao HJ et al., 2008; Ojo et al.,2018). Therefore, we followed the procedure of Yao HJ et al. (2008) and Li MK et al. (2018) to combine the dispersions from TS and noise interferometry. Briefly, for periods shorter than 20 s, we used the dispersions obtained from noise interferometry; for periods longer than 50 s, we used the dispersions from the TS analysis; and for periods between 20 and 50 s, we used the dispersions from the noise interferometry if the bias was ＜0.1 km/s, otherwise,the dispersions from the TS analysis were used. The final dispersion dataset had a relatively denser path coverage and a wider period range. The final dispersion measurements were then used to construct phase and group velocity maps for a set of periods on 0.5°×0.5° grids using a linearized 2D inversion method (Bao XW et al., 2013;Xu Z et al., 2013). We used phase and group velocity images derived from the GDM52 model (Ekström, 2011)as the reference model to stabilize the inversion.Furthermore, for a certain period, by comparing the initial mean square travel time residual and the remaining uncounted residual, we rejected the paths with residual exceeded the unaccounted residual by three times or more.Subsequently, the velocity images were reinverted. The phase and group velocity maps produced after three iterations of this process were taken as the final velocity images.

3.4. Inversion for surface velocity map and S-wave velocity structure

The selected dispersion measurements were used to construct group and phase velocity maps at different periods using the linearized 2D inversion method of Xu Z et al. (2013), which is performed with a linearized inversion by minimizing the following quantity:

whereAis the coefficient matrix;mis the parameter vector; andλandφare the damping and smoothing parameters, respectively. The smoothing operator L is the Laplacian operator. We experimented with different pairs ofλandφ, and selected those that best balanced the data misfit, model damping, and smoothness. Thus,λ= 10 andφ= 300 were used in this study. Resolution tests (Figure S4) showed that a 1° × 1° grid could be recovered in most parts of the study region at periods of 30, 60, and even 90 s;and a 2° × 2° grid could be recovered at periods ＞90 s.Figures S5 and S6 show several representative phase and group velocity maps, respectively. Because we will mainly focus on the 3D shear-velocity model, detailed descriptions of the dispersion maps are not provided in the main text. Nevertheless, after constructing the phase and group velocity dispersion maps, we then extracted 10-120 s phase and 10-140 s group dispersion curves from our model as well as the 130-250 s phase dispersion curve from the GDM52 model (Ekström, 2011) at each grid point. We aimed to construct the shear-wave velocity structure of the crust and upper mantle (down to a depth of～300 km) to reveal the key features of the model. The long-period phase velocity dispersion proved much more sensitive to deeper structures than the group velocity.Therefore, we only used the long-period phase velocity dispersions from Ekström (2011). We merged the dispersion curves and inverted for 1-D shear-wave velocity profile at each grid point.

We jointly inverted the phase and group dispersion curves for each grid using a linearized iterative method(Herrmann and Ammon, 2004). The initial model was modified from the global 1D reference model, AK135(Kennett et al., 1995). Similar to Bao XW et al. (2015), we set the S-wave velocity in the top 100 km to a constant 4.5 km/s to avoid the assumption of a prefixed Moho interface.The deeper part of the initial model was obtained using the AK135 model. The layer thicknesses were 2 km in the top 10 km, 5 km at depths from 10 to 80 km, and 10 km below a depth of 80 km. Because surface wave dispersion curves are primarily sensitive to changes in S-wave velocity(although they also have some sensitivity to P-wave velocity and density; Julià et al., 2000; Bao XW et al.,2015), the S-wave velocity in the top 100 km was linked to the P-wave velocity using a Poisson’s ratio of 0.25 and scaled to obtain density using the empirical relationship of Brocher (2005). Finally, we obtained a 3D crustal and upper-mantle S-wave velocity model for East Asia to a depth of approximately 300 km. The inversion results and data fitting at two selected grid-point pairs are shown in Figure S7. Overall, the data from our dispersion maps and the GDM52 model matched the predicted dispersion matched well.

4. Results

We inverted the newly constructed surface-wave dispersion curve for each grid for the 1D S-wave velocity profile and then constructed a 3D S-wave velocity model beneath East Asia by assembling all the 1D velocity profiles. Based on the 3D model, we also constructed crustal and lithospheric thickness models. In this section,we first briefly describe the crustal thickness model(Section 4.1) and then focus on the 3D S-wave velocity model (Section 4.2) and the lithospheric thickness model(Section 4.3).

4.1. Crustal thickness

Dispersion data can effectively constrain the absolute S-wave velocities at different depths, even though they cannot constrain the depths of the velocity discontinuities(Julià et al., 2000). Xu Z et al. (2013) found that the average shear-wave velocity at the depth of the Moho as defined by the reference model CRUST 2.0 (Bassin et al.,2000) is almost constant at approximately 4.0 km/s, with a slight positive trend as a function of crustal thickness.Based on this, the depth of 4.0 km/s (the value can vary slightly) is a good estimate for the depth of the Moho.Although the method is simplistic, it has been successfully applied in previous studies to provide approximate estimations of crustal thickness (e.g., Bao et al., 2011,2013, 2015; Wang Q et al., 2017).

Figure 2 shows the estimated crustal thickness beneath East Asia based on our models. There are significant lateral variations in thickness ranging from an extremely thin crust (＜30 km) in eastern China and the marginal seas to intermediate thicknesses (～35-50 km) in Central China and the thickest crust (～70-80 km) in the TP region. The thickest crust beneath the TP reflects significant crustal thickening caused by the Indo-Asian collision, while the thin crustal thickness in eastern China and the marginal seas reflects the strong influence of the Mesozoic extension tectonics caused by dramatic lithospheric thermo-tectonic reactivation (Ren JY et al., 2002). Several crustal thickness maps of the Chinese continent from historical deep seismic sounding profiles (e.g., Teng JW et al., 2013), teleseismic receiver functions (e.g., Li YH et al.,2014), and surface wave tomography analysis (e.g., Bao XW et al., 2015) have recently been provided. Our new map shows not only good consistency with these previously published models but also covers the marginal seas, providing a broader-scale map of crustal thickness variations in East Asia. A comparison with the crustal thickness map of the Chinese mainland (Figure S8) shows good consistency between the two models. The largest discrepancy lies in the TP, where the distribution of seismic stations is extremely uneven. Importantly, our new map, based on dispersion inversion, provides a more uniform data coverage than that obtained from receiver function analysis.

Figure 2. Map of crustal thickness of East Asia derived from our 3D shear-wave velocity model.

Our new crustal thickness map also corresponds well with surface geology. For example, the North-South Gravity Lineament (NSGL) is well delineated by a sharp crustal thickness gradient that is compatible with the surface topography gradient. The coastline is also well identified by the crustal thickness transition zone in the sea-land interaction area. The Qilian orogeny in the northeastern TP is clearly identified by increased crustal thickness relative to its surrounding region. The margins of the TP are also well delineated by significant changes in crustal thickness associated with large elevation changes.

4.2. 3D S-wave velocity model of East Asia

We constructed a lithospheric S-wave velocity model from the Earth’s surface to a depth of approximately 300 km, with some selected depth slices shown in Figures 3 and 4. As shown in Figure 3a, the velocity generally correlates well with surface geology and regional tectonics at a depth of 5 km. Large sedimentary basins in Chinese mainland (e.g., the Tarim Basin, Qaidam Basin, Sichuan Basin, North China Basin, and Songliao Basin), the eastern marginal seas of East Asia, and the Japan Sea are dominated by low velocities. The Ordos Block shows relatively complex features in the western and central areas, and the Weihe Rift in the southeastern corner is dominated by low velocities. In comparison, the eastern part of this region is dominated by relatively high velocities suggesting stronger erosion in the eastern part of the Ordos Block. These features are consistent with those of previous studies on continental (Zheng SH et al., 2008;Sun XL et al., 2010; Xin HL et al., 2019; Xu Z et al.,20133; Bao XW et al., 2015; Shen WS et al., 2016) and regional scales (Bao XW et al., 2013; Guo Z and Chen YJ,2017; Wei XZ et al., 2017). Low-velocity zones can also be identified in the northwestern part of the TP and the eastern Himalayan Syntaxis (EHS) region. In contrast,some mountain regions, such as the east Qinling orogen,Changbai Mountain Range, and Da Hinggan Ling, show high-velocity features.

At a depth of 30 km (Figure 3b), the Japan Sea and the East and South China Seas reach the lower crust or uppermost mantle and exhibit high velocities. Western China, except for the Sichuan Basin, is dominated by low velocities, especially in the TP region and the northern part of Yunnan. The Sichuan Basin has a relatively higher velocity than the surrounding mountains. Northeastern China, the eastern part of north China, and southeastern China show intermediate velocities. In addition to the NSGL, there is a clear velocity contrast that correlates well with the sharp crustal thickness and surface topography gradient. In addition, the S-wave velocities in the Tarim Basin exhibit clear heterogeneity, with a small highvelocity zone in the southwestern part.

At a depth of 60 km (Figure 3c), most of the region reaches the uppermost mantle, except for the TP and the Pamirs and Tianshan regions, which remain at the depth of the lower crust. These regions are dominated by conspicuous low velocities. Some other regions, such as the surrounding regions of the Japan Sea and Ryukyu Island, are also clearly imaged as low-velocity zones. A prominent low-velocity anomaly is observed in the Datong volcanic region, which has also been identified in previous regional-scale surface and P-wave tomographic studies(Zhao DP et al., 2011; Chen HP et al., 2014), which may imply the upwelling of melt materials from the deeper mantle.

At depths of 90 and 100 km (Figures 3d and 4a), the Swave velocity images show predominantly low velocities beneath most parts of the TP (especially the northern part),the Songpan-Garzê terrane, eastern China, the Korean Peninsula, and the marginal seas (including the Japan Sea).These low-velocity anomalies are broadly consistent with global- and regional-scale tomographic studies (Shapiro and Ritzwoller, 2002; Friederich, 2003; Huang ZC et al.,2014; Huang JL and Zhao DP, 2006; Priestley et al., 2006;Feng M and An MJ, 2010; Wei W et al., 2012; Li YH et al., 2013; Bao XW et al., 2015; Li HY et al., 2018) but with higher resolutions. The prominent low velocities beneath the TP indicate significant deformation and high temperatures. Other predominant features at these depths are the high-velocity anomalies beneath the western part of the Yangtze Craton (including the Sichuan Basin), the Ordos Block, and the central and eastern parts of the Tarim Basin. The heterogeneities of the S-wave velocities inside these stable platforms imply different deformation states in these regions. High-velocity anomalies are also observed in the Taiwan Strait and the northern part of Japan Island.

Figure 3. Shear-wave velocity maps at 5, 30, 60 and 90-km depths from this study. White lines in (a) represent the locations of the vertical cross-sections in Figures 5—7. The black lines represent the boundaries of geological units as shown in Figure 1.

At a depth of 150 km, the high-velocity anomalies beneath the Tarim Basin, the western part of the Yangtze Craton, and the Ordos Block extend to this depth,suggesting the existence of thick cratonic roots beneath these regions. The low-velocity anomalies in the TP are concentrated to the northwestern part, while most parts of the TP are dominated by high-velocity anomalies, such as the Songpan-Garzê Terrane in the east and other southwestern areas. A clear banded zone with a relatively low velocity is observed in the southern TP region,extending from the Indian Plate to the central TP with a NNE strike direction and extending to a depth of 200 km(Figure 4c) and even 250 km in some places (Figure 4d).This banded low-velocity zone correlates with the wellknown Yadong-Gulu suture. The central and northern parts of Japan Island show a high-velocity anomaly that extends down to a depth of 200 km.

At depths of 200 km and below (Figures 4c-4d), the high-velocity anomalies in the previously mentioned cratonic regions disappear, including the Sichuan Basin,Ordos Block, and Tarim Basin. The western part of the Tarim Basin shows a small-scale high-velocity anomaly from a depth of 150 to 200 km, which may imply the deep subduction of the Indian Plate to the cratonic root in the western Tarim Basin. It is worth noting that the lowvelocity anomalies are concentrated in several banded zones located in the western and central parts of the TP,striking in a NNE direction from the Indian Plate to the northern and northwestern TP. These features may be related to slab tears of the subduction plates, as suggested by Li JT and Song XD (2018) based on P- and S-wave tomography. In general, eastern China and its marginal seas are dominated by widely distributed low velocities below a depth of 90 km (Figure 3d and Figure 4), which can be attributed to the intensive destruction of the lithosphere and the related upwelling asthenosphere.

Figure 4. Shear-wave velocity maps at 100, 150, 200 and 250-km depths from this study. The black lines represent the boundaries of geological units as shown in Figure 1.

Figure 5. Four west-east trending cross-sections across of the shear-wave velocity model from this study. The locations of the profiles are denoted in Figure 3a. The white lines and white crosses denote the estimated Moho interfaces and LABs,respectively.

Several representative cross-sections from the surface to a depth of approximately 300 km are presented in Figures 5-7. The locations of these cross-sections are shown in Figure 3a. These cross-sections have a number of common features. For example, the sedimentary basins and oceanic regions correlate with low-velocity anomalies near the surface of all the profiles, and the boundaries between different geological units are well imaged at the lithospheric mantle depth, i.e., the boundaries between the sedimentary basins and surrounding mountainous regions are well imaged as relatively low-velocity anomalies.

Figure 5 shows four west-east-trending cross-sections.The first two sections (A-A1andB-B1in Figures 5a and 5b,respectively) show a heterogeneous lithospheric mantle structure from west to east. Prominent high-velocity anomalies in the lithospheric mantle and a relatively clear lithosphere-asthenosphere boundary (LAB) are also observed beneath the sedimentary basins, especially those in central and western China, i.e., the Tarim Basin,Qaidam Basin, Ordos Block, and Sichuan Basin. Although the lithospheric mantle of the sedimentary basins in western and eastern China are all imaged as high-velocity anomalies, their velocity values show significant variation.The velocities of the basins in western China are higher than those in eastern China, suggesting changes in lithospheric strength from west to east. The second two sections (C-C1andD-D1in Figures 5c and 5d,respectively) cross the eastern part of the North China Craton (NCC) and northeastern China, respectively, and end in the oceanic regions. Subduction slab of the Pacific Plate to the Eurasian Plate is clearly seen as well as large low-velocity region beneath the Sea of Japan. Another important feature that can be observed from the four sections is the gradual thinning of the lithospheric thickness from west to east, which is consistent with the findings of Feng M and An MJ (2010).

Figure 6 shows four south-north-trending cross-sections.SectionsE-E1andF-F1(Figures 6a and 6b) cross the western and eastern TP, respectively, while sectionsG-G1andH-H1(Figures 6c and 6d) cross central and eastern China, respectively. SectionsE-E1andF-F1show that the underthrusting high-velocity Indian Plate approaches the central-to-northern part of the TP, at 34°N in the western TP and 35°N in the eastern TP. The subduction of the Indian Plate under the TP is associated with severe tectonic destruction. For example, the intrusion of weak materials from the TP to the lithosphere of the Tarim Basin under the extrusion force of the subduction slab also causes the destruction of the lithospheric base in the basin-mountain junction zone (see SectionE-E1, Figure 6a). In sectionFF1(Figure 6b), the lithospheric root of the Qaidam Basin is clearly imaged to have a high velocity; however, there is no indication of southward continental subduction of the Asian Plate under the lithosphere of the TP, which is consistent with the results of seismic tomography and receiver function analyses (e.g., Liang XF et al., 2012;Shen XZ et al., 2017; Wei XZ et al., 2017).

Figure 6. Same as Figure 5 but for four south-north trending cross-sections.

Figure 7. Same as Figure 5 but for a northwest-southeast (a) and a southwest-northeast (b) trending cross-sections.

Two other cross sections,I-I1andJ-J1, are presented in Figure 7. Section ProfileI-I1(Figure 7a) starts from the northern boundary of the Tarim Basin and crosses the Tarim Basin, Qaidam Basin, TP, Sichuan Basin, Yangtze Craton, South China Block, and South China Sea. SectionJ-J1(Figure 7b) starts from the southeastern corner of the TP and crosses the eastern TP, Qinling Orogen, and Ordos Block. Similar to sectionsA-A1andB-B1(Figure 5), these profiles show complex lithospheric mantle structures and significant differences between geological units, i.e., the velocity contrast between the basins and their surrounding regions. The lithospheric mantle is imaged as a highvelocity zone in the basins of these two sections, while the lithospheric roots are also easily recognizable.

4.3. Lithospheric thickness in East Asia

The lithosphere and asthenosphere are the most tectonically active parts of Earth’s interior. A traditionally accepted broad definition of the asthenosphere based on seismic tomography is the low-velocity zone in the upper mantle. In many seismic studies, the high-velocity lid over the asthenosphere has often been equated to the lithosphere; however, the interface between the lithosphere and asthenosphere is often far from a distinct discontinuity, and the transition from the lithosphere to the asthenosphere is typically identified as the lithosphere-asthenosphere boundary (LAB). The distribution of the LAB (or lithospheric thickness) is important for understanding the deformation and dynamic processes in the Earth’s deep interior. In this section, we focus on the implications of seismological studies for modeling lithospheric thickness, which we refer to as the seismic lithosphere, distinguishing it from the thermal lithosphere as inferred from temperature (e.g., An MJ and Shi YL, 2006).

The seismic lithospheric thickness (“lithospheric thickness” hereinafter) is more difficult to define using seismic velocity models (e.g., Li YH et al., 2013; Wang Q et al., 2017). In this study, we attempted to identify a proxy for the LAB by using the following strategies. Considering the relatively low velocity in the asthenosphere, we first discarded the 1D velocity profiles in which the velocity continuously increased with depth below the Moho interface. Several examples of discarded profiles are shown in Figure S9. Second, for each remaining 1D velocity profile, we determined the maximum velocity,v1,from the Moho to a depth of 200 km below, with the corresponding depth denoted as ‘H1’. We then identified the minimum velocity,v2, fromH1to a depth of 150 km below, and recorded the corresponding depth as ‘H2’. The velocity difference (%),R, was then calculated asR= (v1-v2) /v1×100%, and the 1D velocity profiles withRvalues＜2% were discarded. Third, for each remaining velocity profile, we calculated the velocity gradient betweenH1andH2, and took the depth corresponding to the maximum negative gradient as a proxy for the LAB. Examples clearly showing the LAB are provided in Figure S10. For comparison, we also show several 1D profiles with and without the LAB proxy in Figures S11 and S12,respectively.

By following these three steps, we obtained the lithospheric thickness map of East Asia from our S-wave velocity model (Figure 8). Although there are limitations in this approach, such as the measuring artifacts near the basin-mountain junction zones and subduction zones as well as the relatively arbitrary value ofR, the LAB generally indicates the location of a rapid velocity decrease with depth in most parts of our study region (denoted by the white crosses in the cross-sections in Figures 5-7). As such, our approach is useful for examining the general variations in lithospheric thickness at a continental scale.

The lithospheric thickness map of East Asia derived from our velocity model is shown in Figure 8. We compare our map with the results of An MJ and Shi YL (2006), who first converted the 3D shear-wave velocity structure into a temperature structure and then extracted lithospheric thickness. Because of differences in our definitions of the LAB, the values inevitably differ slightly; however, the two maps are generally consistent in absolute values,perturbation patterns, and local features. For example, in both cases, lithospheric thinning is seen in eastern China,and a thick lithosphere is seen in the southern TP region and the Sichuan Basin (Figure S13).

Based on these maps, lithospheric thickness shows significant lateral variation, ranging from extremely thin in the Japan Sea to its thickest in the western part of the Tarim Basin and the southeastern part of the Sichuan Basin. Except for the Yunnan region, western China (west to 110°E) has a thick lithosphere, while eastern China and its marginal seas (east to 110°E) show a relatively thinner lithosphere. Generally, the lithosphere thins from west to east. Eastern and southwestern China have consistent lithospheric thicknesses, while central and western China show large variations. Notably, the boundary between eastern and southwestern China and western China is marked by a sharp change in lithospheric thickness. As shown in Figure 8, the western Yangtze Craton, the majority of the central and eastern TP, and the western Tarim Basin have a thick lithosphere (thicker than 180 km). Some other smaller regions, such as the northwestern and southwestern corners of the Ordos Block and the northern Qilian Orogen, also show a thick lithosphere. The eastern and southeastern Tarim Basin, Qaidam Basin,Qilian Orogen, east Qinling-Dabie Orogen, central and eastern Ordos Block, and the conjunction zone between the Sichuan Basin and TP show a relatively thin lithosphere.

Figure 8. Map of lithosphere thickness derived from our 3D shear-wave velocity model. The white lines represent the boundaries of geological units as shown in Figure 1.

In northern and eastern Tibet, the velocities increase continuously with depth below the Moho interface (Figure S12); thus, we did not determine the LAB as previously described and the lithosphere in these regions may be thicker than the estimated thicknesses. The complex tectonic background in central and western China, including plate subduction and basin-mountain interactions,may contribute to the significant variations in lithospheric thickness. Crucially, despite some discrepancies in LAB values caused by different identification criteria and geophysical techniques, the general trend of our lithospheric thickness map is consistent with that presented in existing continental-scale (An MJ and Shi YL, 2006; Li YH et al., 2013) and regional-scale studies (Zhang XM et al., 2014; Wang Q et al., 2017).

5. Discussion

We generated a 3D lithospheric S-wave velocity model for East Asia via the joint inversion of newly measured Rayleigh wave dispersion data from earthquake-based TS analysis, dispersion datasets from Bao XW et al. (2015) for Chinese mainland and Wang Q et al. (2017) for the marginal seas in East Asia, and the long-period dispersion map from Ekström (2011). This combined dataset broadened the period range of Rayleigh wave dispersion(periods from 10 to 140 s for group velocity and 10 to 250 s for phase velocity), allowing us to investigate the S-wave velocity structure in detail to a depth of approximately 300 km. Our new model has both high lateral and depth resolution; thus, we are able to reveal the key features of lithospheric thickness in East Asia, which are important for studying the geological evolution and tectonics of the region. In this section, after a brief comparison of our new model with some previously published continental-scale models, we focus on two important topics: (1) the lithospheric structure of eastern China and its marginal seas, and (2) subduction and collision interactions in the TP region.

5.1. Comparison with previous studies

Some continental-scale surface-wave tomographic studies of similar regions (e.g., Feng M and An MJ, 2010;Li YH et al., 2013; Pandey et al., 2014; Xin HL et al.,2019), and regional-scale tomographic studies for parts of our study region including Chinese mainland (e.g., Zheng SH et al., 2008; Sun XL et al., 2010; Xu Z et al., 2013;Bao XW et al., 2015; Shen WS et al., 2016), the TP and surrounding regions (e.g., Bao XW et al., 2013; Zhang XM et al., 2014; Wei XZ et al., 2017), and the region’s marginal seas (e.g., Wang Q et al., 2017) have been performed. Here, we compare our model with some selected continental-scale tomographic models (e.g., Feng M and An MJ, 2010; Li YH et al., 2013; Pandey et al.,2014; Bao XW et al., 2015). The S-wave velocity models presented by Feng M and An MJ (2010), Li YH et al.(2013), and Pandey et al. (2014) cover a similar region and have a similar depth range to that of our new model;however, these studies used a few hundred stations and the lateral resolutions of their models are comparatively limited. Therefore, our model not only shows good agreement with these previous studies with respect to large-scale structures but also significantly improves the data resolution for some small-scale structures. For example, our model more clearly reveals lithospheric structures at basin-mountain junction zones as well as lithospheric roots beneath sedimentary basins, such as the Qaidam Basin and Ordos Block.

This work builds on two earlier studies from our research group on Chinese mainland (Bao XW et al., 2015)and the marginal seas in East Asia (Wang Q et al., 2017).By combining our newly measured dispersion data and the dispersion maps from Ekström (2011) with our group’s previous two datasets, we have been able to (1) further increase the path coverage, and (2) extend the effective inversion depth to approximately 300 km. In comparison with the two models from Bao XW et al. (2015) and Wang Q et al. (2017), the effective depth of our new model reaches the asthenosphere, and even deeper parts, which allows us to investigate the features of the entire lithosphere in East Asia, including variations in lithospheric thickness.

5.2. Lithospheric thinning and destruction in east China and marginal seas

East China contains the North China Craton (NCC),South China Block (SCB), and Phanerozoic orogenic belts between the NCC and SCB (Zhao GC and Cawood, 2012).Figure 8 shows that eastern China and the marginal seas have a coherent lithospheric thickness, averaging ～80 km,with a clear thinning trend from west to east (see crosssectionsA-A1andB-B1in Figures 5a and 5b). For the NCC, the lithosphere is still thick in its western part, such as the eastern part of the Ordos Block, and thins gradually to the eastern part. Similarly, for the SCB, the lithosphere is quite thick in the Yangtze Craton (YZC) and becomes significantly thin in its eastern part and marginal seas(Figure 8). Thus, our results show that both the eastern NCC and SCB and their conjunction regions are characterized by lithospheric thinning. This is consistent with global and regional seismic tomographic studies(Feng M and An MJ, 2010; Bao XW et al., 2013, 2015; Li YH et al., 2013; Pandey et al., 2014; Wang Q et al., 2017)and receiver function studies (Chen L, 2009), which also report thinning and destruction of the lithosphere in east China and its marginal seas (e.g., Chen L et al., 2009; Zhu RX et al., 2012; Li YH et al., 2013; Pandey et al., 2014;Bao XW et al., 2015). Nevertheless, the mechanism of lithospheric keel destruction remains unclear (Zhu RX et al., 2012). Previous studies have proposed several tectonic processes that may be responsible for the thinning of the lithosphere in eastern China and its marginal seas including the far-field effect of the collision between the Indian and Eurasian Plates, mantle plume activity, and the westward subduction of the Pacific Plate under the Eurasian Plate [see Zhu RX et al. (2012) for a review].Body-wave tomographic studies (e.g., Zhao L et al., 2012)have confirmed the existence of a deeply subducted slab from the Pacific Plate and its horizontal extension beneath east China, especially beneath Chinese mainland. The initial part of the subduction slab in the upper mantle,previously imaged as a high-velocity zone, is also observed in our model (see the eastern sides of crosssectionsC-C1andD-D1in Figures 5c and 5d). The lowvelocity anomalies observed at depths of ～80-180 km depth beneath eastern China and its marginal seas are particularly notable (Figures 5 and 6d). Similar lowvelocity zones have been observed in recent P-wave tomographic studies (e.g., Jiang GM et al., 2013; Zheng HW and Li TD, 2013) and surface-wave tomographic studies (e.g., Ouyang LB et al., 2014) carried out in the middle-lower Yangtze River based on the joint inversion of the receiver function and surface-wave dispersion in the SCB (Li HY et al., 2018). These low-velocity zones are interpreted as upwelling of the asthenosphere induced by the detached lithosphere, which could be related to the subduction of the Pacific Plate (Jiang GM et al., 2013;Ouyang LB et al., 2014). As shown in Figure 4, the velocities in these large low-velocity zones are significantly lower than those in the surrounding regions in the asthenosphere, and are possibly caused by the partial melting of the subduction slab or the detached lithosphere induced by the subduction slab.

Figure 6d shows the south-north trending velocity cross-section across east China, with distinct low-velocity zones in the asthenosphere that suggest the widespread existence of partial melting. Indeed, the upwelling of the melting materials likely causes a reducing velocity contrast between the high-velocity lithospheric mantle lid and the asthenosphere, and even the complete disappearance of the LAB, especially in oceanic regions (see the eastern part of cross-sectionsC-C1andD-D1in Figure 5). Complete intrusions of partial melting into the lithosphere are also observed in some continental regions, i.e., 114°E in sectionB-B1(Figure 5b) showing the upwelling of melting materials near the NSGL. The upwelling of asthenospheric materials induced by the westward subduction of the Pacific Plate can erode the overlying lithospheric base,resulting in thinning of the lithosphere. Previous work has indicated that these regions have been strongly influenced by the tectonics of the circum-Pacific system (Ren JY et al., 2002). Therefore, the subduction of the Pacific Plate under the Eurasian Plate might be the main cause of the destruction of the lithosphere in eastern China and its marginal seas (Zhu RX et al., 2012).

5.3. Subduction and collision interactions in the TP region

The ongoing collision beginning ～50 Ma ago between the Indian and Eurasian Plates created the TP as a highly deformed zone with thickened crust and elevated topography (Avouac and Tapponnier, 1993; Yin A and Harrison, 2000; Royden et al., 2008; Tapponnier et al.,2001). The complex tectonic evolution and deformation pattern here make it an ideal place to investigate the process of continental collision. Indeed, the subduction of the Indian Plate beneath the Eurasian Plate has played a key role in the tectonic evolution of Tibet. To better understand the formation and tectonic evolution of the TP,numerous geophysical studies based on different techniques have been proposed to investigate the Tibetan velocity structures (Kind et al., 2002; Tilmann et al., 2003;Liang CT and Song XD, 2006; Li C et al., 2008; Feng M and An MJ, 2010; Hung et al., 2010; Liang XF et al., 2012;Li YH et al., 2013; Pandey et al., 2014; Bao XW et al.,2015; Chen M et al., 2017; Li JT and Song XD, 2018). In this section, we discuss the deep velocity structures beneath the TP as revealed by our model alongside their possible tectonic implications.

Our new model reveals a strongly heterogeneous Swave velocity structure at the lithospheric mantle depth beneath the TP (Figures 3 and 4). A distinctive feature is the presence of a large low-velocity zone in the northern and northeastern TP, ranging from ～90 to 150 km. The scale of this zone decreases with depth, its extension covers the Qiangtang terrane and the entire Songpan-Garzê terrane at 90 and 100 km depths, respectively, and it is concentrated in the northwestern part of the TP at a depth of 150 km. A low-velocity feature at the lithospheric mantle depth has also been observed in many previously published seismic studies (Friederich, 2003; Huang JL and Zhao DP 2006; Priestley et al., 2006; Li C et al., 2008; Li YH et al., 2013; Bao XW et al., 2015). The spatial extension of this low-velocity zone corresponds well with the widespread Cenozoic magmatism in the northern part of the TP (Chung et al., 2005), suggesting that the hot and weak materials in the lithospheric mantle may be related to mantle upwelling induced by lithospheric delamination due to collision between the Indian and Eurasian Plates(Chung et al., 2005; Li YH et al., 2013; Bao XW et al.,2015). Furthermore, the large low-velocity zone beneath the TP may suggest lithospheric delamination, as reported by some surface-wave tomographic studies (e.g., Li YH et al., 2013; Bao XW et al., 2015).

Another important issue in the TP region is the collision state of the Indian and Eurasian Plates, such as the northward extension of the subducted Indian lithosphere and the existence of south-verging subduction of the Asian Plate under the TP. Our model shows the subducting Indian slab in the southern part of the TP,which is imaged as having high velocities (Figures 6a and 6b) and is consistent with previous seismic studies(Tilmann et al., 2003; Li C et al., 2008; Feng M and An MJ, 2010; Hung et al., 2010; Li YH et al., 2013; Pandey et al., 2014; Bao XW et al., 2015). SectionsE-E1andF-F1in Figures 6a and 6b indicate that the underthrusting of the Indian slab approaches the middle to north TP regions.SectionE-E1is located in the western part of the TP and shows that the Indian slab approaches ～34°N and is then separated from the high-velocity Tarim cratonic keel by the low-velocity zone. SectionF-F1, on the other hand, is located in the eastern part of the TP region and shows that the Indian slab approaches ～34-35°N and is then separated from the high-velocity Qaidam cratonic keel by a narrow low-velocity zone. Ceylan et al. (2012) observed highvelocity zones within Tibet and inferred that the northern limit of the Indian slab is up to ～34°N. The research region examined by Ceylan et al. (2012) contains ourF-F1section, confirming the consistency of our results.

Figure S1. Illustration of two deviation angles α and β defined in the earthquake-based two-station (TS) method (Yao HJ et al., 2006). The parameters α and β are defined as the azimuthal difference of the earthquake to the two stations, and the azimuthal difference between the earthquake to the nearest station and the nearest station to the other station, respectively.

Figure S2. Histogram of the Rayleigh wave phase velocity difference between the measurements from the earthquake-based TS method and the ANT at the overlapping 15, 20, 25, and 30 s, with the means of 0.013, 0.020, 0.022, and 0.020 km/s,respectively. Letter "T" represents period.

Figure S3. Same as Figure S2 but for 35, 40, 50, and 60 s, with the means of 0.007, 0.007, 0.013, 0.014, respectively. Letter"T" represents period.

Figure S4. Checkerboard tests at periods 30, 60, 90 and 120 s.

Figure S6. Group velocity dispersion maps at periods 10, 30, 60, 80, 100 and 120 s.

Figure S7. Inversion results at grid points (115.0, 31.0) (a) and (100.0, 35.0) (b). The black and blue triangles represent phase dispersion that are extracted from our phase velocity dispersion maps and GDM52 model, respectively. The red triangles represent group dispersion that are extracted from our group velocity dispersion maps. The yellow and green triangles represent the predicted phase and group velocity dispersion, respectively.

Figure S8. Crustal thickness map of Chinese continent from this study (a) and Li YH et al. (2014) (b).

Figure S9. Examples of 1-D velocity profiles having no negative velocity gradient at the depth below Moho interface.

Figure S10. Examples of 1-D velocity profiles having clear LAB.

Figure S11. Typical examples for 1-D velocity profiles having clear identified LAB.

Figure S12. Typical examples for 1-D velocity profiles having no proxy for LAB.

Figure S13. Lithosphere thicknss from An and Shi (2006) (left panel) and this study (right panel).

Inconsistent with some previous geophysical studies(e.g., Zhu JS et al., 2002; Priestley et al., 2006), our work shows that the Indian Plate does not subduct beneath the entire TP, and its northern limits vary from east to west.However, our results are consistent with more recent tomographic studies (Li C et al., 2008; Hung et al., 2010;Zhang H et al., 2012; Li YH et al., 2013), which also support the partial underthrusting of the Indian slab under the TP and the variable subduction northern limits of the Indian slab from west to east. Using high-resolution P- and S-wave tomography, one recent study showed that the subducted Indian mantle lithosphere has variable subduction angles and northern limits in different regions(Li JT and Song XD, 2018).

6. Conclusions

We derived a 3D lithosphere S-wave velocity model for East Asia by combining newly obtained Rayleigh wave dispersions using the earthquake-based TS method, two previously published dispersion datasets (Bao XW et al.,2015; Wang Q et al., 2017), and the long-period dispersion maps derived by Ekström (2011). Our new model provides enhanced lateral and vertical resolutions compared with previous continental models, and provides a new valuable reference model for better understanding tectonics and geodynamic processes in East Asia.

Our new velocity model reveals heterogeneous lithospheric structures beneath East Asia. The enhanced effective depth (down to approximately 300 km) allows us to investigate the velocity structure features of the entire lithosphere, and we derive a lithospheric thickness map from the 3D velocity model. The lithospheric thickness map shows drastic variations from west to east. The boundary between eastern and southwestern China and western China is marked by a sharp change in lithospheric thickness, and East China and its marginal seas clearly indicate lithospheric thinning.

We observed widely distributed low-velocity zones in the asthenosphere above the sub-horizontal slab from the Pacific Plate beneath eastern China and its marginal seas,suggesting the ascent of hot partial melts derived from the subducted slab. Subduction of the Pacific Plate under the Eurasian Plate might be the main cause of the destruction of the lithosphere in eastern China and its marginal seas.We also observe a large low-velocity zone at a depth range of ～90-150 km beneath the northern and northwestern parts of the TP region indicative of lithospheric delamination. Furthermore, our model clearly images the subduction slab of the Indian Plate in the southern TP region. Here, high-velocity zones presumed to be the Indian slab are separated from the cratonic keels of the Tarim and Qaidam Basins by low-velocity zones, implying that the Indian Plate does not subduct beneath the entire TP. Finally, our results also provide evidence that the lateral extent of the Indian Slab varies from west to east.

Data availability

The 3D S-wave velocity, crustal thickness, and lithospheric thickness models as well as the phase and group velocity maps at different periods described in this study are available at https://www.equsci.org.cn/article/doi/10.10 16/j.eqs.2022.05.004?pageType=en in the Supplements.

Acknowledgements

This research was supported jointly by the National Natural Science Foundation of China (Nos. U1939204,41774056, and 41704046) and the Fundamental Research Funds for the Central Universities (No. 2042020kf0010).The seismic waveform data were provided by the Data Management Center of the China National Seismic Network at the Institute of Geophysics, China Earthquake Administration (Zheng XF et al., 2010) and the Data Management Center at the Incorporated Research Institution of Seismology (IRIS). Most of the figures were prepared using Generic Mapping Tools (GMT) (Wessel and Smith, 1998). We appreciate the constructive and detailed comments of the two anonymous reviewers.

Supplements