Shallow seismic structure beneath the continental China revealed by P-wave polarization, Rayleigh wave ellipticity and receiver function

SUMMARY

We construct a high-resolution shallow 3-D seismic model in the top 10 km of the upper crust in the continental China, with constraints of P polarization, Rayleigh wave ellipticity and receiver function obtained from records of 3848 seismic stations. Our 3-D seismic model has a spatial resolution of 0.6–1.2° in the north–south seismic belt and the trans-north China orogen, and 1–2° in the rest of the continental China (except the Tarim basin and the southwest Tibet). The seismic model exhibits low velocity anomalies of deposits in major sedimentary basins and high velocity anomalies of crustal bedrocks in young orogenic belts and old tectonic blocks. The inferred sediment thickness maps display thick deposits in major sedimentary basins, some compacted sediments in the intermontane basins in young orogenic belts and little sediments in old tectonic blocks. We also discuss compaction effects of the sediments and implications of tectonic history and geological evolution of the major basins in the continental China based on the inferred seismic models. This study provides an effective mean of seismic imaging through joint inversion of various seismic constraints and establishes a framework of seismic data sharing for future studies in the seismological community in a first step of developing a China Seismological Reference Model.

Composition and structure of the continental crust, Asia, Joint inversion, Crustal imaging, Seismic noise, Sedimentary basin processes

1 INTRODUCTION

The continental China is an assemblage of many tectonic blocks that were created in various geological time periods, including cratonic blocks of the North China craton and the Tarim craton in the Archean and the South China Block in the Precambrian, and microcontinental blocks of the Junggar, Qiangtang, Lhasa, Qilian, Alxa and Qaidam blocks formed in the Archean or Precambrian (Zhai 2013; Zhao et al. 2018, Fig. 1a and Table 1). In the first order, the assemblage was accomplished by the subduction of the Pacific Plate in the east, and the closure of the Tethys ocean and later continental collisions between the Indian and Eurasian plates in the west (Jia et al. 2013). Those tectonic forces also promoted strong interactions between the internal tectonic blocks in the continental China, generating several orogenic belts in the regions between the tectonic blocks, such as the Qinling range between the North China craton and the South China Block, and the Taihang range between the North China Plain and the Ordos craton (Meng & Zhang 2000a; Wang & Li 2008). Over the geological history, the complex tectonic activities and possibly some upwellings in the upper mantle produced many large sedimentary basins in the continental China, most prominently, the Songliao and Bohaiwan basins in the eastern China that are evolved from a series of rift and subsidence events related to the western subduction of the Pacific Plate (Ye et al. 1985; Ren et al. 2002; Feng et al. 2010), and the Sichuan, Ordos, Tarim and Qaidam basins in the Circum-Tibetan Plateau Basin-Range System that are influenced by the collision between the Indian and Eurasian landmasses (Ye et al. 1985; Xia et al. 2001; Ren et al. 2002; Yang et al. 2005; Wang et al. 2007a, 2016; Feng et al. 2010; Jia et al. 2013; Zhai 2013).

Figure 1.

(a) Geological units and major sedimentary basins in and around China and (b) seismic stations used in the study. Red-crossed regions in (a) are major basins (Deng et al. 2003; Shen et al. 2016), with their numerical labels identifying the basin names in Table 1. The regions enclosed by black lines are geological blocks, with their letter labels identifying the block names in Table 1. (b) Triangles show seismic stations, with blue ones from the China National Seismic Network (CNSN), red ones from the China Seismic Array (ChinaArray) and green ones from the Incorporate Research Institutions for Seismology (IRIS).

Table 1.

Geological units and major basins.

Tectonic blocks (letter label)	Major Basins (number label)
Junggar Block (A)	Junggar Basin (1)
Tarim Block (B)	Tarim Basin (2), Turpan Basin (3)
Qaidam Block (C)	Qaidam Basin (4)
Songpan-Ganzi Terrane (D)
Qiangtang Block (E)	Qiangtang-Tanggula Basin (5)
Lhasa Block (F)	Cuoqing-Lunpola Basin (6)
Qilian Block (G)
Alxa Block (H)	Ejina Basin (7), Jiuquan-Minle-Wuwei Basin (8)
Ordos Block (I)	Ordos Basin (9)
Xingan-Mongolia Block (J)	Erlian Basin (10), Hailaer Basin (11)
North China Plain (K)	Bohaiwan Basin (12), Taikang-Hefei Basin (13)
South China Block (L)	Sichuan Basin (14), Nanyang Basin (15), Jianghan Basin (16)
Chuandian Terrane (M)	Chuxiong Basin (17)
South Yunan Block (N)	Lanping-Simao Basin (18)
East Shandong/Yellow sea Block (O)	Beihuanghai Basin (19), Subei-Nanhuanghai Basin (20)
Greater Khingan Range (P)
Lesser Khingan Range (Q)	Songliao Basin (21)
	Sea of Japan backarc Basin (22), East China sea Basin (23), Pearl river mouth Basin (24), Beibuwan Basin (25)

Note: Letter and number labels are the same as those in Fig. 1(a).

Imaging the shallow seismic structure in the continental China would provide improved understanding of the geological history of the region and a seismic reference model for better assessing earthquake hazards and better constraining the seismic structure in the deeper part of the crust and in the mantle. High-resolution shallow structure could be used to place constraints on the composition of the bedrocks in each tectonic terrane, providing better characterization and identification of tectonic blocks. The inferred seismic properties of the sediments, sediment thickness and geometrical structure of a basin would help better infer geological evolution history of the basin (Liu 1998; Ren et al. 2002; Zhou et al. 2006; Wang et al. 2007b; Li et al. 2012). A more precise shallow seismic model would also allow more accurate estimation of ground motions during earthquakes and, with the seismic wave propagation more accurately quantified, better study of earthquake location and source mechanism. This is especially true in the regions of sedimentary basins where millions of people live in the metropolitans inside and where seismic ground motions are significantly amplified due to geometric focusing of the basin edge, and seismic reverberation and amplification of the sedimentary deposits (Field et al. 1997; Graves et al. 1998; Wald & Graves 1998; Davis et al. 2000). Finally, high-resolution shallow seismic structure would allow the seismic effects of shallow crustal structure more accurately accounted for in the seismic wave propagation, yielding better resolutions in seismic imaging of deeper structures.

Crustal seismic structures in the most areas of the continental China have been extensively studied with the active- or passive-source methods in the past several decades. Many active-source studies provided high-resolution images of crustal structures in many regions of the continental China, providing insights into many geological phenomena at the continental or regional scale (Kan et al. 1986; Auger et al. 2001; Liu et al. 2006; Qiu et al. 2007; Wang et al. 2007a; Guo et al. 2013; Teng et al. 2013; Chen et al. 2017; She et al. 2018; Tian et al. 2018). At the same time, passive source studies, which utilized the seismic data recorded from the earthquake sources, have also provided seismic images of crustal structures in many regions of the continental China, with the resolutions increasing drastically in recent years due to the improvements in spatial coverage and data quality of the permanent and temporary seismic stations in the region. Indeed, various approaches have been adopted in those passive-source imaging studies, including ambient noise tomography method (Yao et al. 2006; Li et al. 2009; Arroucau et al. 2010; Lai et al. 2014; Bao et al. 2015; Shen et al. 2016), Rayleigh wave ellipticity method (Scherbaum et al. 2003; Tanimoto et al. 2013), body wave amplitude ratio analysis (Bao & Niu 2017; Chong et al. 2018; Park & Ishii 2018; Wang et al. 2019) and joint inversion of the surface wave velocity dispersion, Rayleigh wave ellipticity or receiver function (Picozzi & Albarello 2007; Lin et al. 2012; Kang et al. 2016; Li et al. 2016b; Yang et al. 2019). Despite of the abundance of those seismic studies concentrating on the seismic structure in the whole crust, there have been few studies that focused on resolving seismic structure of the shallow crust in the continental China.

In this study, we focus on constructing a high-resolution seismic model in the top 10 km of the crust in the continental China, using constraints of P polarization, Rayleigh wave ellipticity and receiver function (RF) obtained from the seismic data recorded by 3848 seismic stations installed in or around China. Our goals are several folds: (1) to provide seismic results that can be used as constraints on the composition of the bedrocks and sediments in the continental China and as a reference of shallow seismic structure for seismic hazard assessment and future high-resolution studies of deeper seismic structure and (2) to gain insights into the tectonic and geological evolutions of major sedimentary basins in the region. The study is also aimed at developing a mean for jointly inverting various seismic constraints to obtain seismic structure and establishing a framework of seismic data sharing for future studies in the seismological community in a first step of developing a China Seismological Reference Model.

The shallow seismic structure is obtained through this approach: we first search the 1-D Vs structure beneath each station that best fits the P polarization, Rayleigh wave ellipticity and RF observed at the station; we then obtain a shallow 3-D seismic structure in the region through interpolation from the best-fitting 1-D seismic structures beneath the stations. We introduce seismic data in Section 2, present the joint inversion method in Section 3, describe shallow 3-D seismic structure, sediment thickness and their comparisons with published results in Section 4, and discuss geological implications of the inferred seismic results in Section 5.

2 DATA COVERAGE

