Seismic evidence for the stratified lithosphere in the south of the North China Craton

The North China Craton (NCC) formed in the Paleoproterozoic and suffered cratonic destruction in the Mesozoic. It consists of the Western NCC (WNCC), the Trans‐North China Orogen (TNCO) and the Eastern NCC (ENCC). We investigated the upper mantle structures in the southern NCC by using receiver functions analysis. Polarization analysis was applied to increase the quality of receiver functions. The seismic images show significant stratified structures in the upper mantle in the south of the NCC. In the south of the WNCC, the lithosphere is ~280 km thick and contains a low‐velocity zone at ~110–220 km depth. This low‐velocity zone is interpreted as an intracratonic partial melting zone. In the south of the TNCO and ENCC, there is a low‐velocity zone at ~60–200 km depth, which results from the partial melting during the Mesozoic cratonic destruction. Its upper boundary is interpreted as new Lithosphere‐Asthenosphere Boundary at ~60–100 km depth. Below this low‐velocity zone, there are the remains of Archean lithosphere (RAL). Blocks of RAL stagnate in the asthenosphere and reach depths over 300 km. A high‐velocity slab extends from the uppermost mantle of the WNCC into the RAL with a dip‐angle of ~20°. This slab is interpreted as the Paleoproterozoic slab, which indicates an eastward subduction between the WNCC and the ENCC. The WNCC is stable. The lithosphere in the south of the TNCO and ENCC suffered the Mesozoic cratonic destruction, while the destruction is relatively weaker than that in the north of the ENCC.


