The lithospheric thinning of the North China Craton inferred from Rayleigh waves inversion

G J I T ec t o n i c s a nd g e o d y n a m i c s Joint inversion of interstation phase and group velocity dispersion data of Rayleigh waves provides the regional shear-wave velocity structure of the lithosphere across the central and eastern block of the North China Craton (NCC). Both the velocity dispersion curves and V s structure show strong lateral variation across the NCC, which correlates well with the regional tectonic evolution. At the lower crustal level, the characteristics of V s structure reveal differences in composition of the two blocks. In the upper mantle, the lithospheric removal in the eastern NCC is more severe than in the central NCC. The V s differences between the central NCC and the eastern NCC imply that multistage reactivation of the NCC occurred in the Mesozoic and Cenozoic time.


I N T RO D U C T I O N
The North China Craton (NCC) is one of the world's oldest Archean cratons. The NCC is divided into the Western Block, the Eastern Block and Trans-North China Orogen/Central Block based on the age, lithological assemblage, tectonic evolution and P-T-t paths (Fig. 1). The central block is a collisional orogen between the Western and Eastern NCC that formed during the amalgamtion of the NCC. As such, it consists of a number of metamorphic terranes with different petrological compositions (e.g. Zhao et al. 2001). At this time, it is a zone of changes in crustal thickness (from 37 to 42 km, Li et al. 2006), Boguer anomaly (from −25 to −100 mGal, Zuo et al. 2002) and lithospheric thickness (Chen et al. 2008). The eastern boundary of the central block is near the North-South Gravity Lineament (NSGL) (Zhao et al. 2001;Zuo et al. 2002).
A wealth of observations suggests that more than 100 km cratonic mantle lithosphere has been lost, in part or in whole, beneath the NCC during the Mesozoic and Cenozoic era. The mechanisms for such a change beneath the NCC is hotly debated (e.g. Menzies et al. 1993;Deng et al. 1994Deng et al. , 1996Menzies & Xu 1998;Griffin et al. 1998;Zheng et al. 2001;Huang et al. 2003). The contrast in crustal and lithospheric thickness, topographic height, gravity anomaly, shear-wave splitting measurements, lithological and geochronological composition across the NSGL suggests that the NSGL is the westernmost demarkation of Mesozoic-Cenozoic thinning (e.g. Griffin et al. 1998;Menzies & Xu 1998;Zheng et al. 2001;Menzies et al. 2007;Zhao et al. 2007). However, temporal and spatial evolution of Mesozoic and Cenozoic volcanism further implies that the central block also has been involved in the lithospheric removal. The timing responsible for lithospheric removal beneath the NCC still remains controversial (Deng et al. 1994(Deng et al. , 1996(Deng et al. , 2004Zhou et al. 2002;Gao et al. 2004;Xu et al. 2004;Qiu et al. 2005;Xu 2006Xu , 2007. Clearly, detailed knowledge of the lithospheric structure is quite important for understanding the tectonic processes of the NCC. Global seismic tomographic inversions reveal that the lithosphere beneath the NCC is significantly thinner than the lithosphere beneath other cratonic blocks (e.g. Ritzwoller et al. 2002;Lebedev & Nolet 2003;Lebedev & Van der Hilst 2008). However, due to low spatial resolution, the tomographic observations do not provide sufficient constraints on the local structural complexities in this region. Previous seismic studies on this region have indicated a complex and dramatically thinned lithosphere beneath the NCC (e.g. Huang et al. 2003;Hearn et al. 2004, Menzies et al. 2007Xu 2007;Chen et al. 2008). However, they mainly deal with the pronounced differences in the lithospheric structure between the eastern and Western NCC, with little emphasis on the central NCC.
A passive seismic array experiment across the central and eastern block of the NCC is operated from October 2006 to 2009. The network consists of 250 seismic stations, among which 200 are equipped with broad-band sensors (10 CMG-3T with a flat velocity response between 50 Hz and 120 s, and 190 CMG-3ESP with a flat velocity response between 50 Hz and 60 s; see Fig. 1). Fig. 1 shows the main tectonic units of the NCC and 22 broad-band seismic stations used in this study. The dense distribution of these broadband stations in this area provides a unique opportunity to obtain well-constrained local average shear-wave velocity models. Using  Zhao et al. 2001). NSGL indicates the North-South Gravity Lineament (Zuo et al. 2002). The paths along which we measured the interstation phase and group velocities can be found in Table 1. Notice that each set of paths delineates a region corresponding to a single morphotectonic unit. teleseismic surface waves recorded during the experiment, we determine the lateral variations of the shear-wave velocity structure of the lithosphere, through the joint inversion of the interstation group and phase velocity dispersion data of the fundamental mode Rayleigh wave.

DATA S E L E C T I O N A N D D I S P E R S I O N M E A S U R E M E N T S
The methodology used here relies on the measurement of the phase and group variation between two stations for the fundamental mode Rayleigh waves. To ensure a good signal-to-noise ratio, only events with a magnitude greater than 5.0 are considered. The minimum magnitude of the selected events is 5.6 in this study (Table 1). As the energy of higher modes usually increases with the earthquake depth, we only select teleseismic events with depth smaller than 120 km in the USGS NEIC catalogue (80 per cent of the events less than 50 km). Events too close to the array are typically small and poorly dispersed, whereas very distant events exhibit multipathing effects (Lebedev et al. 2006;Laske & Weber 2008). Therefore, events too close (<20 • ) or too distant (>100 • ) from the array are not selected in our study. The two-station great-circle method for surface wave dispersion analysis is based on the principle that the earthquake lies along a common great-circle path joining the two stations. In addition, to improve in the reliability of our Rayleigh dispersion measurement, we must select a large number of earthquakes with epicentres on the interstation great circle path. However, this is usually difficult with temporary installations. The only solution with a limited data set is to consider events located slightly away from the interstation great circle path. In this study, Rayleigh wave phase and group velocities are measured between various station pairs for the selected 29 events with propagation paths deviated by 10 • or less from the great-circle arcs connecting the station pairs and the epicentre (75 per cent of the propagation paths deviated less than 5 • ) ( Table 1). As expected from the distribution of the events, the coverage is excellent with many crossing rays in both the central and eastern block of the NCC (Fig. 1). The density of crossing rays decreases with increasing period because fewer earthquakes generate surface waves with good signal-to-noise ratio at long periods. In addition, the sparsity of data at longer periods also is related to the band limitation of the sensor (i.e. 60 s eigenperiod of the CMG-3ESP sensor).
We use the vertical-component Rayleigh wave seismograms because they are not contaminated by the interference of Love waves and typically have lower noise levels than the horizontal components. Interstation phase velocities of the fundamental mode Rayleigh wave are then measured with the multiple filter analysis method of Herrmann & Ammon (2004). Dispersion data at the defined periods with 2 s interval are interpolated from phase velocities at observed frequencies by a spline interpolation method. We also measure Interstation group velocity dispersion data of the fundamental mode Rayleigh wave by cross-correlation and multiple filtering analysis (MFA) (Dziewonski & Hales 1972;Herrmann & Ammon 2004), whereas some smoothing is implicit in the MFA method.
Interstation paths with similar velocity dispersion curves are shown in Table 1 and Fig. 1. Nineteen and eighteen acceptable paths traversing the eastern and central block form the final data set, respectively. Some interstation paths from the same event are collinear and some for the same two-station pair using different events. It is regrettable that we did not obtain enough Love wave dispersion data to test the presence of radial anisotropy. The velocity models derived in the following should thus be strictly treated as representative of the velocity of subvertically polarized shear waves.
As it is known that the structure beneath the central and eastern block of the NCC shows clear lateral variations (e.g. Zhao et al. 2001;Qiu et al. 2005;Li et al. 2006;Chen et al. 2008). Group and phase velocities from many events are averaged to reduce the influence of lateral heterogeneity along the path from the epicentres to the stations. So we do not claim that the structures we obtain in the following are representative of any one point within each block, they merely present lateral averages. The assignment of the dispersion data to the two blocks is based on the following two reasons. First, many previous studies have pointed out a change in crust-mantle properties between the two blocks (e.g. Zhao et al. 2001;Zuo et al. 2002;Qiu et al. 2005;Li et al. 2006). Second, the dispersion measurements for different paths within the same block are observed to differ less from each other than from the paths in the other block (Fig. 2).

D I S P E R S I O N R E S U LT S
Individual Rayleigh group and phase dispersion curves and their standard deviation are compiled for all interstation paths within one block. The mean of the Rayleigh wave dispersion curves in each unit along with the standard deviation error bars are shown in Figs 2 and 3. We confirmed that neglecting the ray deflection correction does not lead to erroneous results because most of the events and station pairs are well distributed in both epicentral distances and backazimuths (Table 1). It shows that there are differences in Rayleigh wave dispersion curves for different block (Figs 2 and 3). These lateral variations demonstrate that these tectonic units in the NCC are characterized by different shear-wave velocity structures. The interstation group velocity minimum of Rayleigh wave, which corresponds to the continental Airy phase, occurs at a period of 18 and 22 s in the eastern and central block of the NCC, respectively. The Rayleigh group velocities at short periods (<38 s) are higher by approximately 0.1-0.3 km s −1 in the eastern part than those in the central part of the NCC after the group velocity minimum. Similarly, the phase velocities at periods 12-40 s are higher by approximately 0.07-0.2 km s −1 in the eastern part than in the central NCC. However, the phase and group velocities of the central block lie close to the eastern block curve for period over 40 s. The differences in the phase and group velocity dispersions for the various blocks most likely reflect differences in the seismic characteristics of their crust and upper-mantle structures (down to ∼40-50 km depth). These observations indicate that the crustal and uppermost mantle velocities for the eastern part of the NCC are on average faster than that for the central part of the NCC. At periods shorter than 80 s phase and group velocities are sensitive primarily to V s structure in the upper-mantle lithospheric depth range (depths less than 150 km). The phase and group velocities measured from different blocks of the NCC are almost consistently slower than the global averages (The reference Rayleigh curves were computed for the global reference model AK135 (Kennett et al. 1995)) at periods less than ∼80 s, but are very similar at longer periods. It means that V s in the upper-mantle lithosphere of the NCC is slower than that of the global average and continental craton, in agreement with the tomographic images (e.g. Ritzwoller et al. 2002;Huang et al. 2003;Lebedev & Nolet 2003;Lebedev & Van der Hilst 2008).
It is remarkable that both Rayleigh group and phase velocity dispersion data show almost the same scatter feature at short periods (<10 s) due to strong upper crustal scattering of the surface waves. Our dispersion data are insufficient to resolve the upper crustal structure, but enable us to determine the seismic velocity structure of the mid-lower crust and the uppermost mantle with a substantially higher accuracy compared to the previous surface wave studies (e.g. Ritzwoller et al. 2002;Huang et al. 2003), which are absent of short period dispersion data.

S H E A R WAV E V E L O C I T Y M O D E L S
Rayleigh wave dispersion is primarily sensitive to S-wave velocities, much less to density and P-wave velocities. The sensitivity to Pwave velocities is confined primarily to the crust. Therefore we only invert for S-wave velocities by coupling P-wave velocities to S-wave velocities using a constant Poisson's ratio, which is a reasonable approximation for the crust and uppermost mantle material. On the basis of several geophysical investigations conducted in this region (Ji et al. 2009), P-wave velocity in this study is tied to S-wave velocity using a Poisson's ratio of 0.27. Density in each crustal layer is assumed to follow the Birch equation (Birch 1960).
Rayleigh group and phase velocities are jointly inverted to obtain shear-wave velocity models for the central and eastern block of the NCC using the inversion method of Herrmann & Ammon (2004). We use a differential inversion scheme, which minimizes both the magnitude of the error vector between the observed and computed velocities and differences between adjacent layers, to avoid large velocity changes between adjacent layers. To stabilize the inversion, appropriate damping values were chosen at different stages. A slightly higher damping value of 10 was adopted to avoid an overshoot in the initial two iterations, and then a damping value of 0.1 was chosen in the following iterations.
The inversion of surface waves for seismic velocity structure is a typical non-linear problem, in which an initial model needs to be close to the true model to obtain a reliable velocity model. Seismic tomography can provide constraints on the structure of the lithosphere beneath the study region (e.g. Lebedev et al. 2006. The reference models in the initial inversions comprise the structure in the upper mantle from the global 1D reference model AK135 (Kennett et al. 1995) and structure in the crust from the 3D global shear-velocity models by  (http://ciei.colorado.edu/ ∼nshapiro/MODEL/index.html). Layer thickness is set to 10 km below the Moho interface.
The inversion procedure was stopped when the theoretical dispersion curve falls within the error bars of the observed dispersion curve (standard error < 0.1 km s −1 ). To account for the nonuniqueness of inversion, we compute 100 average dispersion curves by randomly excluding ∼1/3 observed curves. Then we use these perturbed dispersions to invert for V s model with the same starting model. Due to the non-uniqueness of the problem, each perturbed dispersion data set results in a slightly different shear velocity structure. If the predicted dispersion for a model is within the error bars of the observations, then this model is retained. All models that fall within the error bars are averaged to produce the final model. We checked that these average models are also solutions to the observed dispersion data within the error bars. Thus these average models are assumed to be the most representative models among our solutions. Their standard deviation is 0.04-0.13 km s −1 depending on depth (Fig. 4).
As the fundamental mode Rayleigh wave phase and group velocities are mainly sensitive to the shear-wave velocity structure around the depth of one-third to half of the wavelength, and the sensitivity decreases as the period increases, the obtained 1-D velocity models are limited to depth about 180 km, which is smaller than half the wavelength.
Our dispersion data are insufficient to resolve the upper crustal structure. However, the obtained velocity structures in the lower curst and upper mantle are robust and well constrained by the observed dispersion curves (not affected by our choices of parametrization and damping). Notice that the crustal thickness is not constrained by our data but based on Shapiro & Ritzwoller's results (2002).
The lower crustal S-wave velocity of the central block has an average value of 3.5 km s −1 , which are 0.2 km s −1 lower than in eastern block (Figs 4 and 5). Our result is in good agreement with the results from recent body tomographic imaging, which show the lower crust velocity beneath the central block is ∼4-5 per cent slower than in the eastern block (Lei et al. 2008;Tian et al. 2008). The lower crust velocity of the NCC is slower than that of typical cratonic regions (∼3.9-4.1 km s −1 ), but closely matches with that of continental rift areas (3.6-3.8 km s −1 ) (Rudnick & Fountain 1995;Hazler et al. 2001;Lebedev et al. 2006. Our average V s model in the crust differs from the deep seismic sounding (DSS) (Li et al. 2006) and tomographic imaging results (Lei et al. 2008;Tian et al. 2008) in which they found the middle-to-lower crust low-velocity zone (LVZ) beneath the NCC. The existence of such a feature is not precluded by our inversions. The reasons that we can not find the LVZ may be because the depth resolution in our inversion is low or the LVZ is locally distributed.
The V s model shows that the average velocities in the upper mantle are ∼2-3 per cent lower than the AK135 model to at least 120 km depth. Below  and 120 km depth), where V s may be as low as 4.15 km s −1 .
The negative velocity-gradient may indicate the transition from lithosphere to asthenosphere (Figs 4 and 5). Our model largely agrees with a recent surface wave study in China by Huang et al. (2003) who require a LVZ in the upper mantle, similar to the one we found here, to fit their data. The velocity decrease of about 3-7 per cent estimated for the asthenosphere by Chen et al. (2006) is also in very good agreement with the model for the eastern NCC shown here. However, our model is clearly different from the shear-wave models of the NCC by , in which velocity increases almost linearly from the Moho to ∼180 km depth. For depths below ∼120 km, our V s are indistinguishable from the velocity models determined by Shaprio & Ritzwoller (2002).  used the fundamental mode Rayleigh wave group and phase velocities with period up to 16 and 40 s, which have poor constraint on the lower crust and upper mantle. It is possible that the lower velocities at ∼40-70 km depth in the models of  are compensated by their higher velocities in the lowermost crust (Fig. S1, see Supplementary Material). That these features in our V s models are required by the data was verified by a series of tests in which inversions using starting models without explicit Moho interface also show a notable LVZ in the upper mantle (Figs S2 and S3, see Supporting information). Therefore these LVZs are unlikely to be the result of poor assump-tions about the Moho depth. In the central and eastern block, the lower crust velocity has a positive gradient of ∼0.12 km s −1 and ∼0.16 km s −1 per 5 km depth, respectively. It also shows that the lower crust velocity in the central block is slower than in the eastern block.

D I S C U S S I O N
This study yields important new insights into lateral heterogeneity of the V s structure across the NCC. In this section, we only focus our discussion on the main features of those regional average V s models that we consider the most robust.
At the crustal level, we find that the lower crust shear wave velocity is by ∼4-5 per cent lower in the central block (∼3.5 km s −1 ) than in the eastern block (∼3.7 km s −1 ). The temperature difference between two blocks probably accounts for less than 1 per cent in this velocity difference (e.g. Wang 2001). The result of a low V s lower crust in the central block agrees well with the DSS V p measurements in this region that are considered to be due to a predominantly felsic crustal composition (Li et al. 2006) in a medium geotherm, but not quite reaching melting conditions (Wang 2001). The result of relatively higher dispersion velocity and V s we find in the eastern block is also similar to the V p obtained from the DSS observations, which implies that its lower crust could have an intermediate to mafic composition (Li et al. 2006) with a slightly higher cratonic geotherm inferred from the Neogene volcanism and heat flow analysis (Griffin et al. 1998;Wang 2001). Based on the geological and geophysical investigation, Qiu et al. (2005) suggest that the continental crust of the central and eastern NCC with an original crustal TTG (trondhjemite-tonalite-granodiorite) component was reconstructed to be granitic crust during the Yanshanian orogenic movement. The continental granitic crust of the eastern NCC was basified again by the eruption of basalt magma along the continental rifting in the Cenozoic Period.
Our V s model shows that the LVZ extends from 65 to 120 km with the lowest velocity of ∼4.15 km s −1 at a depth of 80-90 km. Clearly, the range of LVZ in the upper mantle for the central block (70-130 km) is deeper than that for the eastern block (∼60-110 km). This LVZ underlies a relatively high velocity upper-mantle lid. The velocity contrast between the LVZ and the upper-mantle lid is about 3 per cent. In surface-wave tomography, the base of the lithosphere is sometimes inferred from the strongest negative gradient in Swave velocity marking the base of a high velocity lid (e.g. Debayle & Kennett 2000). If so, our best estimate of the average lithospheric thickness in the central and eastern block of the NCC is about ∼80 km and ∼60 km, respectively. Some other studies also observed similar thickness of the lithosphere in this area. For example, the results of receiver function migration showed that the lithospheric thickness in the central NCC is about 80-100 km (Chen et al. 2008;Wu et al. manuscript in preparation, 2008), whereas the lithospheric thickness in the eastern NCC is about 60-70 km Sodoudi et al. 2006;Wu et al. manuscript in preparation, 2008). Compositional data from the younger mantle samples reveal that the lithosphere-asthenosphere transition beneath the central and eastern NCC is not more than 100 km in depth (e.g. Menzies et al. 1993;Deng et al. 1994;Fan et al. 2000;Shi et al. 2000).
It is obvious that the lithospheric thickness beneath the central and eastern NCC is thinner than that beneath the adjacent Ordos Platform (>200 km) (e.g. Menzies et al. 1993;Griffin et al. 1998;Huang et al. 2003;Deng et al. 1994Deng et al. , 1996Deng et al. , 2004. The combina-tion of thin lithosphere and average shear velocity of only about 4.35 km s −1 in the relatively high velocity mantle lid indicate a partial loss of mantle lithosphere, accompanied by massive rising of asthenospheric mantle beneath the central and eastern NCC. Geological evidences suggest that a thick cratonic lithospheric root existed beneath this region in the Paleozoic. Since the late Mesozoic, the lithosphere has been thinning from the normal cratonic root to the present-day thickness in this region (Menzies et al. 1993(Menzies et al. , 2007Griffin et al. 1998;Deng et al. 1994Deng et al. , 1996Deng et al. , 2004. The direct cause of lithospheric removal in the late Mesozoic is controversial, and various mechanisms have been postulated, including thermal erosion or lithospheric spreading (Xu 2001;Zheng et al. 2001) and delamination (Zhou et al. 2002;Gao et al. 2004). Some authors (Deng et al. 1996(Deng et al. , 2004Qiu et al. 2005) suggest that the lithospheric removal was caused by the delamination beneath the central NCC during the Jura-Cretaceous time, whereas the eastern block of the NCC has been subjected to not only this process, but also initiative rifting and thermal erosion. The latter resulted in the thinning of cratonic lithospheric root in Cenozoic time again. This interpretation is reasonable for our observations beneath the central and eastern NCC.

C O N C L U S I O N
Interstation phase and group velocity dispersion curves of the fundamental mode Rayleigh wave are measured from all the available broad-band teleseismic records by the NCC Array and jointly inverted for the average 1-D V s models for the eastern and central block of NCC. Both the dispersion curves and the resulting velocity models indicate that the nature of the crust and upper mantle changes significantly from the eastern block to the central block across the NCC. The lower crustal S-wave velocity of the central block are 0.2 km s −1 lower than that in eastern block, which can be explained by the differences of their petrologic composition. The average 1-D V s models also reveal pronounced LVZs at the depth range from ∼65 to ∼120 km underneath a lithospheric lid. The average S-wave velocities in the mantle lid are 2-3 per cent slower than the AK135 reference model. If we use the center of the negative velocity gradient to infer the base of the lithosphere, the average lithospheric thickness in the central block is ∼20 km thicker than that in the eastern block, suggesting a more severe lithospheric removal in the eastern part of the NCC than that in the central block of the NCC. The V s differences between the central NCC and the eastern NCC may be related to multistage reactivation of the NCC occurring in the Mesozoic and Cenozoic time.

A C K N O W L E D G M E N T S
We acknowledge the participants of the Seismic Array of IGCEA for collecting the data. We are grateful to Robert B. Herrmann and Charles J. Ammon for providing their Computer Programs in Seismology. We thank Professor Barbara Romanowicz, He Zhengqin and Xu Zhonghuai for comments on an earlier draft. We greatly appreciate many constructive comments by Yao Huajian, Miss Sylvia Hales, an anonymous reviewer and the Editor, Prof. Jeannot Trampert for improving the manuscript. This work is supported by the NSF of China

S U P P O RT I N G I N F O R M AT I O N
Additional Supporting Information may be found in the online version of the article: Figure S1. A comparison of V s and synthetic dispersions in the eastern block of NCC. It shows that the synthetic dispersion computed from model of Shaprio & Ritzwoller (2002) is not consistent to our average dispersions, especially at short period.  used the fundamental mode Rayleigh wave group and phase velocities with period up to 16 and 40 s, which have poor constraint on the lower crust and upper mantle. It is possible that the lower velocities at ∼40-70 km depth in the models of  are compensated by their higher velocities in the lowermost crust. Our tests also show that it happens in the central block. Figure S2. Effect of the initial crustal structure on the velocity model. The inversions with starting models (grey dot lines) without explicit Moho show a notable LVZ in the upper mantle of the resulting model (grey solid lines) .It also show that the resulting model is similar to our final V s model (black solid lines) below ∼60 km depth. Figure S3. Effect of the initial model on the velocity model. In the central and eastern block, the lower crust velocity has a positive gradient of ∼0.12 and ∼0.16 km s −1 per 5 km depth, respectively. It also shows that the lower crust velocity in the central block is slower than in the eastern block.
Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.