We establish databases of continuous waveform and tele-seismic record to constrain the shallow seismic structure beneath the continental China. The continuous waveform database includes three-component waveforms of (1) 929 permanent seismic stations in the China National Seismic Network (CNSN) from 2015 to 2017, (2) 1193 temporary or permanent seismic stations installed in China and nearby regions from 1990 to 2017 with data available from the Incorporated Research Institutions for Seismology (IRIS) and (3) 1726 temporary seismic stations from the China Seismic Array (ChinaArray) from 2006 to 2016. Tele-seismic record database contains three-component waveforms of all 3848 seismic stations from the CNSN (2012–2017), the IRIS (1990–2017) and the ChinaArray (2006–2016) during the intermediate and large tele-seismic events (seismic magnitude M > 5.5 and epicentre distance from |${30^\circ }$| to |${95^\circ }$|⁠) occurring during the operating time of seismic stations (catalog from the IRIS and with a total of 9691 seismic events). Most stations (3520) are equipped with broad-band seismometers with a sampling rate of 100 Hz, and some of stations (328) from the IRIS and the CNSN belong to various types of seismometers with a sampling rate above 10 Hz (Fig. 1b).

3 JOINT INVERSION FOR 1-D SEISMIC MODELS BENEATH STATIONS

The joint inversion for the best-fitting 1-D seismic structure beneath each station consists of two steps: (1) we first infer the near-surface average Vs structure from polarization analysis of the first arriving tele-seismic P waves and (2) we then obtain the shallow seismic structure based on joint inversion of both Rayleigh wave ellipticity and P wave RF extracted from the seismic data recorded at the station, with the constraint of the near-surface Vs structure inferred in the P polarization analysis. Rayleigh wave ellipticity is chosen in the frequency range that is sensitive to Vs profiles at the shallow depth, while RF is retained to high frequency contents as it is sensitive to velocity contrasts across the internal seismic discontinuities. Specifically, Rayleigh wave ellipticity is represented by the ZH ratio, defined as the spectral amplitude ratio between the vertical and horizontal component waveforms of the Rayleigh wave dominated seismic noise, that is the recording that has a phase shift similar to Rayleigh wave (from 60° to 120°) between two components (Tanimoto et al. 2013).

3.1 Near-surface Vs structure from P polarization analysis

The body-wave polarizations carry significant information about near-surface velocity structures (Park & Ishii 2018). In this study, we obtain the near-surface shear-wave velocities (⁠|$V_s^{ns}$|⁠) beneath seismic stations based on polarization analysis of the direct P waves radiated by tele-seismic events.

We use three components of P waveforms recorded in the database for events occurring in an epicentre distance range from |$30^\circ $| to |$95^\circ $|⁠, with an event depth greater than 60 km and an event magnitude larger than 5.5. All recorded three-component waveforms are corrected with the instrumental response, rotated to the great circle plane and bandpass filtered in a frequency range of 0.04–2 Hz. The filtering scheme is chosen based on the best quality of P wave data observed in that frequency range. P wave arrivals are picked with an autoregressive Akaike information criterial (AR-AIC) function method (Chen & Holland 2016). Waveform segments from 1 s before and 4 s after the picked arrivals are used for polarization analysis. We adopt the principal component analysis algorithm (Jurkevics 1988) to measure the P wave apparent polarization angle (APA) and polarization degree, defined as the ratio of the primary eigenvalue over summation of two eigenvalues. We also measure the signal-to-noise ratio (SNR), defined as the ratio of P phase amplitude over the standard deviation (s.d.) of the ambient noise amplitude in a time window 10–5 s before the picked P arrival. Only the high-quality APA measurements with a SNR > 3 and a polarization degree |$\ge $|0.7 are retained for further analysis.

We invert the APAs measured from seismic records of a seismic station to obtain the near-surface shear velocity (⁠|$V_s^{ns}$|⁠) beneath the station. A best-fitting |$V_s^{ns}$| is searched through minimizing the following misfit function:

$$\begin{eqnarray*} f \left( m \right) = \frac{{\mathop \sum \nolimits_{i = 0}^N {{\left( {\bar{\theta }_i^{obs} - \bar{\theta }_i^{syn}\left( m \right)} \right)}^2} \times {\varphi _i}}}{{\mathop \sum \nolimits_{i = 0}^N {\varphi _i}}}, \end{eqnarray*}$$

(1)

where |${\boldsymbol{\ }}m$| represents the |$V_s^{ns}$| with a model space of 0.5–5 km s^–1, |$\bar{\theta }_i^{obs}$| and |$\bar{\theta }_i^{syn}( m )$| are the measured and predicted APAs for ith event, respectively, and |${\varphi _i}$| is the polarization degree of the ith measurement used as a weight for the inversion. |$\bar{\theta }_i^{syn}( m )$| is related to m by:

$$\begin{eqnarray*} \bar{\theta }_i^{syn} \left( m \right) = 2\arcsin \left( {m\times{p_i}} \right), \end{eqnarray*}$$

(2)

where |${p_i}$| is the P wave ray parameter of the ith event, determined based on the Preliminary Reference Earth Model (PREM, Dziewonski & Anderson 1981) (Fig. 2a).

Figure 2.

Near-surface Vs inversion with P wave polarization data for an example seismic station FJ.MXXF in the South China Block. (a) Data fitness of the observed P wave apparent polarization angles (APAs) in one of 5000 near-surface Vs inversions. Observed high-quality APAs (dots) are shown as a function of theoretical ray parameter calculated using PREM (Dziewonski & Anderson 1981), and coloured with the polarization degree or inversion weight, along with the synthetic APA curve (red line), predicted with the inverted near-surface Vs (3.42 km s^–1) that minimizes the misfit function. (b) Histogram of inverted near-surface Vs in the 5000 times of inversions based on the bootstrapping technique (blue histogram), with the mean (black dashed line) and s.d. (grey-shaded area) taken as the final near-surface Vs (3.30 km s^–1) and associated uncertainty (0.1 km s^–1) beneath the station, respectively.

We use a bootstrap technique (Efron 1992) to estimate the final best-fitting |$V_s^{ns}$| and its uncertainty beneath a station. Inversions are performed 5000 times on measurement subsets randomly selected from all the measurements, resulting in an assemblage of |$V_s^{ns}$| from 5000 inversions. The average of the assemblage is regarded as the final best-fitting |$V_s^{ns}$| and the associated s.d. is defined as its associated uncertainty (Fig. 2b).

This polarization analysis method is applied to the seismic data recorded at all stations (Fig. S1a), and a near-surface Vs structure beneath the continental China is obtained by interpolating the near-surface Vs inferred from P polarization analysis of 3705 seismic stations (Fig. 3a, see support material for details of the interpolation scheme). The near-surface Vs structure geographically correlates with the geological units. It exhibits low velocities in major sedimentary basins, such as the Bohaiwan, Sichuan, Nanyang, Subei-yellow sea, Songliao, Ordos, Qaidam, Junggar and Tarim basins, and high velocities in the tectonic regions with little sediments, such as the South China Block, Greater Khingan range, Qinling range and Changbaishan range (Fig. 3a). In addition, the subduction fronts of the Pacific Plate in Taiwan and the Indian Plate in Nepal also exhibit low velocity patterns, indicating sediment accumulation in those subduction prisms (Fig. S1a).

Figure 3.

(a) Near-surface Vs and (b) its associated uncertainty in the continental China, inferred from the Pwave polarization analysis. The north–south seismic belt and the trans-north China orogen are enclosed with green-dashed lines. Black lines are the block boundaries and black-crossed regions are the major sedimentary basins.

Mean and s.d. of |$V_s^{ns}$| beneath all seismic stations are 3 and 0.70 km s^–1, respectively (Fig. S1a). Near-surface Vs obtained at most of the stations (>93.8 per cent) have an uncertainty smaller than 0.3 km s^–1, or about 10 per cent relatively (Fig. S1b). Uncertainties are larger in the basins with thick Quaternary sediment, such as the Bohaiwan basin and the Taikang-Heifei basin, than in the tectonic regions with little sediment, such as the South China Block, or basins with compacted sediment, such as the Sichuan basin, due to the waveform distortion induced by the soft Quaternary sediment (Fig. 3b).

3.2 Estimations of isotropic ZH ratio and receiver function

3.2.1 Estimation of isotropic Rayleigh wave ZH ratio

We infer ZH ratios from the continuous three-component waveforms of 3848 seismic stations that have at least 1-yr of recording. The seismic data are first corrected for their respective instrumental responses and then bandpass filtered in a frequency range of 0.08–0.3 Hz for the estimation of ZH ratios in 0.125–0.25 Hz, a frequency range that ZH ratios are sensitive to Vs profiles in the top 10 km of the crust. We cut the continuous waveforms into 1-hr segments with a 180 s sliding step between the consecutive segments, and only retain the earthquake-free segments. For each segment and frequency, we define the backazimuth with the maximum projected horizontal spectral amplitude as quasi-radial direction and the horizontal direction perpendicular to the quasi-radial direction as quasi-transverse direction. ZH ratio is obtained by dividing the spectral amplitude of the vertical component by that of the quasi-radial component.

We estimate ZH ratios at each frequency independently. For each frequency, ZH ratios are first estimated based on the seismic data recorded in the trimester centred with each month of the year, yielding isotropic ZH ratios as a function of month. The seasonal component is then removed from the monthly ZH ratio measurements, resulting in a time-independent isotropic ZH ratio. To ensure the stabilization of measurement, the amplitude ratio at a particular frequency is calculated based on the average spectra of the vertical and quasi-radial components over a narrow frequency range of 0.01 Hz centred at that particular frequency, rather than the spectra at the frequency (Tanimoto et al. 2013). We retain the measurements of ZH ratio for further analysis from Rayleigh-wave dominated segments, which are defined as those having a phase shift angle between the quasi-radial and vertical components ranging from |$60^\circ $| to |$120^\circ $|⁠, and an amplitude ratio of the quasi-radial component over the quasi-transversal component (H/T) higher than 3 (Fig. 4a). To avoid contamination of outliers, we compute mean and s.d. of the retained ZH ratios and remove the measurements outside the twofold s.d. regions, with the procedure repeated three times. The retained ZH ratios are then bin stacked with a bin width of |$1^\circ $| in backazimuth to equalize the weight of measurements from different directions. We remove station data that exhibit strong azimuthal variations of ZH ratio, an indication of presence of strong 3-D site effects at the station. The mean of the stacked ZH ratios is defined as the ZH ratio of the month and the s.d. of the stacked ZH ratios is regarded as the uncertainty (Fig. 4b). Generally, the resultant ZH ratios exhibit a seasonal variation (Tanimoto et al. 2006). We remove the seasonal component of the ZH ratios through curve fitting and obtain a time-independent isotropic ZH ratio (Fig. 4c, see support material for details of removing the seasonal component of the ZH ratios).