Introduction
[2] Cratons are characterized by tectonic inactivity and stable lithospheres over billions of years. Recent studies bring new understanding about cratons. The Proterozoic and Archean tectonics of cratons was revealed by subductionrelated fossil structures [Bostock, 1998;Calvert, 1995;Cook et al., 1998]. Cratons cannot be as stable as previously thought, and the mantle root of many cratons has been removed or modified [Griffin et al., 1998;Lee et al., 2001]. The mantle lithosphere of cratons may be stratified [Moorkamp et al., 2010;Yuan and Romanowicz, 2010]. The North China Craton (NCC), the oldest and largest craton in China (Figure 1), is perhaps the most suitable example for studying the characters of cratons.
[3] The NCC consists of the Eastern NCC (ENCC), the Western NCC (WNCC) and the Trans-North China Orogen (TNCO). The NCC formed by the assembly of the ENCC and the WNCC at~1.8 Ga, and the related collision belt is the TNCO [Zhao et al., 2005]. The direction of the related subduction is still under debate [Faure et al., 2007;Kröner et al., 2005;Trap et al., 2009;Zhao et al., 2005;Zheng et al., 2009].
[4] The cratonic lithosphere of the ENCC used to be very thick [Gao et al., 2002;Xu, 2001]. However, its root is onlỹ 60-80 km thick at present, because most of the root was removed during the Mesozoic cratonic destruction [Chen et al., 2009;Wu et al., 2005;. The magmatism related to cratonic destruction lies mostly in the north and east of the NCC, while it is rare in the south of the ENCC (Figure 1b)  . Moreover, the seismic tomography shows significant longitudinal variation in the upper mantle beneath the ENCC [Tian et al., 2009;Xu and Zhao, 2009]. A high-velocity anomaly beneath the ENCC varies from~100 km thick in the north to~400 km thick in the south. We infer that the cratonic destruction is of relatively low degree in the south of the ENCC, and the evidence of fossil structures could be preserved.
[5] The WNCC is regarded as stable, while its lithospheric thickness is still under debate. It is over 200 km thick in the north of the WNCC from Receiver Functions (RF) analysis [Chen, 2010;Chen et al., 2009], and~150 km from surface wave tomography studies [Feng and An, 2010;Huang et al., 2009]. Additionally, a mine blasting seismic study reveals two velocity boundaries in the upper mantle beneath the southeastern WNCC, at~110-120 km depth and~200-220 km depth, respectively [Zhao et al., 2008]. These results indicate complex deep structures beneath the WNCC.
[6] The previous seismic studies mainly focused on the northern and central ENCC, while a few studies were carried out in the south of the ENCC. Thus, in this study, we deployed portable broadband seismic array in the south of the NCC (along Lat. 36 N; Figures 1 and 5a) and investigated the lithospheric structures by using RF analysis. The seismic profile crosses the western boundary of the WNCC and reaches the western margin of the ENCC. This profile almost overlaps a crustal section published recently [Wei et al., 2011;Zheng et al., 2012]. Here we focused on the refined structures of the lithosphere and its lateral variation, which provide evidence for understanding the structures and evolution of the NCC.

Methods and Data
[7] The RF technique uses converted seismic waves to detect velocity discontinuities beneath seismic stations [Langston, 1979]. The compression-to-shear wave RF (PRF) and shear-to-compression wave RF (SRF) are widely used to image the crustal and upper mantle structures [Farra and Vinnik, 2000;Kumar et al., 2006;Niu et al., 2005;Yuan et al., 2006]. The SRF is especially efficient to detect the discontinuities at 50-250 km depth, for which strong signal-generated interfering phases are present in the PRF images. However, the S wave usually contains several cycles with the period of 3-5 s or even longer, and the multiples SpPp and SsSp probably arrive close to or before the end of S wave. Therefore, to ensure the quality of source function for deconvolution, it is necessary to increase the precision of the data cutting. Another disadvantage of SRF is that, for the same depth and the same event, the distance from the converted point to the receiver of SRF is much farther than that of PRF. Thus, some converted points may be far away from the target region, and will interfere with seismic images, if the deep structures are laterally heterogeneous. It is necessary to consider the structural variation when applying Common Conversion Point (CCP) stacking.

Polarization Analysis
[8] Polarization analysis is a general data processing for decomposing multiple axes seismograms into the polarization directions. In this study, we applied polarization analysis to decompose seismograms into the P, SV, and SH systems before data cutting. The S wave can be distinguished from its multiples SsSp and SpPp by different polarization patterns. This processing increases the precision of the data cutting and obtains high quality source functions. As the seismograms are decomposed into the P and SV components, the signal-to-noise (S/N) ratio increases as well [Kennett and Engdahl, 1991;Kumar et al., 2006Kumar et al., , 2007.
[9] According to the work of Kanasewich [1981] and Jackson et al. [1991], time domain polarization analysis decomposes seismograms by the statistical polarization characteristics in a time window. For each sampling time, selecting certain time window and computing the covariance matrix; calculating the eigenvalues and the eigenvectors of the covariance matrix; figuring out the polarization directions; decomposing seismograms into the P, SH, and SV components. An example of the polarization analysis is given in Figure 2 based on seismic synthetics. In this study, we improved polarization analysis for teleseismic receiver functions.
[10] The converted waves only exist in the P and SV components. Thus, we only decompose the Z and R components into the P and SV components. The time window selection is the key to the success of the polarization analysis, and three requirements were mentioned by Jackson et al. [1991]. First, the window should contain only one arrival; second, the ratio of S/N ratio of the window should be maximized; and third, the window should be as long as possible to allow discrimination of random noise from signals.
[11] It is difficult to satisfy the first requirement in the RF study. While a time window contains several arrivals, the polarizations of which parallel or alternatively perpendicular to each other, the statistical characteristic of this window is elliptical polarization. The axes of this polarization ellipse correspond to the polarization directions. In the RF study, the arrivals are direct waves, converted waves and multiples. Their polarization directions are approximately parallel, or alternatively perpendicular to each other. Therefore, we weakened the first requirement: the window could contain several arrivals, the polarizations of which are parallel or alternatively perpendicular to each other. Based on many tests, the optimal window is about three times of the wavelength of direct P or S.
[12] The modified polarization analysis is as follows: first, the seismograms were band-pass filtered with corner frequencies of 0.03 and 1 Hz, and then rotated into the Z, R, and T components. For each sample t, we give a time window [t -W/2, t + W/2]. Then calculate the covariance matrix V in the window.
where, 1/Δt is the sample ratio, N is the sample number of the time window, and m i is the mean of the window.
[13] The eigenvalues of V are l 1 , l 2 (| l 1 | > | l 2 |). The corresponding eigenvectors of l 1 , l 2 are e 1 , e 2 , respectively (e 1 = [e Z 1 , e R 1 ] T ; e 2 = [e Z 2 , e R 2 ] T ). The P wave vector is e 1 or alternatively e 2 , and the SV wave vector is the other one (Figures 3b-3e). The polarization characteristics could distinguish them, and help decompose the waveform into the P and SV components. The rectilinearity for time t is identified as RL(t) = || l 1 | -| l 2 || / | l 1 |, which indicates elliptical polarization when RL(t) < 1 and linear polarization when RL(t) = 1. There are five different polarization patterns (Figures 3a-3). A function PP(t) is used to represent the polarizing pattern, where PP(t) = sign(e Z 1 • e R 1 ) • RL(t). Essentially, for teleseismic waves, P waves polarize in the first and third quadrants, while S waves polarize in the second and fourth quadrants. Therefore, when PP(t) > 0, e 1 is the P wave vector and e 2 is the SV wave vector; when PP (c) Synthetics of Z, R, and T components generated by using RAYSUM code [Frederiksen and Bostock, 2000]; the pulse is 2 s wide, and the ray comes from the northeast; as the dipping Moho deflects the ray, all phases can be detected in the T component [Langston, 1977]. (d) Waveforms are decomposed into P, SV, and SH components by using polarization analysis according to the works of Kanasewich [1981] and Jackson et al. [1991]; the time window used is 6 s wide.
(t) < 0, e 1 is the SV wave vector and e 2 is the P wave vector. Two operators OP and OS were defined to carry out the signal decomposition: [14] The P and SV components are calculated by [15] A slight modification is required in the practical polarization analysis, because the polarization directions could shift into a wrong quadrant (Figure 3f) for following reasons. First, for the teleseismic events with an epicentral distance over 80 , the polarizations of direct P and direct S would be very close to the Z and the R axes, respectively. Second, sometimes the seismometer might not be correctly oriented due to accidental errors. Thus, there is slight difference between the actual axes of the seismometer (Z' -R', the dashed lines in Figure 3f) and the real axes on ground (Z -R, solid lines in Figure 3f). Third, the ray path could be deflected by complex deep structures [Langston, 1977]. Therefore, we set a tolerance angle a (Figure 3f, E. 8-10) for the polarizations, which are very close to the axes but in the wrong quadrant, decomposed into the right components. The optimal value of a is 5 . Then, the function PP(t) and the operators OP and OS are improved as [16] Figure 4 presents the polarization analysis of a teleseismic event recorded in the western margin of the WNCC. The time window W is 30 s wide. The P and SV components and the function PP(t) are shown in Figure 4b. The PP(t) with a value about -1 indicates the direct S. The slight increasing of the PP(t) indicates the end of the direct S. Contrasted to the seismograms in the Z component, the direct S is suppressed and the multiples are clear in the P component. Figure 4c shows the SRF of this event computed in the Z-R and the P-SV systems.

Receiver Functions
[17] We deployed 42 stations in the south of the NCC (Figure 5a). Each station was equipped with a GMT-40 three-component seismometer, a Reftek 72A/130 type digital acquisition, and a GPS timing system. These data were recorded at many stations over a period of several years with most of the individual stations operating about one year. We visually examined the teleseismic events with magnitude greater than 5.0 and epicentral distances in the range of 30 -110 , and choose the events with high S/N ratio. Among these high quality events, we choose P wave with the epicentral distance of 30 -90 for PRF, S wave of 60 -90 and SKS wave of 90 -100 for SRF. For both PRF and SRF, the polarization analysis was applied before the deconvolution. We employed maximum entropy deconvolution [Wu et al., 2003] to estimate RF. Auto-correlation and cross-correlation are determined by the maximum entropy method. Then the Toeplitz equation and the Levinson algorithm are used to calculate the iterative formula of error-predicting filter. RFs are estimated with a 1.5 Hz Gaussian filter. Finally, we obtained 2989 PRFs and 628 SRFs after quality control.
[18] The S/N ratios for the PRFs in the Shanxi rift are very low (Figures 6a-6c). The PRF studies of Langston and Hammer [2001] revealed that the complex structures beneath sedimentary basins may induce strong P wave coda and interfere with the deconvolution. The strong multiples (especially PpPp) generated by dipping interfaces would (f) Z -R is the coordinate on ground; Z' -R' is the actual coordinate of the seismometer; tolerance angle a is introduced to make the polarizations, which shift into wrong quadrant, decomposed into right components (E. 8-10).
result in the failure of the deconvolution as well . The Shanxi rift is covered by thick sediments, and the complex structures below cause the failure of calculating PRF (Figure 6b). Therefore, we excluded the PRFs in the Shanxi rift while applying CCP stacking.
[19] In contrast with the PRF, the quality of the SRF in the Shanxi rift is almost as good as that outside the rift (Figures 6d-6f). We infer that the complex structures   2004-2005; yellow, 2005-2006; green, 2006-2007; red, 2007-2008; white, 2008-2009; the purple and the blue spots are the predicted Sp piercing points at 120 and 250 km depth, respectively; black rectangles are the sliding windows for CCP stacking. (b) Sketch of CCP stacking beneath a single rectangle. (c) Epicenters of the high quality RF; blue circles for PRF, red crosses for SRF. influencing the P wave coda are probably not sensitive to S wave, because S wave attenuates significantly in sediments and is relatively lower in frequency contrasted to P wave.

Sliding Window CCP Stacking
[20] The projections of Sp piercing points at the depth of 120 km (purple) and 250 km (blue) are computed based on the IASP91 model ( Figure 5a). They cover a very large area. We only stack the CCP in the target region to reduce the interference of the lateral variations of deep structures. However, this reduces the data coverage of CCP stacking. Another problem is that our latitudinal profile does not fit the surface geological strike. The geological strike is NE-SW in the TNCO and ENCC, while it is NNW-SSE in the western margin of WNCC. In the general CCP stacking, the length of the stacking windows (or bins) is fixed to perpendicular the profile. The results will be better, if the lengths of the windows are flexible to fit the geological strike.
[21] Therefore, to ensure the data coverage and make the CCP stacking consistent with the geological strikes, we introduced a sliding window method for the CCP stacking. Based on the IASP91 model, we calculated the predicted piercing points from the surface to 450 km depth with an interval of 0.2 km for all the RFs (Figure 5a). Then, a rectangle (sliding window) in the target region is set for CCP stacking, the length of which parallels the strike of the surface structures ( Figure 5b). The size of the rectangle on ground is 150 km by 30 km, and it could be even larger to cover sufficient piercing points or smaller to increase the lateral resolution. As the rays become sparse with increasing depth, we set the width and length of the rectangle to increase with depth ( Figure 5b). In this study, it is about 150 km by 30 km on ground and 200 km by 40 km at 450 km depth.
[22] For the predicted piercing points at depth, only those within the rectangle are stacked and averaged (Figure 5b). Repeating this for each depth, we obtain a curve reflecting the velocity discontinuities beneath the rectangle. The CCP image was obtained by shifting the rectangle along the profile. The shifting step was~20 km, which was about two thirds of the width of the rectangle. In Figure 5a, the piercing points at 120 and 250 km depth for SRF are plotted. These indicate good data coverage from top to bottom. The epicentral distribution of the high quality RFs shows good distance and azimuth coverage (Figure 5c). The PRF and SRF CCP images are presented in Figures 7a and 7b. The RF number at special depths is shown in Figures 7c-7h.

CCP Stacking Results
[23] The CCP profile is~1000 km long, which crosses the WNCC and reaches the western margin of the ENCC (Figures 1 and 7). This profile images the velocity discontinuities from the surface to~450 km depth. There are many traceable signals from the surface to the bottom. We present the stacked RFs at the 150, 250, and 350 km depths to demonstrate the reliability where the seismic rays become sparse. The data coverage of the PRF image (Figures 7c-7e) is much better than that of the SRF image (Figures 7f-7h). Though the data coverage of the SRF image is not as good as expected, it shows traceable signals in the upper mantle (Figure 7b). To insure the reliability of these signals, we present the SRF moveout corrections without stacking in Figure 8.
[24] All the SRF are moveout corrected to the case of ray-parameter p equaling zero. The SRF piercing the stacking window at the 250 and 350 km depths are selected, then are stored by the piercing points' projections on the profile. The moveout corrections results are shown in Figures 8b and 8c.
[25] The lithosphere of cratons is characterized by highvelocity anomalies in seismic tomography. The Lithosphere-Asthenosphere Boundary (LAB) of cratons is over 200 km depth, and presents negative velocity contrast. LAB could be identified by the negative signals of RF. Here in Figure 8a, the blue lines N1-N4 are the traceable negative signals on the SRF CCP image (Figure 7b). The N2 and N4 are divided in to two parts. For the moveout corrections (Figures 8b and 8c), only the negative signals are plot in black to detect the LAB.
[26] In Figures 8b and 8c, the original signals beneath the WNCC align and indicate high reliability of the traces N2w and N3. Beneath the TNCO and ENCC, the signals are more complicate than those beneath the WNCC. The original signals for traces N2e and N4, which are traceable after stacking, do not align as well as those for N2e and N3. We thought the velocity discontinuities N2e and N4 do exist, while the original signals could not well align due to significant lateral velocity variation in the upper mantle. The lateral variation from crust would introduce several kilometers offset, while this offset could not be distinguished due to the 1.5 Hz Gaussian lowpass. And the significant lateral velocity variation in upper mantle is detected by seismic tomography [Tian et al., 2009;Xu and Zhao, 2009]. Thus, beneath the TNCO and ENCC, the anomalies in the upper mantle mostly interfere with the alignment of the original signals. In the SRF CCP image, the traces N2e and N4 exist, but with some ambiguity in the depth.

Interpretations and Discussions
[27] Our interpretation of CCP stacking images (Figures 7a and 7b) is shown in Figures 9a and 9b. Red lines and blue lines are the traceable positive and negative signals, respectively. The Moho discontinuity is clear in the PRF CCP image at the depth of~30-60 km ( Figure 9a). As the total SRFs are fewer than the PRFs, the Moho converted signals in the SRF CCP image (Figure 9b) are not as distinct as those in the PRF. Therefore, we present an integrated interpretation of the Moho based on PRF and SRF. The dotted Moho means poor data coverage. There is another latitudinal profile along Lat. 36 N [Wei et al., 2011;Zheng et al., 2012]. Our results are consistent with that profile, which indicate the Moho depth is~60 km beneath Liupan Orogen,~40 km in Ordos Block, and~30 in the ENCC.
[28] In the SRF CCP image (Figure 9b), the positive and negative traces beneath the Moho are marked as P1-P5 (red lines; velocity increasing) and N1-N5 (blue lines; velocity decreasing), respectively. Because there are no multiples in SRF, these discontinuities indicate the velocity discontinuities in the upper mantle. The velocity discontinuities in the uppermost mantle and Moho generate multiples in PRF, which correspond to the signals at~50-250 km depth after time-depth migration. These multiples make it impossible to distinguish the signals for the real velocity discontinuities at~50-250 km depth in the PRF CCP image (Figure 9a). Only the signals with the depth over 250 km, such as N4e, P5, N5, and part of the N2, could be traced in both PRF and SRF CCP images. There are some short and weak traces with thin dashed lines beneath the ENCC and TNCO, which are lower in reliability than the N1-5 and P1-5. The velocity discontinuities indicate complex and stratified upper mantle in the south of the NCC. Thus, together with previous seismic studies, our results provide an opportunity to identify the refined cratonic structures in the south of the NCC.

The Lithospheric Structures Beneath the WNCC
[29] Beneath the WNCC, tomography results revealed a high-velocity anomaly (H5, Figure 9c) extending from the surface to~400 km depth in the south, and to the~250 km depth in the north [Tian et al., 2009]. A recent RF study in the north of the WNCC detected a negative velocity discontinuity at~250 km depth, which corresponds to the bottom of the H5, and the lithosphere is over 250 km thick [Chen et al., 2009]. In the south of the WNCC, another tomography study revealed a high-velocity anomaly thinner than H5, the thickness of which was over 250 km [Xu and Zhao, 2009]. Here we suggest the negative trace N2w at~280 km depth is the LAB in the south of the WNCC (Figures 9b and 9d). The lithosphere in the south of the WNCC is~280 km thick.
[30] The N3 and P3 traces beneath the WNCC (Figure 9b) indicate a negative velocity discontinuity at~110-160 km depth and a positive velocity discontinuity at~180-220 km depth, respectively. The N3 and P3 constrain a low-velocity zone, which is~80 km thick. In the south of the WNCC, the mine blasting seismic study reveals two velocity discontinuities at~110-120 km depth and~200-220 km depth, with low velocities between them [Zhao et al., 2008]. Our results are consistent with the mine blasting seismic study. The low-velocity zone constrained by N3 and P3 also is consistent with the result of the surface wave tomography, which revealed a low-velocity zone at~150-250 km depth beneath the WNCC [Huang et al., 2009]. This low-velocity zone is very similar to the intracratonic low-velocity zones in other stable cratons, which were detected at~120-160 km depth [Moorkamp et al., 2010]. Moreover, the amplitude of the N3 is smaller than the LAB (N2w) below. Therefore, we suggest the N3 and P3 are the intracratonic structures. [31] The velocity discontinuities detected by receiver functions are sharp boundaries, and the negative discontinuity in the upper mantle results from partial melting [Kawakatsu et al., 2009;Yuan et al., 2006]. Thus, we interpret the lowvelocity zone, constrained by N3 and P3, as a partial melting zone in the cratonic mantle lithosphere (Figure 9d). The negative discontinuity N3 is the upper boundary of this partial melting zone, while the positive discontinuity P3 is the lower boundary possibly resulting from the solidifying. This partial melting in the cratonic mantle lithosphere can be explained by the relationship between the mantle melting curve and the geotherm.
[32] The melting curve in the upper mantle strongly depends on the reduction-oxidation (REDOX) conditions [Brey et al., 2008;Dasgupta and Hirschmann, 2006;Foley, 1988Foley, , 2008Green and Falloon, 1998;Taylor and Green, 1988]. The actual melting curve varies between the limit of oxidated conditions [Brey et al., 2008;Dasgupta and Hirschmann, 2006] and the limit of reduced conditions [Green and Falloon, 1998;Taylor and Green, 1988] ( Figure 10a). The deep cratonic lithosphere is dominated by reduced conditions [Foley, 2008], while it becomes more reduced with the depth increasing [Woodland and Koch, 2003]. Thus, the cratonic melting curve is close to the reduced limit. The reduced limit has a significant backbend at~100 km depth. Another backbend is at the LAB over 280 km depth, for the contrast in REDOX status [Foley, 2008]. Finally, we present the inferred cratonic melting curve in Figure 10b. The geotherm exceeds the melting curve around the backbend at~100 km depth and results in the partial melting. As the depth increases, the geotherm is exceeded by the melting curve and results in solidification. This partial melting zone could result in a low-velocity zone in the cratonic mantle lithosphere. Beneath the WNCC, the upper and the lower boundaries of this intracratonic low-velocity zone correspond to velocity discontinuities N3 and P3, respectively.
[33] Intracratonic low-velocity zones have been detected in stable cratons in the previous study. Beneath the Slave Province  craton and the Kaapvaal craton, a low-velocity zone was detected at~120-160 km depth by using joint inversion of receiver functions, surface wave dispersion, and magnetotelluric data [Moorkamp et al., 2010]. Thus, we suggest that the lithosphere of stable cratons is stratified. It probably consists of three layers: an upper unmelted layer, an Intra-Cratonic Partial Melting Zone (ICPMZ), and a lower unmelted layer (Figures 9d and 10b). Obviously, the ICPMZ would contribute a relatively low-viscosity zone in the cratonic lithosphere.
[34] The low-velocity zone beneath the WNCC is~80 km thick with a surface heat flow of 60 AE 8.7 mWm -2 [Hu et al., 2000]. The low-velocity zone beneath the Slave Province Craton is~40 km thick with the surface heat flow of 43 AE 5 mWm -2 [Mareschal et al., 2004;Moorkamp et al., 2010]. An elevated geotherm will enlarge the thickness of the ICPMZ (Figure 10b). Generally, in the south of the WNCC the cratonic lithosphere is~280 km thick with an ICPMZ at~110-220 km depth.

The Lithospheric Structures Beneath the TNCO and the ENCC
[35] In the south of the TNCO and ENCC, the SRF CCP image is very complex (Figure 9b). The negative trace N1 is at~60-100 km depth, and terminates below the western boundary of the TNCO. There are two significant eastward-dipping traces (P2 and N2e). The positive P2 extends from the Moho of the WNCC to~240 km depth beneath the ENCC, while the negative N2e links to the LAB of WNCC. The negative trace N4 varies from~300 km depth beneath the TNCO to~380 km depth beneath the ENCC. These traces indicate significant stratified structures.
[36] The geological studies indicated that the cratonic destruction from late Carboniferous till early Cenozoic resulted in the removal of the cratonic lithosphere of the ENCC, the peak of which was the Mesozoic [Chen et al., 2009;Wu et al., 2005;. However, most of these geological studies were in the north of the NCC, because the related magmatism was mainly observed in the northern and central ENCC but rare in the south of the NCC (Figure 1b). The distribution of the related magmatism indicates that the Mesozoic cratonic destruction spreads from the northeastern edge to the center of the NCC . Seismic tomography provided evidence for the longitudinal variation of the lithospheric structures beneath the ENCC, which is consistent with the geological studies.
[37] In the south of the TNCO and ENCC, the tomography results reflected a low-velocity anomaly (L3) at thẽ 50-200 km depth and a high-velocity anomaly (H1) extending from the surface to~400 km depth (Figure 9c) [Tian et al., 2009;Xu and Zhao, 2009]. The L3 extends northward and increases in volume. In the north of the NCC, it was interpreted as the asthenosphere. The H1 horizontally tapered off northward and was surrounded by L3 in the central NCC. The H1 terminated in the northern ENCC and was interpreted as the remains of Archean lithosphere (RAL), which stagnates in the asthenosphere [Tian et al., 2009;Xu and Zhao, 2009]. In the south of the ENCC, the upper mantle lithosphere has been affected by the cratonic destruction, while blocks of the RAL were reserved. The refined structures of L3 and H1 provide more information about the degree of the Mesozoic cratonic destruction.
[38] In the south of the TNCO and ENCC (Figures 9b and  9c), a low-velocity zone at~60-200 km depth is imaged by the N1, P1, and partial P2. It terminates beneath the western boundary of the TNCO. This low-velocity zone corresponds to the low-velocity anomaly L3. Obviously, it is connected to the asthenosphere beneath the northern ENCC. We suggest that this low-velocity zone results from the Mesozoic cratonic destruction and its partial melting is due to the lateral eroding of the asthenosphere (Figure 9d). The N1 is the upper boundary of this partial melting zone. N1 is consistent with the LAB in the north of the ENCC, which was detected at~60-100 km depth by receiver functions and deep seismic sounding [Chen et al., 2009;Li et al., 2011]. Beneath the ENCC, this negative discontinuity is coherently detected Figure 10. (a) Melting curves for upper mantle peridotite in different REDOX status, modified from Foley [2008]; the limits of oxidized solidus (with H 2 O and CO 2 [Brey et al., 2008;Dasgupta and Hirschmann, 2006]), and reduced solidus (with H 2 O and CH 4 [Green and Falloon, 1998;Taylor and Green, 1988]) are presented. (b) Inferred melting curve (blue solid line), geotherm (red dashed line) of cratons, and an elevated geotherm (red dotted line) are presented; LAB, Lithosphere-Asthenosphere Boundary; ICPMZ, Intra-Cratonic Partial Melting Zone; layered craton is presented.
from the north to the south [Chen et al., 2009;Li et al., 2011;Sodoudi et al., 2006]. Thus, the N1 is probably the new LAB due to the Mesozoic cratonic destruction (Figures 9b and 9d).
[39] The negative discontinuity N4 at~300-380 km depth corresponds to the bottom of H1, which equal to the lower boundary of the RAL (Figures 9b and 9c). Obviously, the RAL contains several traces like P2 and N2e, which indicate significant inner structures.
[40] The positive trace P2 extends from the uppermost mantle of WNCC to the~240 km depth beneath the ENCC. The negative trace N2e links to the LAB of the WNCC and parallels P2 with a distance of~60 km. The traces P2 and N2e image a relatively high-velocity slab extending from the uppermost mantle of the WNCC into the RAL beneath the ENCC. The geological studies suggest that there was an oceanic subduction between the WNCC and the ENCC in Paleoproterozoic, but the subduction direction is still under debate [Faure et al., 2007;Kröner et al., 2005;Trap et al., 2009;Zhao et al., 2005;Zheng et al., 2009]. Here this high-velocity slab imaged by P2 and N2e is interpreted as the slab of the Paleoproterozoic subduction; meanwhile, it indicates that the subduction was eastward with a dip-angle of~20 (Figures 9b and 9d). The thickness of this slab is more or less consistent with the thickness of the subducted oceanic slab [Kawakatsu et al., 2009]. The relatively highvelocity of this slab is consistent with the basalt-eclogite transformation of the subducted oceanic slab, which results in an increase of the velocity and the density [Anderson, 2006;Ruff and Kanamori, 1983;Wang et al., 2005].
[41] Additionally, because the melting temperature of eclogite is lower than peridotite [Jull and Kelemen, 2001;Rapp et al., 1999], in the south of the ENCC, the preservation of the Paleoproterozoic slab indicates that the RAL is stable in thermal condition and inactive in tectonics (Figure 9d). In contrast, the Archean lithosphere beneath the northern ENCC has been removed during the Mesozoic cratonic destruction.
[42] At the bottom of the RAL, there are short discontinuities like N5 and P5 (Figure 9b). The N5 and P5 constrain a low-velocity belt connecting to the asthenosphere directly. This low-velocity belt is probably caused by partial melting as well (Figure 9d), based on the melting pattern at the base of cratonic lithosphere [Foley, 2008]. Here this limited erosion is chemical eroding-REDOX melting, rather than thermal eroding. REDOX melting means that, without temperature increasing, the reduced peridotite of cratonic lithosphere would melt while it is oxidated by the melts from asthenosphere [Foley, 2008]. The intensive tectonic events like oceanic subduction and continental collision would result in incisions and shear zones in the lithosphere. These incisions and shear zones are weak belts in the lithosphere. The REDOX melting would take place along these weak belts, and result in short partial melting belts in the lithosphere. Due to these short partial melting belts and the Paleoproterozoic slab, RAL is divided into blocks.
[43] Generally, the deep structures in the south of the TNCO and ENCC are very complex (Figures 9b and 9d). The new LAB resulting from the Mesozoic destruction ranges from~60 km to~100 km, which terminates beneath the border between the TNCO and WNCC. There is a partial melting zone at~60-200 km depth. Below this partial melting zone, the RAL stagnates in the asthenosphere.
Considering the ambiguity of N4's depth, we thought that the RAL reaches the depth over 300 km. The Paleoproterozoic slab extends from the uppermost mantle of the WNCC into the RAL beneath the WNCC, indicating the subduction between the WNCC and the ENCC was eastward direction with a dip-angle of~20 . Below the slab, there are several short partial melting belts at the bottom of the RAL.
[44] Based on the analysis of the lithospheric structures in the south of the NCC, we present a possible pattern of the Paleoproterozoic subduction ( Figure 11). Before the oceanic subduction (>1.85 Ga), the cratonization of the ENCC should have completed. Its cratonic lithosphere was thick with an ICPMZ. The relatively low viscosity of the ICPMZ made it possible that the slab subducted into the ENCC with a low dip-angle (Figure 11a). The basalt-eclogite transform would make the slab denser than peridotite, so that it subducted into the lower lithospheric layer. Then the Paleoproterozoic slab was stuck in the ENCC because of the relatively high viscosity of the lower lithospheric layer. As the slab subducted into the cratonic lithosphere and the intensive continental collision, the lithosphere of the ENCC thickened. Additionally, the lithosphere would deform and result in shear zones, which are weak belts in the lithosphere (Figure 11b). While the following cratonic destruction took place, the chemical eroding of asthenosphere resulted in the new LAB at~60-100 km depth and the erosion of the weakened belts in the lithosphere (Figure 9d). In the south of the TNCO and ENCC, the Paleoproterozoic slab and blocks of RAL stagnate in the asthenosphere. We suggest that this pattern of subduction that a slab subducted into the ICPMZ might play an important role on the growth and thickening of cratons. Chemical eroding (REDOX melting) plays an important role in the Mesozoic cratonic destruction.

The Degree of the Cratonic Destruction in the South of the NCC
[45] In the south of the NCC, the imaged lithospheric structures (Figure 9d) indicate that the cratonic destruction only occurred beneath the TNCO and ENCC. The WNCC is stable with about 280 km thick cratonic lithosphere. The Figure 11. The inferred pattern of the Paleoproterozoic subduction during the assembly of the NCC. (a) Oceanic subduction; the slab subducted eastward into the ICPMZ of the ENCC with~20 dip-angle. (b) Shear zones in the lithosphere result from the continental subduction; the slab was stuck in the lower lithosphere of the ENCC.
Archean lithosphere in the south of the TNCO and ENCC suffered the Mesozoic cratonic destruction. The new LAB N1 is at~60-100 km depth. Blocks of RAL stagnate in the asthenosphere and reach the depth over 300 km. The Paleoproterozoic slab stuck in the RAL.
[46] Although the entire TNCO and ENCC suffered the Mesozoic cratonic destruction, the degree of the destruction in the south is quite different from that in the north. In the north of the TNCO and ENCC, there are neither RAL nor Paleoproterozoic slab, and the LAB of the ENCC directly connects the LAB of the WNCC by a steep negative discontinuity [Chen et al., 2009]. In contrast, in the south of the ENCC, the LAB is disconnected from the LAB of the WNCC by the Paleoproterozoic slab and the RAL. Combined with the magmatism distribution  and the tomography of the upper mantle [Tian et al., 2009;Xu and Zhao, 2009], we suggest that, beneath the ENCC and TNCO, the Mesozoic cratonic destruction is weaker in the south than that in the north.

Conclusions
[47] The North China Craton consists of the WNCC, the TNCO, and the ENCC. In this study, we carried out seismic investigations in the south of the NCC, and obtained a series of velocity discontinuities from the surface to 450 km depth by using receiver functions technique. We applied polarization analysis before deconvolution and introduced sliding windows to improve the quality of the receiver functions and the common conversion point stacking images. Combining the RF results with previous geological and geophysical evidence, we obtained the stratified lithospheric structures in the south of the NCC.
[48] The lithosphere of the WNCC is~280 km thick in the south, with an intracratonic partial melting zone at 110-220 km depth. In the south of the TNCO and ENCC, there is a partial melting zone at~60-200 km depth resulting from the Mesozoic cratonic destructions. Its upper boundary is the new LAB at the~60-100 km depth. Below this low-velocity zone, blocks of the RAL stagnate in the asthenosphere and reach a depth over 300 km. A high-velocity slab extends from the uppermost mantle of WNCC into the RAL corresponding to the Paleoproterozoic slab. This slab indicates that the Paleoproterozoic subduction between the WNCC and ENCC was eastward direction with a dip-angle of 20 . The stratified lithosphere indicates that the WNCC is stable. In the south of the TNCO and ENCC, the lithosphere suffered the Mesozoic cratonic destruction, while the destruction is weaker than that in the north of the TNCO and ENCC.