Figure 4.

Measurement procedure of isotropic ZH ratio for an example seismic station FJ.MXXF in the South China Block, marked as red triangle in top-right insert. (a-c) Isotropic ZH ratio measurement at an example period of 7 s. (a) (Left bottom) Measured phase shift angles and H/T ratios (amplitude ratio of quasi-radial component over quasi-transversal component) at a period of 7 s (circles) from all hourly seismic segments from May 2016 to July 2016, and (top and right) measurement histograms of phase shift angle and H/T value. Only the measurements that are above a H/T boundary (H/T = 3) and Rayleigh-wave dominated, defined as those with a phase shift angle from 60° to 120°n two vertical red dashed lines (blue circles), are retained for ZH ratio analysis. (b) (left-hand bottom) ZH ratios and particle rotation backazimuths of the segments accepted in (a) (blue circles), and average ZH ratios within |${1^\circ }$| bins in backazimuth (red dots) over the selected segments and threefold standard deviation (s.d.) interval of the binned ZH ratios (between red horizontal dashed lines), (right-hand side) histograms of measured ZH ratio (blue) and of binned ZH ratio (red), and (top) histogram of backazimuth. (c) Measured ZH ratios (blue dots) and their associated uncertainties quantified by s.d. of binned ZH ratios (bars) between January 2015 and December 2017 at a period of 7 s, along with estimated seasonal function (red line) and time-independent isotropic ZH ratio (black line). (d) ZH ratio dispersion curve measured from the seismic observations at station FJ.MXXF (blue dots for the measured values and error bars for the uncertainties) composed of isotropic ZH ratios measured in the periods of 4–8 s following the procedure outlined in (a–c). The uncertainty of ZH ratio of each period is the summation of seasonal variation magnitude, s.d. of fitness residual and uncertainty of time-independent isotropic ZH ratio.

We estimate the isotropic ZH ratios in a period range from 4 to 8 s, with an interval of 0.2 s from 4 to 6 s and an interval of 0.5 s from 6 to 8 s. The obtained isotropic ZH ratios over the periods constitute the ZH ratio dispersion curve that is used in the joint inversion for 1-D Vs structure beneath the seismic station (Fig. 4d).

3.2.2 Estimation of isotropic receiver function

We obtain RFs from all tele-seismic event waveforms in the database that are recorded in an epicentre distance range from |${30^\circ }$| to |${95^\circ }$|⁠. Waveforms are corrected for their respective instrumental responses, rotated to the great circle plane and bandpass filtered in a frequency range of 0.05–2 Hz. The filtering scheme is chosen based on the best quality of RFs observed in that frequency range. RFs are computed with waveforms in a time window from 30 s before the predicted P arrival based on PREM and 120 s after the prediction. RFs are obtained by deconvolving the vertical component of the seismic data from the radial component through a time domain iterative deconvolution method (Juan & Charles 1999) with a Gaussian width factor of 5.

We apply following rules for RF quality control and selection. We discard (1) seismic records that have a SNR, defined as the ratio of s.d. of the waveform within 30 s after PREM-predicted P arrival over that 30 s before, smaller than 1.5, (2) RFs obtained for the seismic stations in CNSN and ChinaArray that have a cross-correlation coefficient lower than 0.7 with the mean of all RFs of the station, (3) RFs obtained from the IRIS stations that have a cross-correlation coefficient lower than 0.6 with the mean of all RFs of the station (A relatively high correlation coefficient threshold is used for stations from the CNSN and ChinaArray in the data selection, because of the availability of a large number of high-quality RFs due to the long-time employment of the networks.) and (4) RFs that exhibit strong variations with the backazimuth of the incoming seismic P waves and are likely affected by local 3-D structures.

RFs of a station are amplitude-corrected and bin stacked over all backazimuths. Before stacking, the RF amplitudes are corrected to a reference ray parameter |$\ {p_r} = \ 0.06\ {\rm s\,km}^{-1}$| by multiplying an amplitude correction factor |${f_{\rm amp}}$|  

$$\begin{eqnarray*} {f_{\rm amp}} = {\rm{\ }}\frac{{\tan \left( {2\arcsin \left( {\beta {p_r}} \right)} \right)}}{{\tan \left( {2\arcsin \left( {\beta p} \right)} \right)}}, \end{eqnarray*}$$

(3)

where |$\beta $| is the inverted near-surface Vs from the P polarization analysis and p the theoretical ray parameter predicted with PREM. For seismic stations without reliable estimate of near-surface Vs, |$\beta $| is set to be 3.0 km s^–1. The amplitude-corrected RFs are bin stacked over backazimuth, with a bin width of |${30^\circ }$| and an increment of |${10^\circ }$|⁠, yielding average RFs and their associated s.d.s as a function of binned backazimuth. Finally, an isotropic RF is estimated by stacking mean RFs of bins with a RF number larger than 5, and its uncertainty is quantified by twofold s.d. of those stacked mean RFs. Because this study aims at resolving shallow seismic structure, we only utilize RFs and their uncertainties before the PmS (the converted phase at the Moho) arrivals. We illustrate the RF stacking procedure with an example seismic station FJ.MXXF in the South China Block (Fig. 5a) and another example station T0.K038 in the Bohaiwan basin (Fig. 5b).

Figure 5.

Receiver function (RF) stacking procedure for example seismic stations (a) FJ.MXXF in the south China block and (b) T0.K038 in the Bohaiwan basin. (a) For FJ.MXXF. Top left-hand panel: RFs as a function of backazimuth of incoming P wave. RFs are amplitude-corrected utilizing the near-surface Vs inferred from P wave polarization analysis. Top right-hand panel: bin stacked RFs and their associated uncertainties (black lines) as a function of backazimuth, along with raw RFs in bins (grey lines). Stacking is performed with a bin width of 30° and an increment of 10° in backazimuth, and the uncertainties are estimated from the mean and s.d. of RFs in the bins. Number of RFs in each bin stacking is listed in the right. Bottom left-hand panel: location of seismic station (red triangle) and tele-seismic events (red circles) associated with the RFs in the top left-hand panel. Bottom right-hand panel: the stacked RF (red line) and its uncertainty region (grey area), obtained by stacking the bin-stacked RFs in the top right-hand panel. Green dashed lines mark zero time in all panels. (b) Same as (a), except for station T0.K038.

3.3 Joint inversion of ZH ratio and receiver function for 1-D seismic model

We only search the best-fitting 1-D isotropic seismic structure in the top 20 km of the crust based on joint inversion of the isotropic ZH ratio dispersion curve and RF, with constraint of the near-surface Vs inferred from the P-wave polarization analysis. The sensitivity kernels of ZH ratios indicate that the resolvable depth of Vs structure is about 10 km with a maximum period of ZH ratio of 8 s (Fig. 6b). We parametrize Vs profiles in the top 20 km of the crust in the inversion, and only adopt the results in the top 10 km of the crust. The 1-D profiles of density and P-wave velocity are mapped from Vs profile based on the Brocher's empirical functions (Brocher 2005).

Figure 6.

(a) Model parametrization and (b) sensitivity kernels of ZH ratio. (a) Three-layer, 12-parameter parametrization of 1-D Vs profile (red thick line) beneath each station, with top two layers representing sediments and the third layer representing the upper crust. The velocity profile in the first top layer is parametrized as a linear gradient and those of other two layers are represented by a linear combination of 4 B-spline basis functions (four coloured thin lines at the bottom left). 12 independent parameters include two thicknesses of the top two sedimentary layers, the topmost and bottommost Vs of the first soft sediment layer, four B-spline function coefficients describing Vs profile in the compacted sediment layer and another four representing that of the upper crust layer. (b) Sensitivity kernels of ZH ratio (defined as |$\frac{1}{{ZH}}| {\frac{{\delta ZH}}{{\delta {V_s}}}} |$|⁠) as a function of δVs at depth for four example periods, calculated based on PREM (Dziewonski & Anderson 1981).

3.3.1 Parametrization of shallow velocity profile

A typical 1-D Vs profile is assumed to consist of three isotropic planar layers: a top soft sediment, a bottom compacted sediment and an upper crust layer below the sedimentary layers. The Vs profile in the soft sediment layer has a constant gradient, parametrized with three free parameters including thickness, the bottommost Vs and the topmost Vs. The Vs profiles in the compacted sediment layer and upper crust layer are represented by a combination of four B-spline functions with five independent parameters for the compacted sediment layer and four for the upper crust layer. The five independent parameters for the compacted sediment layer include the thickness of the layer and coefficients of four B-spline functions, while the four free parameters for the upper crust layer are the coefficients of four B-spline functions. The bottom of the upper crust layer is fixed to be 20 km (Fig. 6a). These 1-D layered Vs models are further discretized into 0.5-km-thick sublayers for numerical computations of ZH ratio dispersion curve and RF.

3.3.2 Joint inversion for 1-D Vs structures beneath stations

We adopt Markov Chain Monte Carlo (MCMC) method (Press 1968; Mosegaard & Tarantola 1995; Sambridge 1999; Shen et al. 2013; Goutorbe et al. 2015) to search for 1-D Vs model beneath a seismic station based on joint inversion of the estimated isotropic ZH ratio and RF, starting with an initial model derived from the Rayleigh wave phase/group velocity dispersion data (Shen et al. 2016). MCMC is a global optimization algorithm and is well suited for solving nonlinear inversion problems. Compared with traditional linearized algorithms, MCMC has no need of computing partial derivatives, which is a time-consuming task for the numerical computation of ZH ratio. The method also provides an assemblage of sampled models acceptably fitting the data and properly smoothed, from which the ultimate inverted model and associated model uncertainty can be estimated. In addition, MCMC is nearly independent of the initial model used in the inversion, in contrast to other linearization inversion algorithms in which the inversion results could be greatly affected by the choosing of the initial model (Ojo et al. 2019).

The MCMC inversion is to seek a best-fitting model based on the minimization of an object function defined as:

$$\begin{eqnarray*} {S_{joint}}\left( {\boldsymbol{m}} \right) &=& \alpha {S_{ZH}} + \beta {S_{RF}} + \gamma {S_{SM}} \nonumber \\ &=& \frac{\alpha }{N} \mathop \sum \limits_{i = 1}^N \frac{{{{\left( {ZH_i^{syn}\left( {\boldsymbol{m}} \right) - ZH_i^{obs}} \right)}^2}}}{{\delta _i^2}} \nonumber \\ && + \, \frac{\beta }{M}\mathop \sum \limits_{j = 1}^M \frac{{{{\left( {RF_j^{syn}\left( {\boldsymbol{m}} \right) - RF_j^{obs}} \right)}^2}}}{{s_j^2}} \nonumber \\ && +\, \frac{\gamma }{{L - 1}}\mathop \sum \limits_{k = 1}^{L - 1} \frac{{H({V_k} - {V_{k + 1}})\left\| {\left( {{V_k} - {V_{k + 1}}} \right)} \right\|}}{{0.02}} \end{eqnarray*}$$

(4)

and

$$\begin{eqnarray*} {\rm{H }}\left( {\rm{n}} \right) = \left\{ \begin{array}{ll} 0 & \quad n \le 0 \\ 1 & \quad n \gt 0, \end{array} \right. \end{eqnarray*}$$

(5)

where |${S_{ZH}}$| is the norm-2 misfit between predicted ZH ratio (⁠|$Z{H^{syn}}( {\boldsymbol{m}} )$|⁠) and observed ZH ratio (⁠|$Z{H^{obs}}$|⁠) weighted by the uncertainties (⁠|${\rm{\delta }}$|⁠) of ZH ratio measurements over |${\rm{N}}$| sampled periods, |${S_{RF}}$| is the norm-2 misfit between predicted RF (⁠|$R{F^{syn}}( {\boldsymbol{m}} )$|⁠) and observed RF (⁠|$R{F^{obs}}$|⁠) weighted by the uncertainties of RF data (s) over |${\rm{M}}$| time sampling points, and |${S_{sm}}$| represents an a priori model smoothness term that penalizes presence of low-velocity zones in the model, with 0.02 km s^–1 for normalization. As models consist of 0.5-km-thick finer sublayers after discretization, |${V_k}$| indicates Vs of the kth sublayer from surface. |${\rm{\alpha }}$|⁠, |${\rm{\beta }}$| and |${\rm{\gamma }}$| are the relative weighting parameters between ZH ratio misfit, RF misfit and the model smoothness term, chosen as 100, 100 and 100 in this study by trial and error.

ZH ratio dispersion curve and RF are theoretically computed using a propagator-matrix method (Herrmann 2013) and |${\rm{H\kappa }}$| package (Zhu & Kanamori 2000), respectively, based on the discretized 1-D homogenous isotropic plane layered models.

MCMC sequentially samples models to estimate the posterior probability distributions of the inverted parameters. All plausible models require that (1) parameter values are in the predefined bounds (Table 2) and (2) Vs increments across layer boundaries (the boundary between the soft and compacted sediment layers, and the boundary between the compacted sediment layer and the upper crust layer) are not negative. One iteration of the inversion algorithm involves two steps. (1) The algorithm generates a candidate model |${{\boldsymbol{m}}_{\boldsymbol{j}}}$| near the current model |${{\boldsymbol{m}}_{\boldsymbol{i}}}$| with a uniform probability in the model space. Values of two randomly selected parameters of the candidate model (one thickness and one shear wave velocity parameter) are sampled by random walks with the pre-defined step size and maximum jump size in Table 2, while the values of the rest parameters remain the same as those of the current model |${{\boldsymbol{m}}_{\boldsymbol{i}}}$| (Mosegaard & Tarantola 1995; Goutorbe et al. 2015). (2) The algorithm determines model move |${{\boldsymbol{m}}_i} \to {{\boldsymbol{m}}_j}$| with an acceptable probability |${P_{\rm accept}}$| defined as

$$\begin{eqnarray*} {P_{\rm accept}} = \ {\rm{min}}\left\{ {\frac{{{\cal L}\left( {{{\boldsymbol{m}}_j}} \right)}}{{{\cal L}\left( {{{\boldsymbol{m}}_i}} \right)}},\ 1} \right\}, \end{eqnarray*}$$

(6)

where |${\cal L}\ ( {{{\boldsymbol{m}}_j}} ) = {\rm{\ exp}}( { - {S_{joint}}( {{{\boldsymbol{m}}_j}} )} )$| is likelihood function quantifying the probability of obtaining the measured data with a given specific model |${{\boldsymbol{m}}_j}$|⁠. We generate a random number in the range of |$0 - 1$| with a uniform probability. If the acceptable probability |${P_{\rm accept}}$| is larger than the random number, we accept model move and initiate the next iteration with the candidate model |${{\boldsymbol{m}}_j}$|⁠. Otherwise, we reject model move and start the next iteration with the current model |${{\boldsymbol{m}}_i}$|⁠.

Table 2.

Markov Chain Monte Carlo parameters in joint inversion for 1-D Vs profile.

Parameters	Bound	Max. jump	Step
Soft sediment thickness (km)	0–3	1.5	0.06
Compacted sediment thickness (km)	0–10	2	0.08
Soft sediment topmost Vs (km s–1)	|${m_0} \pm 0.5{m_0}$|	|$0.2{m_0}$|	|$0.01{m_0}$|
Soft sediment bottommost Vs (km s–1)	|${m_0} \pm 0.5{m_0}$|	|$0.2{m_0}$|	|$0.01{m_0}$|
Coefficients of four B-spline functions representing Vs profile in the compacted sediment layer (km s–1)	|${m_0} \pm 0.5{m_0}$|	|$0.2{m_0}$|	|$0.01{m_0}$|
Coefficients of four B-spline functions representing Vs profile in the upper crust layer (km s–1)	|${m_0} \pm 0.5{m_0}$|	|$0.2{m_0}$|	|$0.01{m_0}$|

Note: m₀ is the value of that parameter in the initial model.

Beneath a specific seismic station, the inverted 1-D Vs model and its uncertainty are estimated through statistical analysis of the model assemblage selected from sampled models in 100 000 iterations except for the first 200 iterations (burn-in phase) (Fig. S2a). In order to construct the model assemblage, we find sampled models with the minimal ZH ratio data misfit (⁠|${{\boldsymbol{m}}_{zh}}$|⁠), minimal RF data misfit (⁠|${{\boldsymbol{m}}_{rf}}$|⁠), minimal model smoothness (⁠|${{\boldsymbol{m}}_{sm}}$|⁠) and minimal object function value (⁠|${{\boldsymbol{m}}_{all}}$|⁠), respectively, and compute a normalized object function value |${\varepsilon _{norm}}( {\boldsymbol{m}} )$| for a sampled model |${{\bf m}}$| defined as:

$$\begin{eqnarray*} {\varepsilon _{norm}} \left( {\boldsymbol{m}} \right) = \frac{1}{3} \left\{ {\frac{{{S_{ZH}}\left( {\boldsymbol{m}} \right)}}{{{S_{ZH}}\left( {{{\boldsymbol{m}}_{{\boldsymbol{zh}}}}} \right)}} + \frac{{{S_{RF}}\left( {\boldsymbol{m}} \right)}}{{{S_{RF}}\left( {{{\boldsymbol{m}}_{{\boldsymbol{rf}}}}} \right)}} + \frac{{{S_{SM}}\left( {\boldsymbol{m}} \right)}}{{{S_{SM}}\left( {{{\boldsymbol{m}}_{{\boldsymbol{sm}}}}} \right)}}} \right\}. \end{eqnarray*}$$

(7)

Sampled models with a normalized object function value smaller than |${\varepsilon _{norm}}( {{{\boldsymbol{m}}_{all}}} ) + 1$| are retained as members of the model assemblage. All these models are regarded as acceptably satisfying the constraints of the observed ZH ratios, RF data sets and model smoothness. Mean Vs profile and model uncertainty (quantified by s.d. of the accepted models) are computed among the model assemblage (Figs S2b–S2d). Similarly, we then estimate the depths and associated uncertainties of two discontinuities (bases of the top two layers) with the model assemblage (Figs S2b, S2e and S2f). The sampled model that is the nearest to the mean Vs profile, is taken as the inverted Vs model beneath the station (Fig. S2b).

The joint inversion resolves seismic structures well for stations located in various geological settings. We show examples in four typical geological settings: (1) a basin region with thick soft sediment where RF has a strong reverberation (Fig. 7a), (2) a basin region with soft and compacted sediments where RF has a small-amplitude direct P phase at time zero, followed by a positive phase and a negative phase (Fig. 7b), (3) a basin region with compacted sediment where RF has a large-amplitude direct P phase at time zero, followed by a positive phase and a negative phase (Fig. 7c) and (4) a hard-rock region with no sediment where RF has a single direct P phase at time zero (Fig. 7d).

Figure 7.

1-D Vs inversion results for example seismic stations in four geological units: (a) T0.K038 (Bohaiwan basin), (b) X2.289 (Ordos basin), (c) CQ.FUL (Sichuan basin) and (d) FJ.MXXF (South China Block). (a) Inverted Vs models, ZH ratio fitting and RF fitting for T0.K038. Left-hand panel: initial model (green line), inverted model (red line), average model (black dashed line) and threefold s.d. interval (grey shaded area) of the acceptable models, along with the inferred soft sediment base (purple line) and its uncertainty region quantified by s.d. of the soft sediment base depths of the acceptable models (purple shaded area). Right-hand top panel: estimated RF region from the observation (area enclosed by blue lines), synthetic RFs computed with the inverted (red) and initial (green) models and threefold s.d. interval (grey area) of synthetic RFs from all acceptable models and their average (black dashed line). Right-hand bottom panel: observed isotropic ZH ratios along with uncertainties (blue dots with bars), predicted ZH ratio curves based on the inverted (red) and initial (green) models and threefold s.d. interval of synthetic ZH ratio curves (grey shaded area) of all acceptable models and their mean (black dashed line). (b–d) Same as (a), except for stations X2.289, CQ.FUL and FJ.MXXF, respectively. Specifically, depth of the compacted sediment base (blue line) and associated uncertainty (blue shaded area) for X2.289 are shown in the left of (b). (e) Station locations (blue triangles) labeled with names. Symbols have the same meanings as those in Fig. 1(a).

Inverted Vs models exhibit major characteristics that fits the specific features of the RF and ZH ratio (Fig. 7). For examples, (1) a near-surface Vs gradient is present in most of the inverted velocity profiles, corresponding to the relatively large observed width of the first energy pulses in RFs and small ZH ratio amplitude at the relative short periods. (2) The average Vs and thickness of the sedimentary layer govern the apparent time lag from the time zero of the first RF energy pulse (a supposition of direct P wave and Ps conversion from the base of the sedimentary layer). Models with a larger average Vs and/or a thinner sediment thickness correspond to the RF observations with a smaller time-shift of the first energy pulse. (3) The Vs contrasts across the sediment-crust discontinuity contributes to the positive (direct phase) and negative (multiple phase) secondary phases in RFs. Models of larger Vs contrasts correspond to the RF observations with larger amplitudes of the secondary phases (Figs 7a and c). Models with a very large Vs contrast and a thin sediment thickness correspond to the RFs with strong reverberations (Fig. 7a). (4) Models with a higher Vs at the Earth's subsurface correspond to larger RF amplitudes of direct P phases observed at time zero. (5) The inferred Vs gradients at different depth correspond to the observed ZH ratio trends in different period ranges. For instance, Vs profile beneath seismic station X2.289 has different Vs gradients in the soft sediment layer (0–2.6 km) and the compacted sediment layer (2.6–3.9 km), in accordance with the slightly downward (4–6 s) and upward (6–8 s) trends in the ZH ratio dispersion curve (Fig. 7b). (6) ZH ratios of longer periods (⁠|$\ge $|7 s) are not fitted as well as those of shorter periods (⁠|$\le $|6 s) due to larger measurement uncertainties and fewer sampling points (Figs 7a and d).

We perform resolution tests to check the robustness of the inversion results and the dependency of the inversion results on the initial model used in the inversion. We show an example using a basin station T0.K001. We use the inverted 1-D Vs model of the station as input model and forward calculate synthetic ZH ratio and RF. Gaussian noises are added into the synthetic ZH ratio and RF. Inversions of the synthetics with different levels of noise yield consistent inverted models that are nearly identical to the input model, indicating that our inversion approach is robust in the presence of noise up to a level of 20 per cent (Fig. S3). Inversions with two distinctly different initial models, a constant velocity model and a depth-varying model, also yield consistent final models that are similar to the input model, showing that our inversion approach is nearly independent of the initial model (Fig. S4).

4 SHALLOW SEISMIC STRUCTURE AND SEDIMENT THICKNESS BENEATH THE CONTINENTAL CHINA

The shallow seismic structure and sediment thickness in the study region are obtained through interpolation from the best-fitting 1-D seismic structures beneath the stations.

4.1 Shallow 3-D Vs model

We obtain a 3-D Vs model in the continental China by interpolating those inverted 1-D Vs profiles beneath 3350 seismic stations with a simple kriging algorithm (Schultz et al. 1999), resulting in gridded 2-D Vs maps at various depths (see support material for details of the interpolation scheme). The interpolation is performed with a grid size of |$\ {0.3^\circ } \times {0.3^\circ }$| in two regions with dense station coverage including the north–south seismic belt and the trans-north China orogen, and a grid size of |${0.5^\circ } \times {0.5^\circ }$| in the rest of the continental China (except the Tarim basin and the southwest Tibet where no station coverage is present, Fig. 8).

Figure 8.

Left-hand panel: shallow 3-D Vs model and (right-hand panel) its associated uncertainty in the continental China at four selected depth intervals. (a) Lateral variation of the Vs averaged in a depth interval of 0.5–1.5 km, along with block boundaries (thin black lines), major sedimentary basins in China (black-crossed regions). The north–south seismic belt and the trans-north China orogen are enclosed with green-dashed lines. (b) Uncertainty of Vs in (a), along with locations of the stations with reliable inverted 1-D Vs profiles (blue triangles). (c,d), (e,f), (g,h), same as (a,b), except for the averages in the depth intervals of 2.5–3.5 km, 6.5–7.5 km and 9.5–10.5 km (shown at bottom left of all panels), respectively.

In the shallow 3-D Vs model, high Vs is evident at all depths for the regions with little sedimentary deposits (Fig. 8). These regions include some old tectonic blocks with relatively higher velocities such as the South China Block, and orogenic belts with relatively lower velocities such as the Greater Khingan range, Lesser Khingan range, Changbaishan range, Dabie-Sulu orogenic belt, Qinling range and western Sichuan plateau (Fig. 8). The relative velocity differences of the bedrocks between two types of tectonic regions are likely caused by their formation processes, with the old tectonic blocks consisting of ancient microcontinents and the orogenic belts induced by the interactions between the continental blocks (Li et al. 1999; Zhou & Li 2000; Meng & Zhang 2000b; Royden et al. 2008; Jia et al. 2013).

The most prominent velocity features in the shallow crust are low Vs anomalies in the sedimentary basins, manifesting the seismic signatures of sedimentary deposits in the basins. At 0.5–1.5 km depth, low Vs anomalies are evident in major sedimentary basins, including the Junggar, Tarim, Qaidam, Qiangtang-Tanggula, Bohaiwan, Taikang-Hefei, Lanping-Simao, Chuxiong, Sichuan, Jiuquan-Minle-Wuwei, Jianghan, Songliao, Ordos, Subei-Yellow sea and Erlian basins (Fig. 8a). The sediments deposited in different geological periods exhibit different magnitudes of low-velocity anomalies. The Sichuan basin (excluding the Chengdu plain), where sediments were mainly deposited during the Mesozoic (Wang et al. 2016), exhibits a relative high average subsurface Vs of about 2.4 km s^–1, while other basins with Quaternary deposits, such as the Bohaiwan basin, have an average subsurface Vs lower than 2 km s^–1 (Fig. 8a). This is likely because the shallow sediments in the Mesozoic basins are more compacted than those of the Quaternary basins. At about 3 km depth, low Vs velocities appear in major sedimentary basins, such as the Songliao, Bohaiwan, Sichuan, Ordos, Qaidam, Junggar and Tarim basins, while other basins exhibit Vs of crustal bedrocks (Fig. 8c). At ∼7 km depth, parts of the Sichuan, Ordos, Qaidam, Songliao, Junggar and Tarim basins still manifest sedimentary structures with an average Vs of about 3 km s^–1, while the Bohaiwan basin show high Vs anomalies corresponding to its cratonic basement (Fig. 8e). At about 10 km depth, only the northern Tarim, southern Junggar and northwest Sichuan basin still exhibit seismic structure of a sedimentary basin with Vs lower than 3 km s^–1, while other basins display velocity structures of the crystalline crustal bedrocks (Fig. 8g).

4.2 Thickness of sediments

We construct maps of sediment thickness in the continental China following a similar procedure. We first identify the soft and compacted sediment layers in each of the 3350 inferred 1-D Vs models beneath the stations, and then interpolate the inferred sediment thicknesses into 2-D grids in the continental China. We define the soft sediments as young unconsolidated deposits and the compacted sediments as old partially consolidated deposits. While the top two layers of the 1-D models are designed for characterizing seismic structures of the soft and compacted sediments respectively, we need a separate judgment to confidently identify the soft and compacted sediments in the top two layers of the inverted seismic models. In reality, each of the top two layers of an inverted seismic model could be a soft sediment layer, a compacted sediment layer or a crustal bedrock layer (Fig. 7). We adopt this approach to classify the characteristics of the top two layers. We first derive empirical separation lines between the layers of the soft sediment, compacted sediment and crustal bedrock based on the inferred seismic characteristics beneath the stations where the site geology is known (Figs 9a and b). We then apply those empirical separation lines to classify the seismic characteristics of the top two layers of the inferred models beneath other stations. Our analysis indicates that layers of the soft sediment, compact sediment and crustal bedrock are well separated by the equivalent Vs at the surface (⁠|$V_s^0)\ $| and the mean Vs gradient inside the layer.

$$\begin{eqnarray*} V_s^0 = \overline {{V_s}} \ - \bar{h}\frac{{{\rm d}{V_s}}}{{{\rm d}h}}, \end{eqnarray*}$$

(8)

where |$\overline {{V_s}} $| is the mean Vs, |$\bar{h}$| is the mean depth and |$\frac{{{\rm d}{V_s}}}{{{\rm d}h}}$| is the mean Vs gradient of the layer. We derive the linear relationships with a Support Vector Machine (SVM) method (Pedregosa et al. 2011), which correctly identify the sedimental materials in a 94 per cent accuracy (Figs 9c and d). The separation lines are:

$$\begin{eqnarray*} V_s^0 + 0.06{\rm{\ \ }}\frac{{{\rm d}{V_s}}}{{{\rm d}h}} = \ 1.47\ \left( {\rm km\,s}^{-1} \right) \end{eqnarray*}$$

(9)

separating layers of the soft and compacted sediments and

$$\begin{eqnarray*} V_s^0 + 0.12{\rm{\ \ }}\frac{{{\rm d}{V_s}}}{{{\rm d}h}} = \ 2.3\ \left( {\rm km\,s}^{-1} \right) \end{eqnarray*}$$

(10)

Figure 9.

Known material types in the depths of top two layers and classifiers of material types in the top two layers of inverted 1-D Vs models. (a–b) Maps of known soft sediment (orange dots), compacted sediment (green dots) and crustal bedrock (blue dots) at the depths of the (a) first and (b) second top layers of inverted models, plotted on the sites of seismic stations, along with the block boundaries (solid grey lines), plate boundaries (blue lines) and major sedimentary basins (grey-crossed areas). (c–d) Classifiers (colour dashed lines) of soft sediment, compacted sediment and crustal bedrock in the (c) first and (d) second top layers of inverted models (dots coloured with material type), based on attributes of Vs at the surface (⁠|$V_s^0$|⁠) and Vs gradient of the layer (⁠|$\frac{{{\rm d}{V_s}}}{{{\rm d}h}}$|⁠). The classifier line for soft and compacted sediments is |$V_s^0 + 0.06{\rm{\ \ }}\frac{{{\rm d}{V_s}}}{{{\rm d}h}} = \ 1.47\ ( {\rm km\,s}^{-1} )$| (red dashed line), while that for compacted sediment and crustal bedrock is |$V_s^0 + 0.21{\rm{\ \ }}\frac{{d{V_s}}}{{dh}} = \ 2.30\ ( {\rm km\,s}^{-1} )$| (black dashed line).

separating layers of the compact sediment and crustal bedrock. After the layers of the sediments are identified in those individual models beneath the seismic stations, maps of the thicknesses of the soft sediment (Figs S5a and S5b), compacted sediment (Figs S5c and S5d) and total sediment (Figs S5e and S5f) are generated in the continental China. Models of the sediment thicknesses are produced through interpolation from the inferred sediment thicknesses beneath the 3350 seismic stations with the simple kriging algorithm (Fig. 10, see support material for more details of the interpolation scheme).

Figure 10.

Inferred thicknesses of soft sediment, compacted sediment and total sediment in the continental China. Lateral variation of (a) thickness and (b) its associated uncertainty of the soft sediment, along with major sedimentary basins (black-crossed regions) and block boundaries (solid lines). The north–south seismic belt and the trans-north China orogen are enclosed with the green-dashed lines. (c, d), (e, f), same as (a, b), except for the compacted sediment and total sediment (including both the soft and compacted sediments).

The soft sediments are concentrated in central regions of the rift basins such as the Bohaiwan, Songliao and Qaidam basins with a thickness greater than 2 km, and in regions of the foreland basins such as the northern margin of the Ordos basin with a thickness greater than 3 km (Fig. 10a). Those soft sediments likely reflect the unconsolidated deposits of Cenozoic age. Thick compacted sediments are distributed in the old craton basins such as the Sichuan basin and the Ordos basin, and the foreland basins such as the northern and southwestern margins of the Tarim basin and the southern margin of the Junggar basin, with a thickness greater than 5 km (Fig. 10c). Those compacted sediments probably reflect the partially consolidated sediments of Mesozoic age. Overall, thick sediments are present in major sedimentary basins with the largest thickness reaching about 10 km in the southwest Tarim basin, some sediments in the intermontane basins in young orogenic belts such as the Taihang range with a thickness of about 2 km, and little sediments in the old tectonic blocks such as the South China Block with a thickness smaller than 1 km (Fig. 10e).

4.3 Seismic effects of basin compaction and sediment composition

Seismic properties of the sediments manifest the compaction effect of basin environment and the compositional difference of depositions. Compaction effect in a basin may result in a change of intrinsic Vs and velocity gradient with depth, while the compositional difference of the sediments during deposition will result in different intrinsic seismic velocities and velocity gradients. The compaction effects are evident for the soft and compacted sediments in the seismic model. For example, the soft sediments in the Bohaiwan basin transit from a small surface Vs and a large Vs gradient in the first top layer to a larger equivalent surface Vs and a smaller Vs gradient in the second top layer (Fig. 11a), whereas the compacted sediments in the Sichuan basin exhibit a similar transition of seismic characteristics from the first to the second top layer (Fig. 11b). Composition effects of the sediments are also evident in the seismic model and can be best illustrated between the soft sediments in the Chengdu plain in western Sichuan basin which are dominated by conglomerates of igneous and metamorphic clasts from Precambrian sources (Burchfiel et al. 1995) and those in the Bohaiwan basin which mainly include the silty clay and sandstone (Ye et al. 1985). Although both are of Quaternary age, the soft sediments in the Chengdu plain exhibit a larger surface Vs and a larger Vs gradient than those in the Bohaiwan basin, due to their compositional difference (Fig. 11c).

Figure 11.

Examples of compaction and composition effects on seismic structure of sediments. (a) Compaction effect of the soft sediments in the Bohaiwan basin. Left-hand panel: map of the soft sediments at the depths of the top two layers in inverted models, plotted on the sites of stations (orange dots). Right-hand panel: the surface Vs and Vs gradient of the soft sediments at the depths of the first (red hollow circles) and second (blue hollow circles) top layers of inverted models beneath seismic stations shown in the left-hand panel, with their medians and s.d.s indicated by solid circles with error bars (red one for the first top layer and blue one for the second top layer). (b) Same as (a), except for the compaction effect of the compacted sediments in the Sichuan basin. (c) Composition effect of the soft sediments of Quaternary age in the Bohaiwan and Sichuan basins. Left-hand panel: map of the soft sediments at the depths of the first top layers in inverted models, plotted on the sites of stations in the Bohaiwan (orange dots) and Sichuan (orange triangles) basins. Right-hand panel: the surface Vs and Vs gradient of the soft sediments in the first top layers of inverted models beneath stations in the Bohaiwan (red hollow circles) and Sichuan (red hollow triangles) basins shown in the left-hand panel, along with their medians and s.d.s (red solid circle with error bars for the Bohaiwan basin and red solid triangle with error bars for the Sichuan basin). (a–c) In the left-hand panels, regions of the Bohaiwan and Sichuan basins (red-crossed areas) are coloured with the surface topography and indicated by the blue enclosed region in the inserts. Classifiers (colour dashed lines) of the soft sediment, compacted sediment and crustal bedrock in the top two layers of inverted models are shown in the right-hand panels, which are the same as those in Fig. 9(c).

4.4 Uncertainties and resolution of seismic models

Uncertainties of the seismic models are estimated by an interpolation procedure similar to that in producing the 3-D Vs model and maps of sediment thicknesses. Uncertainties of the 1-D Vs profiles and sediment thicknesses are firstly inferred from joint inversion results of seismic stations, and then interpolated into 2-D grids with a simple Kriging algorithm (see support material for more details of the interpolation scheme). Uncertainties of the 3-D Vs model are mostly smaller than 300 m s^–1, or less than 10 per cent, over all depths (Figs 8b, d, f and h), while those of the sediment thicknesses are also less than 10 per cent in most areas (Figs 10b, d and f).

The spatial resolution of the seismic models is determined to be 0.6–1.2° in the north–south seismic belt and the trans-north China orogen, and 1–2° for those in the rest of the continental China (except the Tarim basin and the southwest Tibet with no station coverage). The resolution is quantified as the minimum distance that two |$\delta $|-shape anomalies can be resolved in the interpolated seismic models. The resolution of the seismic models varies from region to region, depending on the station geometry around the model grid. Poor station coverage leads to some model artifacts. For example, the low velocity feature appearing in the Tianshan range is resulted by the smearing effect in the interpolation of the seismic profiles beneath the nearby basin stations (Fig. 8a).

4.5 Comparisons with results from other studies

For comparisons, we collect velocity profiles measured in drilling wells and sediment basements inferred from active-source studies across the continental China in the published literatures. We are able to find 1-D P-wave velocity (Vp) profiles in seven drilling wells located in the Bohaiwan basin (Zhang et al. 2018), Ordos basin (Li et al. 2016a), Sichuan basin (Yao & He 2006), North China craton (NCC) (Zhu et al. 2008), Songliao basin (Li et al. 2011), Qaidam basin (Liu et al. 2019) and Tarim basin (Li et al. 2015, labelled 1–7, respectively, in Fig. 12a) and four seismic profiles in active-source studies located in the North China craton (Jia et al. 2014), Sichuan basin (Wang et al. 2016), Songliao basin (Feng et al. 2010), and Tarim basin (Zhang et al. 2001, labelled as a–a’, b–b’, c–c’ and d–d’, respectively, in Fig. 12a). We estimate Vs profiles in these drilling wells from the Vp measurements based on an empirical relation between the two velocities (Brocher 2005).

Figure 12.

Comparisons of the inferred 3-D Vs model and sediment basement with published results. (a) Locations of drilling wells (red squares) with borehole Vs estimated from Vp measurements presented in (b) and active source profiles (bold black lines) with results of sediment basement presented in (c). Blue triangles indicate seismic stations nearest to the drilling wells, while grey triangles show stations within a distance of 1° from the seismic profiles. Other symbols have the same meanings as those in Fig. 1(a). (b) Comparison of 3-D Vs model at the sites of drilling wells (black lines), inverted 1-D profiles beneath the nearest seismic stations (blue lines) and estimated Vs profiles in the drilling wells (red lines), with panels labeled accordingly with the drilling wells in Fig. 12(a). (c) Comparison of the inferred sediment basements in this study (black bold lines) with published results (red dashed lines), with panels labelled accordingly with profiles in Fig. 12(a). Grey-shaded areas indicate topography along profiles. Note the disproportional scales in horizontal distance, depth and elevation.

The shallow 3-D Vs model is in accordance with those Vs profiles in most of the drilling wells (Fig. 12b). For the wells in the Ordos basin, Sichuan basin, North China craton and Tarim basin, the inferred Vs profiles from the Vp measurements are nearly identical at a depth larger than 0.5 km with the inverted Vs profiles beneath their nearest seismic stations with a distance smaller than 40 km, and are similar with the 3-D Vs model at the drilling sites. For the wells in the Qaidam basin and Songliao basin, the Vs profiles in the wells are comparable at a depth larger than 2 km to the inverted Vs profiles beneath their nearest seismic stations with a distance larger than 60 km, and to the shallow 3-D Vs model at the drilling sites. One notable discrepancy between the results lies in the Nanpu 5–98 well located in the Bohaiwan basin. The estimated Vs profile is systematically smaller than the 3-D Vs model at the site, but it is comparable to the inverted Vs profile beneath its nearest seismic station (∼20 km away) at a depth range of 3–6 km. Such consistency and difference are likely caused by the interpolation of the seismic results in a region with rapid change of geological structure, but with insufficient seismic station coverage.

Our inferred sediment thicknesses are comparable to those from the published active-source seismological studies in the regions with dense seismic station coverage and smooth change of geological structure (Fig. 12c). The differences between the inferred sediment thicknesses and those from the active-source studies are smaller than 1.5 km in the central regions of the Bohaiwan basin, Ordos basin, Sichuan basin, and Songliao basin, and the young orogenic belts of the Taihang range and the Ludong-jiaobei uplift (Fig. 12c). Larger differences of inferred sediment thickness exist between our study and the published results in the Kashgar depression in the southwest marginal regions of the Tarim basin where our inferred sediment basement is about 3 km shallower than the published result of 11 km and in the southwest Sichuan basin where our inferred sediment basement is about 2 km shallower than the published result of 6 km (Fig. 12c). While some of the discrepancies may be explained by the sparse data coverage and rapid change of geological structure in those regions, it could also be possible that the sediment basement is defined differently between the studies.

5 IMPLICATIONS OF BASIN EVOLUTION FROM SHALLOW SEISMIC STRUCTURE

Our inferred seismic structures in the basins would be expected to provide significant insights into the evolution of the major basins in the continental China. While it is beyond the scope of this paper to use our seismic results to quantitatively constrain the basin evolution, we discuss possible relationship of our inferred seismic results in the major basins with their tectonic history and geological evolution.

The Songliao basin is a typical rift basin developed inside an old craton and has been accumulating deposits since the Jurassic. The Songliao basin is covered by a layer of soft sediment (∼1.5 km thick) and a layer of compacted sediment (∼4 km thick). Basement of the sediments exhibits a classical dish-like shape of a rift basin, with the maximum thickness of the sediments reaching ∼6 km at the central downwarp zone (Figs 10 and 13a). These estimates are similar with the results from previous active-source seismological study (Fig. 12c, Feng et al. 2010). The sediments in the basin likely comprise mudstone, siltstone and sandstone of lacustrine, fluvial and deltaic origins, which were deposited in the episode of thermal subsidence since the Lower Cretaceous, induced by subduction of the Okhotsk seafloor. The soft and compacted sediments are probably separated by an unconformity formed between the early and the late period of thermal subsidence in the Upper Cretaceous (Fig. 13a, Feng et al. 2010).

Figure 13.

Cross-sectional views of shallow shear velocity structure and typical 1-D Vs models of three geological units in the eastern continental China. (a) For the region of the Songliao basin. Top left-hand panel: the region of the Songliao basin (red-crossed area) is coloured with the surface topography and indicated by the blue enclosed region in the insert. Locations of a demo site (green triangle) and two cross-sections A—A’ and B—B’ (black dashed lines), with the 1-D seismic model of the demo site shown in the top right-hand panel and the shallow velocity profiles of the cross sections shown in the bottom panels. Top right-hand panel: the inferred 1-D seismic model of the demo site (a green triangle in the top left-hand panel), including Vs profile (blue line) with its associated uncertainty range (red dashed lines), soft sediment basement (green line) with its associated uncertainty range (green-shaded area), and total sediment basement (black line) with its associated uncertainty range (black-shaded area). Bottom-hand panel: the shallow seismic structures along cross-sections A—A’ and B—B’, showing Vs structure (coloured background), and bases of soft sediment (bold green line) and total sediment (bold black line), along with the surface topography (grey area) and location of the demo site (green triangle) on the top of the profiles. Major geological units are labelled. (b–c) Same as (a), except for (b) the region of the Bohaiwan basin and (c) the region of the South China Block. Note the disproportional scales in horizontal distance, depth and elevation.

The Bohaiwan basin is also a rift basin inside an old craton receiving deposits since the Mesozoic. The Bohaiwan basin also has a dish-like basement of the sediments, with the maximum thickness of the sediments approaching ∼7 km at the Bozhong sag (Zhang et al. 1994; Zhao et al. 1999). The soft and compacted sediments in the basin probably comprise mudstone, siltstone and silty clay of fluvial and lacustrine origins, deposited in several cycles of rifting and thermal subsidence in the Cenozoic era (Ye et al. 1985). Rifting of the craton basement was possibly driven by both the asthenospheric upwelling and the interactions between the Eurasia, the Pacific and the Indian plates during the late Mesozoic (Ren et al. 2002). The large Vs velocities in the crustal layer indicate a cratonic bedrock of the basin, also supported by a large Vs velocity contrast between the sediments and the underlying crustal bedrocks (Fig. 13b).

The Sichuan basin is a basin inside a craton with a long complex evolution history, with its earlier basin formation in a time period from the Precambrian to the Triassic and a later foreland basin deposit since the Cretaceous. The compacted sedimentary layer of the seismic model likely reflects the earlier formation deposits before the Triassic, while the soft sedimentary layer likely reflects the deposits of the foreland basins since the Cretaceous. The sediment thickness of the Sichuan basin decreases eastward from ∼7 km at the Longmenshan depression to ∼3 km at the central uplift (Figs 10 and 14a), comparable to the results from the drilling wells and seismic reflection profiles (Wang et al. 2016). Most areas of the basin are covered by compacted sediments with thickness reaching ∼5 km (Fig. 10), indicating extensive shallow-marine and non-marine sediments of the basin before the Triassic. Soft sediments of the Sichuan basin are concentrated at east of the Longmenshan fault (the Chengdu plain) and south of the Micang range (northwest Sichuan basin) with a thickness of ∼1 km. This geographic distribution of soft sediments is consistent with a scenario that the foreland basins are developed above foreland depressions, caused by loadings of the western Sichuan plateau and the Qingling range (Figs 10 and 14a). Furthermore, the soft sediments in the Chengdu plain likely reflect the fluvial-lacustrine conglomerate, dominated by igneous and metamorphic clasts from Precambrian sources in the Quaternary (Burchfiel et al. 1995), while those in the northwest Sichuan basin are probably deposits of the fluvial-lacustrine sedimentation in the Cretaceous (Meng et al. 2005).

Figure 14.

Cross-sectional views of shallow shear velocity structure and typical 1-D Vs models of four major sedimentary basins in the Circum-Tibetan Plateau Basin-Range System. (a–d) Same as the (Fig. 13a), except for the regions of (a) the Sichuan basin, (b) the Ordos basin, (c) the Tarim basin and (d) the Qaidam basin.

The Ordos basin is a basin inside the Ordos craton that was developed in multicycle depositions since the Palaeozoic. The Ordos basin has a saucer-shaped base, with the maximum thickness of the sediments reaching ∼7 km at the centre (Figs 10 and 14b). The old compacted sedimentary layer (∼4 km) and the young soft sedimentary layer (∼2 km) reflect continuous deposits of the basin since the Palaeozoic. The soft sediments around the Ordos basin are likely related to the deposits in the foreland depressions caused by the stretching of the nearby tectonic regions during the Jurassic, such as the Yinshan range (Figs 10 and 14b), while the compacted sediments likely reflect deposits inside the Ordos craton since the Palaeozoic. The crustal bedrocks of the Ordos basin exhibit relatively larger Vs velocities than the crustal bedrocks in some young tectonic regions such as the Qinling range, manifesting its cratonic nature (Fig. 14b, Yang et al. 2005; Zhang et al. 2007).

The Tarim basin is a sedimentary basin developed inside a Precambrian craton that has been accumulating deposits since the Palaeozoic, with its marginal regions reformed by foreland basins since the Mesozoic. Limited by the poor spatial coverage of seismic stations, only the seismic structures in some marginal regions of the basin are resolved. Among the regions are the Kaqu depression with a sediment thickness of ∼8 km, the Kashgar depression with a sediment thickness of ∼7 km and the Keping uplift with a sediment thickness of ∼2 km (Figs 10 and 14c). These two depressions in the foreland basins and the Keping uplift in the centre are the results of the loadings of the Kunlun range and the Tianshan range since the Mesozoic (Sobel 1999; Zhang et al. 2011; Jia et al. 2013). The Kaqu depression consists of a thin layer of soft sediment (∼0.5 km thick) and a thick layer of compacted sediment (∼8 km thick), likely reflecting deposition since the Mesozoic (Yan et al. 2003). The Kashgar depression also consists of layers of the soft and compacted sediments (Fig. 14c), likely representing the deposits since the Mesozoic (Cheng et al. 2012; Tang et al. 2015). The relative low Vs velocities of the upper crust layer in the Kaqu and Kashgar depressions possibly reflect the marine sediments deposited during the Palaeozoic (Figs 8g and 14c, Tang et al. 2015). The Keping uplift consists of layers of soft sediment (∼0.5 km thick) and compacted sediment (∼2 km thick, Fig. 14c). The compacted sedimentary layer probably reflects the marine sediments deposited during the Palaeozoic, while the soft sedimentary layer likely indicates deposits inside the thrust nappe structure since the Cenozoic. The relative high Vs velocities of the bedrocks in the Keping uplift suggest a Precambrian crystalline basement (Fig. 14c, Zhang et al. 2001).

The Qaidam basin was developed in a foreland depression in the Jurassic and later experienced a rifting process in the Cenozoic possibly caused by an upwelling in the upper mantle (Xia et al. 2001). Sediment thickness of the Qaidam basin increases southward and reaches ∼5 km before the east Kunlun range, consistent with the scenario that the basin was formed in a foreland basin above a depression caused by loading of the Kunlun range in the Jurassic (Fig. 14d). The soft sedimentary layer likely reflects sediments of the rift basin in the Cenozoic, while the compacted sedimentary layer probably indicates deposits of the foreland basin in the Jurassic (Figs 10 and 14d). Crustal bedrocks of the Qaidam basin exhibit Vs velocities comparable to those of the nearby orogenic belts such as the east Kunlun range and smaller than those of the cratonic areas such as the South China Block, indicating the basin was developed above a non-cratonic crust (Figs 8g, 13c and 14d, Xia et al. 2001; Yang et al. 2004).

6 CONCLUSIONS

We construct a high-resolution shallow 3-D seismic model in the top 10 km of the upper crust in the continental China, using constraints of P-wave polarization, Rayleigh wave ellipticity and receiver function obtained from the seismic data recorded by 3848 stations in multiple networks deployed across China since 1990 (the CNSN, ChinaArray and networks from the IRIS). Our shallow 3-D seismic model has a spatial resolution of 0.6–1.2° in the north–south seismic belt and the trans-north China orogen, and a resolution of 1–2° at the rest of the continental China (except the Tarim basin and the southwest Tibet). Our shallow 3-D seismic model exhibits low velocity anomalies of deposits in major sedimentary basins such as the Songliao, Bohaiwan, Sichuan, Ordos, Tarim and Qaidam basins and high velocity anomalies of crustal bedrocks in young orogenic belts such as the Taihang range and the Greater Khingan range and old tectonic blocks such as the South China Block and the Ordos craton. Our inferred sediment thickness maps display thick deposits in major sedimentary basins such as the soft sediments of Cenozoic age in the Bohaiwan basin with a thickness of about 4 km and compacted sediments of Mesozoic age in the Tarim basin with a thickness of about 8 km, some compacted sediments in the intermontane basins in young orogenic belts such as the Taihang range with a thickness of about 2 km and little sediments in the old tectonic blocks such as the South China Block with a thickness lesser than 1 km. Compaction effects are evident in the sediments of the major sedimentary basins, with their seismic profiles transit from a small intrinsic velocity and a large velocity gradient in shallow depth to a larger intrinsic velocity and a smaller velocity gradient in large depth. We also discuss implications of tectonic history and geological evolution of the basins based on the inferred geometries from the seismic model, such as dish-like bases of the Songliao and Bohaiwan basins, a westward dipping base of the Sichuan basin, a southward dipping base of the Qaidam basin, and a saucer-shaped base of the Ordos basin.

Our seismic results can be used as constraints on the composition of the bedrocks and sediments in the continental China and as a reference of shallow seismic structure for seismic hazard assessment and future high-resolution studies of deeper seismic structure. This study would provide an effective mean for the joint inversion of various seismic constraints to obtain the seismic structure and establishes a framework of seismic data sharing for future studies in the seismological community in a first step of developing a China Seismological Reference Model.

SUPPORTING INFORMATION

Figure S1. (a) Inverted near-surface shear wave velocity and (b) its associated uncertainty, obtained from P-wave polarization analysis in individual stations and plotted at the sites of seismic stations. Grey shaded areas represent major basins, while the regions enclosed by black lines are geological blocks.

Figure S2. Statistical analysis of seismic models sampled in the Markov Chain Monte Carlo (MCMC) inversion for an example seismic station FJ.MXXF in the South China Block. (a) Misfit function of 100 000 sequentially sampled models (red line). Models sampled in the burn-in phase (blue-shaded area), from the 1st to 200th iteration, are discarded in the statistical analysis. (b) Inverted 1-D seismic model (red line), the mean (black dashed line) and threefold s.d. region (grey-shaded area) of all acceptable models, along with the mean (purple diamond) and s.d. (purple error bar) of the base depth of the first top layer, and the mean (green diamond) and s.d. (green error bar) of the base depth of the second top layer. (c) Distributions of the a prior probability (red stair-steps) and the a posteriori probability (blue histogram) of Vs at 3 km in depth estimated from the model assemblage, along with the mean (black diamond) and threefold s.d. (error bar) of Vs at 3 km in depth estimated from the distribution of the a posteriori probability. (d) Same as (c), except for Vs at 7 km in depth. (e) Distributions of the a prior probability (red stair-steps) and the a posteriori probability (blue histogram) of the base depth of the first top layer, along with the mean (purple diamond) and s.d. (purple error bar) calculated from the distribution of the a posteriori probability. (f) Same as (e), except for the base depth of the second top layer.

Figure S3. Model recovery tests based on synthetics with noise levels of (a) 5 per cent, (b) 10 per cent, (c) 15 per cent and (d) 20 per cent. In each subfigure, (left) inverted model (red line) and its associated uncertainty (grey-shaded area), along with the input model (green line), (right-top) synthetic RF calculated with the input model (green line) and its uncertainty range with noise added (two blue lines), predicted RF based on the inverted model (red line), and threefold s.d. interval of RFs predicted with acceptable models (grey-shaded area), and (right-bottom) synthetic ZH ratios with noise added (blue dots) and their uncertainties (error bars), predicted ZH ratio curve based on the input (green line) and inverted model (red line), and threefold s.d. interval of ZH ratios predicted with the acceptable models (grey-shaded area).

Figure S4. Tests of dependency of joint inversion results on the initial models used in the inversion, based on the inverted seismic model of a basin station (T0.K001). Left-hand panel: the initial (solid lines) and inverted models (dashed lines) in two inversion tests (models colour-coded with each test), along with the input model (blue line). Right top-hand panel: synthetic RF computed with the input model (blue line) and its uncertainty (two black lines), along with RFs predicted with the inverted models of the two tests (dashed lines, colour-coded with the tests). Right bottom-hand panel: synthetic ZH ratios (black dots) predicted with the input model and their uncertainties (error bars), along with synthetic ZH ratio curves computed with the inverted models of the tests (dashed lines, colour-coded with the tests).

Figure S5. Left-hand panel: thicknesses and (right-hand panel) their associated uncertainties of (a, b) the soft sediment, (c, d) compacted sediment and (e, f) total sediment, obtained in individual stations and plotted on the sites of seismic stations in the continental China. Block boundaries, plate boundaries and major sedimentary basins are marked by solid grey lines, blue lines and grey-crossed areas, respectively.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

We are grateful to the colleagues in the China Earthquake Administration and other related institutions on the deployment and operation of the seismic stations, and the China Seismological Reference Model Project for its efforts in collecting and sharing the seismic data (http://chinageorefmodel.org/). We thank editor Lapo Boschi, Walter Mooney and an anonymous reviewer for their constructive reviews that have significantly improved the paper. This work is supported by the Key Research Program of Frontier Sciences of Chinese Academy of Sciences (QYZDY-SSW-DQC020, CAS) and the National Natural Science Foundation of China under grants NSFC41874002 and NSFC41674045. We also thank the use of GMT (Wessel et al. 2013), SOD (Owens et al. 2004), TauP (Crotwell et al. 1999) and ObsPy (Beyreuther et al. 2010) software packages.

DATA AVAILABILITY

Seismic waveforms used in this study consist of three parts: (1) CNSN data is provided by the China Earthquake Networks Center and AH, BJ, BU, CQ, FJ, GD, GS, GX, GZ, HA, HB, HE, HI, HL, HN, JL, JS, JX, LN, NM, NX, QH, SC, SD, SH, SN, SX, TJ, XJ, XZ, YN, ZJ Seismic Networks, China Earthquake Administration. (2) ChinaArray data is provided by the China Seismic Array Data Management Center at Institute of Geophysics, China Earthquake Administration (doi: http://dx.doi.org/10.12001/ChinArray.Data). These waveforms are accessible from these data centres following their respective regulations. (3) IRIS data are openly accessible from the Incorporate Research Institutions for Seismology (IRIS) Data Services (ds.iris.edu, last accessed December 2018).

The estimated P-wave polarization, isotropic ZH ratio dispersion curves and receiver functions in this paper are available in the level 2 data repository of the China Seismological Reference Model Project at http://chinageorefmodel.org/en/data/level2, and can be accessed following their regulation. The constructed 3-D Vs model, sediment thickness maps and near-surface Vs model in this article are available in the member-contributed model repository of the China Seismological Reference Model Project at http://chinageorefmodel.org/en/china_models_individual/crust_models.