Moho Depth Variations From Receiver Function Imaging in the Northeastern North China Craton and Its Tectonic Implications

A detailed knowledge of the crustal thickness in the northeastern North China Craton (NCC) is important for understanding the unusual Phanerozoic destruction of the craton. We achieve this goal by employing a 2‐D wave equation‐based migration method to P receiver functions from 198 broadband seismic stations, using Ps conversions and surface‐reflected multiples. By combining receiver function images along 19 profiles, we constructed a high‐resolution Moho depth model for the northeastern NCC. The results present dominant E‐W Moho depth variations similar to previous observations and new regional N‐S variations beneath both sides of the North‐South Gravity Lineament. To the west, while a deeper Moho (∼42 km) appears in the interior of the Trans‐North China Orogen, a relatively shallow Moho (∼38 km) exists in the northern margin of the Trans‐North China Orogen to western NCC. To the east, the crust beneath the Yan Mountains in the marginal area is thicker (∼32–40 km) than that (∼26–32 km) beneath the Bohai Bay Basin in the craton interior, and the Moho further shallows from NE (∼32 km) to SW (∼26 km) within the basin. Along with other observations, we conclude that the dominant E‐W difference may have been associated with the Paleo‐Pacific plate subduction under eastern Asia since the Mesozoic. The newly observed complex N‐S variations may have reflected the structural heterogeneity of the cratonic lithosphere inherited since the formation of the NCC in the Paleoproterozoic, or spatially uneven effects on the cratonic lithosphere of subsequent thermotectonic events during the long‐term evolution of the craton, or both.


Introduction
Ubiquitous structural heterogeneities in the continental lithosphere are believed to have played essential roles in the tectonic evolution and dynamics of continents, particularly in the stabilization/reactivation of cratons (e.g., Begg et al., 2009;Chen, 2010;Chen et al., 2014). The North China Craton (NCC), one of the oldest cratons in the world, is an ideal place for understanding lithospheric heterogeneities and their influences on the tectonic evolution of continents, because Phanerozoic tectonic activities have led to spatially uneven reactivation and destruction of the cratonic lithosphere and pronounced structural heterogeneities of both lateral and vertical extents in the NCC (e.g., Chen, 2010;Zhu et al., 2011).
The NCC is composed of two Archean blocks, namely, the western NCC and eastern NCC, separated by the Paleoproterozoic Trans-North China Orogen (TNCO, central NCC; Figure 1; Zhao et al., 2005). The whole craton was tectonically stable from its final cratonization ∼1.85 Ga ago (Zhao et al., 2005) to the early Mesozoic, before its collision with the Yangtze Craton to the south (Yin & Nie, 1993;Zhang, 1997). In the late Mesozoic and Cenozoic, the NCC experienced significant tectonothermal reactivation and fundamental destruction of the cratonic lithosphere in its eastern part (Chen, 2010;Fan et al., 2000;Griffin et al., 1998;Menzies et al., 1993Menzies et al., , 2007Wu et al., 2005;Xu, 2001), whereas only localized lithospheric modifications occurred at the rift or marginal areas in its central and western parts Tang et al., 2006Xu, 2007;Xu et al., 2004). In particular, the main portion of the western NCC, that is, the Ordos Block, has remained stable and retained the nature of a typical craton (Kusky et al., 2007;Zhai & Liu, 2003), although lithospheric modifications have been revealed in the circum-Ordos rift systems, that is, the Shaanxi-Shanxi rift to the east and southeast and the Yinchuan-Hetao rift to the northwest (Figure 1b; e.g., Chen, 2010;Xu, 2007). Together with the rapid changes in surface topography and gravity field across the NNE-trending North-South Gravity Lineament (NSGL, Figures 1b and 2), these findings indicate that the different parts of the NCC may have evolved differently in the Phanerozoic (e.g., Chen, 2010;Xu, 2007).
Three major tectonic units in the region, the Yin-Yan orogenic belt in the north, the Bohai Bay Basin in the east, and the Taihang Mountains within the TNCO (Figure 1b), have developed as the direct consequence of the Phanerozoic cratonic rejuvenation. The E-W trending Yin and Yan Mountains were proposed to have undergone multiple phases of the N-S directed contraction in the early-to-middle Jurassic and early Cretaceous, and the regional extension and intense magmatism soon afterward (Davis et al., 2001;Meng, 2003). The Bohai Bay Basin consists of a series of NE-SW oriented uplifts and depressions ( Figure 1b) and is generally believed to have experienced a regional extension with two-phase rifting before the Oligocene and subsequent post-rift phase of thermal subsidence that has lasted until the present day (Allen et al., 1997). The evolution of this basin is thought to have been tectonically coupled with that of the NE-SW trending Taihang Mountains, which probably occurred in conjunction with the NW-SE tectonic extension in the late Mesozoic to Cenozoic (Liu et al., 2000;Ren et al., 2002).
The dynamic causes of the Phanerozoic lithospheric reactivation of the NCC and the variable responses of different parts of the NCC during the process are currently still under debate (e.g., Wu et al., 2008, and references therein). Several mechanisms have been proposed to be responsible for the destruction of the eastern NCC, including the Paleo-Pacific plate subduction (Griffin et al., 1998;Sun et al., 2007), the North-South China collision (Yin & Nie, 1993;Zhang, 1997), the collision of an amalgamated North China-Mongolian plate with the Siberian plate (Davis et al., 2001), the India-Eurasia collision (Menzies et al., 1993), and an enhanced mantle temperature associated with plumes . Even with these prior studies, how the various tectonic events influenced the major tectonic units and boundary zones of the NCC is still not fully understood. Understanding the complex tectonic evolution of this craton demands detailed knowledge about the deep structure.
The crust-mantle boundary, known as the Moho, is one of the main first-order discontinuities in the Earth's interior and may exhibit complex evidence of past tectonic processes. Accurate Moho depths provide key information for studies such as tectonic and geodynamic modeling and further high-resolution tomography. In the NCC, detailed knowledge of the Moho depths can also provide important constraints on the spatial extent of lithospheric thinning and destruction, that is, whether the crust has also been thinned during the cratonic reactivation and how the thinned crust distributes in association with regional tectonics and lithospheric deformation.
Previous studies on the deep structure of the NCC primarily focused on the distinctly different characteristics in the east and west or the E-W structural variations in the NCC, many of which were limited along linear arrays/profiles with an assumption of two-dimensional structure beneath the region (e.g., Chen et al., 2006a, 10.1029/2018JB016122 2006bChen, 2010;Deng et al., 2007;Jia et al., 2001Jia et al., , 2014Jia & Zhang, 2005;Li et al., 2011;Wang et al., 2013;Wei et al., 2011;Wu et al.,2011;Zheng et al., 2006Zheng et al., , 2009 or to relatively large-scale structures (>100 km; e.g., Duan et al., 2016;He et al., 2014;Huang et al., 2009;Li et al., 2006;Li et al., 2014;Wei et al., 2016). N-S structural changes or small-scale features within or between major blocks of the NCC have not been thoroughly investigated. There are limited constraints on N-S Moho variations and small-scale features (e.g., Deng et al., 2007;He et al., 2014;Jia et al., 2001Jia et al., , 2014Jia & Zhang, 2005;Li et al., 2006;Li et al., 2014;Wang et al., 2013;Wei et al., 2011Wei et al., , 2016Wu et al., 2011;Yang et al., 2018;Zheng et al., 2006Zheng et al., , 2009 in the area of interest (Figure 1). The receiver function (RF) approach is the major tool in mapping the Moho in passive-source studies. However, RFs are notoriously challenging to interpret in regions of large sediment basins, for example, the Bohai Bay Basin, due to the potential contamination of the Moho signal by strong reverberations from sedimentary layers (e.g., Chen et al., 2006a;Zheng et al., 2006). On the other hand, active-source studies can potentially image the Moho with higher resolution because of the higher frequencies of data. Duan et al. (2016) integrated 42 previous deep seismic sounding profiles and provided a unified Moho depth map in the central-eastern NCC. The map is, however, relatively smooth over the entire region as a result of the interpolation of nonuniformly distributed 2-D profiles, mostly showing the apparent E-W variations in the Moho depth. Furthermore, detailed comparisons with crosschecks and evaluations between passive-and active-source studies on the Moho depths are still not available, due largely to their incoordinate data coverage and inherent method limitations.
In this study, the densely distributed stations in the northeastern NCC and the resulting dense data coverage ( Figure 2) allow us to construct a high-resolution Moho depth map for the region. We adopted a 2-D wave equation-based migration method (Chen et al., 2005a(Chen et al., , 2005b with an integrated 3-D crustal velocity model for both Ps conversions and surface-reflected multiples to better constrain the absolute depths of the Moho. We then systematically compared our Moho depth map with previous results from both passive-and active-source studies to verify the observations. Finally, we investigated the correlations of the Moho depth with surface geology, gravity anomalies, and lithospheric thickness in detail, and we discuss their tectonic implications on the Phanerozoic destruction of the NCC.

Data and Methods
In this study, we used approximately 2 years of teleseismic records (from December 2006 to September 2008) collected by 198 broadband stations from a temporary dense seismic array (North China Seismic Array; Chi-nArray, 2006) in the northeastern NCC ( Figure 2). The North China Seismic Array consists of two dense linear arrays with an interstation spacing of ∼10 km in the north and south and an average spacing otherwise of ∼35 km between nearly uniformly distributed stations ( Figure 2). We selected 341 high-quality teleseismic events with magnitude above 5.5 and epicenter distance between 30 • and 95 • . Almost 70% of the events are concentrated within the back-azimuth range of 100 • -220 • , while others fall within the back-azimuth ranges of 0 • -100 • and 220 • -360 • separately (inset Figure 2). We then computed P RFs by iterative deconvolution in the time domain (Ligorria & Ammon, 1999) with a Gaussian width of 2.5 (the corner frequency is ∼0.8 Hz). After careful visual inspection, 28,220 RFs with high signal-to-noise ratio were manually selected for imaging. Examples of the RFs are shown in Figure S1.
According to the overall strike of the structures in the study region and the coverage of piercing points of P to S conversions at depths close to the Moho, we selected 10 profiles perpendicular (P1-P10) and 9 profiles parallel (P11-P19) to the NE-SW trending direction of surface structures for imaging ( Figure 2). The profiles P1 and P2 lie along the two linear dense arrays in the north and south, respectively ( Figure 2). We then employed a 2-D wave equation-based RF poststack depth migration method for Ps (Chen et al., 2005a) and PpPs phases (Chen et al., 2006a), respectively, to construct the crustal structural images along the 19 profiles.
Following Chen et al. (2005a), the Ps migration of RFs was performed with the time-domain common conversion point (CCP) stacking and subsequent poststack depth migration. The first-step stacking builds on the idea of common middle point stacking of reflected data in exploration seismology, which was pioneered and modified to stack RFs for the depth of one target discontinuity (e.g., Moho) at a time by Dueker and Sheehan (1997). Here we processed the RF data by moveout correction, binning, and stacking for different depths simultaneously in the time domain, analogous to the space-domain CCP stacking introduced by Zhu (2000). Therefore, we termed the processing as time-domain CCP stacking (Chen et al., 2005a). The major procedures are as follows: First, each RF was moveout corrected for all the depths of target discontinuities to the case of horizontal slowness p = 0, based on the Ps delay times of individual discontinuities calculated from a 1-D reference model; Next, all RFs were segmented and sorted into different bins and subsequently stacked with respect to each depth according to both the horizontal locations and depths of their sampling points. The resultant CCP-stacked traces can be regareded as a good approximation to a zero-offset wavefield data set on the surface, that is, Ps converted phases propagating vertically from the sampling points of the discontinuities to the surface directly above them (Chen et al., 2005a). This stacked RF data set was then backward propagated to the image space through frequency-domain wave equation-based wavefield extrapolation (Chen et al., 2005a).
Above procedures were applied to the migration of surface multiples in RFs with small modifications, in a similar way as the previous study by Chen et al. (2006a). For the case of PpPs migration, the delay time-distance relation for PpPs was used in moveout correction and the sampling locations of PpPs signals were considered in binning and time-domain CCP stacking. During wavefield extrapolation, the migration velocity for PpPs (Chen et al., 2006a) was adopted to backward propagate the stacked section to the image space.
We carefully selected the major parameters, that is, frequencies and bin sizes, to achieve the best performance of migration with this data set. The frequency content of the RFs was set to 0.03-1 Hz for Ps migration and 0.03-0.35 Hz for PpPs migration. We considered rectangular stacking bins with a fixed length of 30 km perpendicular to each profile. Therefore, data overlaps exist between profiles, which contribute to the mutual verifications of neighboring profiles and thus help to construct a reliable 3-D Moho depth map from 2-D images. The lateral half bin width for the rectangular bins varies at different depths (with wider bin sizes at greater depths) and also changes at the same depth based on the data density. The minimum number of RFs in each bin was set to be 200 based on the data density at depths of ∼30-40 km for effective noise suppression. For the sake of balancing uneven data coverage along 19 profiles, the final lateral half bin width along the profiles is approximately 34 km on average at 32 km depth. Such choice leads to construct continuous and coherent structural images even for areas where data are very sparse (e.g., NE-SW profiles). However, the observations with quite dense data coverage along profiles P1 and P2 were more laterally smoothed at the same time, where the spatial resolution of the images for the target was decreased. At the edges of profiles, gradually widened bins were adopted to compensate the reduced number of data. We are aware of the problem that the very large bin size at the edges is likely to cause aliasing and strong smoothing effects, so the results there were excluded from our interpretation. Figure 3 shows an example of CCP-stacked sections for Ps conversions and PpPs multiples along profile P1, with associated data density and bin sizes at 32 km depth. The positive Ps signals (∼5 s, Figure 3b) and PpPs multiples (∼15-20 s, Figure 3c) from the Moho are clearly visible in the stacks and can be more distinctly identified in migration images (at ∼30-40 km depth, Figure 4). As expected, the PpPs signals is also present in the stacked section for Ps (Figure 3b), but is more pronounced and continuous in the one for PpPs (Figure 3c). These results confirm that the data quality is sufficient and that the selected parameters are reasonable for stacking. More details about the method and parameter setting can be found in Chen et al. (2005aChen et al. ( , 2005b. We used both a 1-D average velocity model and an integrated 3-D crustal velocity model for RF imaging to better constrain the absolute depths of the Moho. Initially, we adopted the 1-D IASP91 velocity model (Kennett & Engdahl, 1991) superposed by a 1-D crustal model of eastern China ( Figure S2; Ai & Zheng, 2003) as an average velocity model for RF migration. We then performed moveout correction of RFs by considering the individual 1-D crustal model for each station. The station-based crustal models were constructed by first incorporating the 3-D crustal Vs model (0-50 km) constrained by surface wave dispersion (periods from 4 to 30 s) from ambient noise in the study region (Fang et al., 2010). Bulk crustal Vp/Vs ratios were set (1) based on the results of H-stacking of RFs for stations outside the Mesozoic-Cenozoic basins, (2) to be 1.76 for stations in the sedimentary basins in the TNCO and western NCC, which is an average value of all bedrock stations in these areas, and (3) to be 1.78 for stations located in the Bohai Bay Basin based on the previous Poisson's ratio study in this region (Ji et al., 2009;Wang et al., 2009). Second, we superposed the sedimentary velocity structures obtained by waveform modeling analyses of P RFs ) on the crustal model for each station. In particular, some sedimentary structures were slightly corrected based on our P RF migration images and previous local geological surveys, oil or mineral exploration, and geothermal studies (Lin & Gong, 2005;Ma, 1989), such as in the Datong Basin and the Yanqing-Huailai Basin (see locations in Figure 1b). Finally, we used empirical relationships given by Brocher (2005) to compute densities from velocities for the crustal model. An example of model corrections is presented in Figure 4 and discussed in the next section.

Factors Affecting Migration Results
In this study, we considered not just Ps conversions but also surface-reflected PpPs multiples in the migration. The large Bohai Bay Basin in the eastern NCC and several small basins (Xuanhua Basin, Yanqing-Huailai Basin, Yangyuan-Yuxian Basin, and Datong Basin) located in the Shaanxi-Shanxi rifts in the TNCO are covered by thick sediments (Figure 1b), resulting in strong sediment reverberations in the RFs (e.g., stations L232, A003, and A602, Figure S1). Complex intracrustal structures, such as the well-known low-velocity zone beneath the Tangshan area (Fang et al., 2010), also severely disrupt the Moho signals in the RFs (e.g., stations K047-K051). In some cases, intermittent or even invisible Ps signals of the Moho make it very challenging to trace the Moho continuously, particularly for profiles that transverse the Bohai Bay Basin (>550 km distance along profile P1 in Figures 4b and 4d, >400 km distance along profile P2 in Figure 6. Map of the Moho depth derived from receiver function imaging. Seismicity is represented by red circles. The red hollow stars mark the locations of mantle xenolith localities. The yellow diamonds denote four cities including Zhangjiakou, Jining, Datong, and Yuxian, which are also marked in Figure 2. The large yellow dashed circle represents the ZJDY area. The white solid lines denote the four representative profiles for discussion. The gray shadow zone is a noninterpreted region with imaging results of lower confidence (see the last paragraph in section 3.1). Other symbols and labels for tectonic units are the same as in Figure 1. Figure 5b, >200 km distance along profile P17 in Figure 5h, and the corresponding portions along profiles P5-P10, P16, P18, and P19 in Figure S3). However, reverberations from sedimentary structures seem not to significantly interfere with the PpPs signals, as shown by the RF waveforms of stations A003 and L232 in the Bohai Bay Basin ( Figure S1) and from the stacked section for profile P1 ( Figure 3) and migration images for profiles P1 (Figure 4), P2, and P17 ( Figure 5). Compared with the Ps images (Figures 4b and 4d), the PpPs images (Figures 4c and 4e) show a weaker sensitivity to basin structures. Previous RF studies in the NCC revealed the similar superiority of the PpPs phase over the Ps phase in subsurface structure imaging at basin areas covered by thick sediments (Chen et al., 2006a;Zheng et al., 2006). Therefore, we conducted poststack depth migration on both Ps and PpPs phases for all the profiles and identified the Moho by analyzing both images. Figure 4 shows the imaging sections for profile P1 as an example of model correction. Comparing the migration images using the 1-D average velocity model (Figures 4b and 4c) with the images corrected by 3-D crustal velocity model (Figures 4d and 4e), we observe apparent improvement of the imaging results, particularly in the areas covered by sediment (e.g., at distances of ∼100, ∼250, ∼350, and >400 km along the profile, Figure 4).
Generally speaking, the Moho can be traced coherently and robustly on both Ps and PpPs images but may appear unclear beneath thick sediments (e.g., at 300-400 km distance along profile P17 in Figures 5g-5i) or complex intracrustal structures (e.g., at ∼600 km distance along profile P1 in Figure 4). We conducted PpSs + PsPs migration for profiles traversing the Bohai Bay Basin ( Figure S4) to further verify the final Moho depths. The detailed processing is very similar to the case for PpPs migration (refer to the last section) and follows the previous publication by Wang et al. (2018). Additionally, uneven ray coverage due to the events used mainly from SE (inset Figure 2), data absence locally (e.g., at Datong Basin along profiles P9 and P10 in Figure 2), or sharp changes of the Moho depths (e.g., at 300-500 km distance along profile P2 in Figures 5b and 5c) make it challenging to trace the Moho continuously. Given that converted and topside reflected phases from a dipping discontinuity exhibit distinct behaviors for different incident P wave directions (Chen et al ., 2006b), we migrated these two kinds of phases with different back-azimuth ranges of data sets and compared the corresponding images in detail to better constrain the geometry of the Moho (e.g., see Text S4 and images of profile P6 in Figure S5). Moreover, we calculated theoretical delay times for all Moho phases using the established station-based models and adjusted the Moho depths to fit the observed delay times ( Figure S1). Finally, we picked the consistent Moho depths from the images by combining the above analyses for each profile (e.g., Figures 4, 5, and S3) and further integrated the results of individual profiles into a 3-D map of the Moho depth for the northeastern NCC ( Figure 6).
In the following discussion, we focus only on the results that were cross-checked by two profiles or with relatively dense stations. We do not consider the following situations: (1) low-quality RF data (e.g., from station K051 in the Tangshan area); (2) both Ps and PpPs images at Moho depth strongly contaminated by known crustal multiples with no solid confirmations from other tests (e.g., at >600 km distance along profile P1 in Figure 4); (3) unconstrained crustal velocity structure for a distance of >50 km along a profile (e.g., at 0-100 km distance along profile P6 in Figure S3); (4) obvious aliasing and smoothing effects at the edges of a profile ( Figure S3).

Moho Depths Along Profiles
The Moho can be coherently identified at depths of 25 to 45 km along 19 profiles through the integration of migration images using Ps conversions and surface-reflected multiples (Figures 4, 5, and S3). Both the Ps and PpPs images are dominated by strong positive Moho signals at consistent depths (∼32-44 km) beneath the mountain ranges, which can be clearly observed and continuously traced with little multiple interference or noise-induced artifacts (e.g., profiles P3, P4, P11-P15 in Figures 5 and S3). The Ps images at Moho depth, however, exhibit large ambiguity with intermittent or even invisible signals beneath the basin areas, in particular the Bohai Bay Basin. However, these disruptions can be well compensated by using much more continuous Moho signals in the PpPs images (e.g., profiles P1, P2, P5-P10, P16-P19 in Figures 4, 5, and S3). Note that the Moho topography in the PpPs images appears to be smoother than in the Ps images due to laterally wider sampling and lower frequencies of PpPs multiples than Ps conversions.
The Moho along the NW-SE trending profiles (P1-P10) shows an eastward shallowing, with the most rapid depth variations (up to 10 km over a distance of ∼100 km) coincident with the boundaries of the TNCO (Figures 4, 5a-5c, and S3). We consider principally the two profiles P1 and P2 with the densest coverage of RF data. While a gentle reduction in Moho depth is observed along profile P1 from the Yin and Yan Mountains (∼40 km) to the northern margin of the Bohai Bay Basin (∼30 km; Figure 4), the Moho depth along profile P2 in the south varies much more rapidly from ∼44 km under the Ordos Block to ∼26 km within the Bohai Bay Basin (Figures 5a-5c).
The observed Moho in the NE-SW oriented profiles (P11-P19), however, are much more complex, with more localized undulations in the Moho topography, including a deepening from NE (∼38 km) to SW (∼44 km) on the western side of the NSGL (P11-P13) in addition to a southward shallowing on the eastern side of the NSGL (P14-P19; Figures 5d-5i and S3). We focus on the two NE-SW trending profiles P12 and P17 with very high-quality images and assembling major structural features discussed in the interpretation part. For the NE-SW trending profile P12 in the TNCO, the Moho is nearly flat (∼40 km depth), with only a slight shallowing of ∼2 km in the middle beneath the Shaanxi-Shanxi rift systems (at distance ∼250 km in Figures 5d-5f). In contrast, profile P17 with the same orientation in the eastern NCC is characterized by a much thinner crust and a gradual decrease of Moho depths from NE (∼32 km beneath the Yan Mountains) to SW (∼26 km beneath the Jizhong Depression; Figures 5g-5i). Figure 6 shows lateral variations in the Moho depth in the study region. In general, the Moho depths are distinctly different on the opposite sides of the NSGL, as deep as >40 km in the northwestern mountain ranges but as shallow as <30 km in the southeastern basin areas, with a rapid Moho uplift of ∼10 km beneath the NNE-trending Taihang Mountain ranges in the middle. Sudden changes of crustal thickness are also observable near the boundaries of the TNCO on both sides, suggesting that the Moho depth variations are consistent with the threefold subdivision of the craton (Zhao et al., 2005). The Moho depths beneath the Yin and Yan Mountains in the north gradually decrease from ∼40 km in the west to ∼30 km in the east, with variation much smoother than that (from ∼44 to ∼26 km) to the south.

The 3-D Moho Depth Variations
In addition to the dominant E-W (NW-SE) changes, small-scale N-S (NE-SW) variations in the Moho depth also appear in Figure 6. In the central to western NCC (mainly to the west of the NSGL), the Moho depths to the south of Yin-Yan Mountains display a prominent increase from ∼38 km north of Fanshi (red star in Figure 6) to >42 km farther south. In the eastern NCC, the crust beneath the Yan Mountains in the NE (∼32-40 km) is thicker than that beneath the Bohai Bay Basin in the SW (∼26-32 km). Moreover, the Moho 10.1029/2018JB016122 within the Bohai Bay Basin further shallows from the NE (∼32 km) to SW (∼26 km), showing apparent variations from the margin to the interior of the basin.
Specifically, the imaged Moho undulations in Figure 6 show a strong correlation with the geological and tectonic divisions in the NCC. The thickest crust (>44 km) appears in the northeastern margin of the Ordos Block and slightly thins to ∼42 km toward to its eastern boundary. A local Moho uplift (∼38 km) appears around the Zhangjiakou-Jining-Datong-Yuxian area (yellow dashed circle in Figure 6) in the central to western NCC, mainly confined to the west of the NSGL. For simplicity, we refer to this region as the ZJDY area hereafter. This local Moho undulation also correlates with distributions of the Cenozoic basalts (Jining, Hannuoba, Yangyuan, and Datong, red stars in Figure 6). In the eastern NCC, the crustal thickness of the Bohai Bay Basin is ∼26-32 km. Within the Bohai Bay Basin, Moho undulations of ∼3 km also exist between the Jizhong Depression and Cangxian Uplift. The Moho under the southern Jizhong Depression is the shallowest, reaching ∼26 km, where is covered by the thickest sediment layer in the Bohai Bay Basin (∼10 km; Fang et al., 2010;Jia et al., 2014;Zheng et al., 2006).

Comparisons with Previous Studies
The pattern of variation in the Moho depth observed in Figure 6 generally agrees with previous passive-source observations of the region (e.g., Chen et al., 2006a;Wang et al., 2013;Wei et al., 2011;Wu et al., 2011;Zheng et al., 2006Zheng et al., , 2009, but these changes are now shown with unprecedented detail. First, the migration images in this study show improvements in constraining the Moho depths relative to the conventional CCP-stacked Ps images using the same data sets (e.g., compare our images for profile P1 in Figures 4d and 4e and for profile P2 in Figures 5b and 5c with the corresponding CCP images shown in Figures S6a and S6b, respectively; Wang et al., 2013;Wu et al., 2011). These improvements arise because we employed the advanced migration technique and considered a 3-D crustal model while integrating analyses on images of multiple phases. Our results also provide a robust 3-D Moho image that is comparable to previous results but with denser sampling (beneath the region with an average station spacing of ∼35 km in this study but ∼10 km in the previous study, Figure S7; Zheng et al., 2006Zheng et al., , 2009. Refer to Text S5 for more details. The imaged Moho variations across the study region ( Figure 6) also generally agree with the Moho map from the active-source study in Figure 7a (also see results along individual profiles in Figures 4, 5, and S3). Duan et al. (2016) integrated 42 deep seismic sounding profiles to build a 3-D crustal velocity model for the central-eastern NCC, namely, HBCrust1.0. Their study region spatially overlaps that of this study, which brings the possibility to compare Moho depths from the two independent data sets of different nature. Figure 7a shows the derived Moho depth distributions from HBCrust1.0, which are relatively smooth, probably due to the interpolation method (Duan et al., 2016;Waldhauser et al., 1998). Our results coincide better with the reported Moho depths by individual active-source profiles in most geological units (e.g., Deng et al., 2007;Jia et al., 2001Jia et al., , 2014Jia & Zhang, 2005), particularly beneath the southern Jizhong Depression where the Moho depth is found to be ∼26 km in this study ( Figure 6) but ∼28 km in the previous one (red circle in Figure 7a). Different sensitivities to the structure of RF and active-source data should be taken into account regarding the differences between two results (Figures 6 and 7a). The high-frequency reflected phases (e.g., PmP) provide excellent resolution of crustal Vp reflectors, whereas RF (e.g., Ps) is primarily sensitive to changes in the Vs structure but depends weakly on the crustal Vp structure, and a gradual velocity transition is also resolvable with RF data. The ZJDY area that shows prominent disagreements of Moho depths in Figures 6 and 7a probably underlies a ∼10 km thick crust-to-mantle transition zone (at depths of ∼32-42 km) inferred both seismically  and petrologically (Chen et al., 2001). In this study, the local Moho undulations (yellow dashed circle in Figure 6) and diffused and weak PpPs signals (e.g., stations A503 and A602 in Figure S1) in the region might also indicate this structural anomaly (discussed in the next section). Furthermore, we speculate that the depth of ∼38 km in Figure 6 corresponds to the average depth of the gradual crust-to-mantle transition zone, while ∼42 km in Figure 7a is likely its bottom boundary beneath the ZJDY area. More details can be found in Text S6.

Correlations of the Moho Depths With the Surface Geology, Gravity Field, and Lithospheric
Thickness and Tectonic Implications The lateral variations in the Moho depth observed in this study correlate with the geological and tectonic units in the NCC (Figure 6), suggesting that the Moho undulations may be directly related to the complex regional tectonics, particularly the Phanerozoic lithospheric destruction of the NCC. To examine this idea,  (Duan et al., 2016). The red lines are two separate DSS profiles published by Jia et al. (2014; The Wendeng-Alashan Profile) and Deng et al. (2007; The Zhucheng-Dingxian-Tuoketuo Profile). Both profiles cross the Jizhong depression and the reported Moho depth there is approximately 28 km. (b) LAB depth map from S receiver function study in the study region (Chen, 2010). The blue lines are the four representative profiles for discussion.
we selected four individual imaging profiles that traverse different parts in the study region and further analyzed the Moho depth variations with the topography, geology and tectonics, gravity field (Ma, 1989), crustal shear wave velocity structures (Fang et al., 2010;Wu et al., 2014), and LAB depth (Chen, 2010) from the surface to the bottom of lithosphere ( Figure 8; the profile locations are marked in Figures 6 and 7b).

E-W Difference of the Moho Depths Across the Northeastern NCC
A dominant feature revealed in this study is the E-W (NW-SE) variations in the Moho depth across the region, which are approximately concordant with the changes in the surface geology, crustal structure, and LAB depths (Figures 6, 7b, 8a, and 8b). Such observations have been reported by previous studies (e.g., Chen, 2010, and references therein), indicating that the first-order feature of E-W variation in the crustal structure is coupled with the lithospheric processes, that is, fundamental destruction involving both the crust and lithospheric mantle and the complex but more localized lithospheric modifications to the east and west of the NSGL. Additionally, the eastward Moho shallowing is generally coincident with the E-W changes of gravity field from the strong negative Bouguer anomalies (< −100 mGal) in the central and western NCC (at distance of <400 km) to the weakly negative to positive Bouguer anomalies (> −20 mGal) in the eastern NCC (at >400 km distance; Figures 8a and 8b). However, there are secondary lateral variations of the two, which are mutually inconsistent in the Bohai Bay Basin (Figures 8a and 8b). This finding suggests either  (Ma, 1989), surface altitude, tectonic divisions and major cities, Moho depths (this study), crustal shear wave velocities (Fang et al., 2010;Wu et al., 2014), and LAB depths (Chen, 2010)  lateral variations in the crustal density or contributions from the mantle to the gravity observations beneath the Bohai Bay Basin.

N-S Difference of the Moho Depths to the West of the NSGL
On the western side of the NSGL, the relatively shallow Moho (∼38 km, Figures 6 and 8c) in the ZJDY area may reflect a local structural anomaly in the northern margin of the central to western NCC. This area geographically correlates with a ∼10 km thick crust-to-mantle transition zone inferred by P RF waveform inversions in the junction of the Taihang Mountains and the Yin-Yan Mountains . This inference was also demonstrated by observations on basalt-borne xenoliths in the Hannuoba area (see Figures 1 and 6 for its location), in which spinel lherzolites and mafic granulites are intermixed at consistent depths of ∼32-42 km (Chen et al., 2001) with seismological constraints . Along profile P12 in Figure 8c, the Vs model (the fourth panel) shows a very gradual velocity change close to the depths where the Moho is expected in the region (at a distance of ∼200-400 km), contrasting with the sharp velocity jump in either the northern (beneath the Yan Mountains at <200 km distance) or southern segments (beneath the interior of the TNCO at >400 km distance). This region is also characterized by a local increase of ∼40 mGal in gravity relative to the adjacent regions (first panel in Figure 8c). All these observations, considered together with the inferred distinctly anisotropic lowermost crust (Cheng et al., 2013;Fu et al., 2015), the presence of consecutive low velocities from the lower crust to upper mantle (Fu et al., 2015;Jiang et al., 2013), low densities , and high-conductivity anomalies , support the inference that complex crust-mantle interactions may have taken place and perhaps are related to magmatic underplating and hot mantle upwelling processes. Given the spatial coincidence of this region with the distributions of Cenozoic basalts (Figure 6; Tang et al., , 2013 and the locally relatively high heat flow observations (Hu et al., 2000), we suggest that lithospheric reworking and related dynamic processes of this region are likely ongoing and even currently active.
In the interior of the TNCO, however, the Moho depths are deeper (>42 km) than that to its north (∼38 km), which may indicate different subsurface structural features from the ZJDY area (approximately bounded by Fanshi in Figure 6 and at a distance of 400 km along profile P12 in Figure 8c). According to Zheng et al. (2008), the crust-to-mantle transition zone is estimated as ∼5 km thick in the interior of the TNCO, which is much thinner than that to the north (∼10 km). The crust in the same region was revealed to be of high velocity by surface wave tomography (at ∼400 km distance, fourth panel in Figure 8c; Fang et al., 2010;Fu et al., 2015), showing the different structural features relative to that beneath the northern margin of the TNCO. Cenozoic rifting near this region is also weaker than in the northern extensional regions (Zhang et al., 2003). No direct evidence on crust-mantle interactions has been reported by petrological or geochemical studies in the interior regions such as Fanshi and Hebi. Recent dating and isotopic analysis of granulite xenoliths attests to the existence of Archean lower crust beneath the TNCO Zhang et al., 2012). Based on these observations, we suggest that the intense crust-mantle interactions are absent in the interior of the TNCO or that the crust may have retained its Archean nature and still be relatively stable, with only slight thinning resulted from tectonic extension in the region. Note that an analysis of the finer-scale structure of the crust-to-mantle transition zone is outside the scope of this paper, but the general features revealed here do not conflict with the results of previous studies.
To the west of the NSGL, the LAB depths increase from ∼110-120 km north of Fanshi to >130 km farther south (Figure 7b), showing similar N-S variations as in the Moho depth ( Figure 6). The uppermost mantle beneath the central part of the TNCO was detected to be of high velocity  and high electrical resistivity (Wei et al., 2008) as opposed to the low velocity  and high conductivity  characterizing the northern margin of the TNCO. Observations from the peridotite xenoliths from Cenozoic basalts also reveal isotopic and chemical heterogeneities between these two regions. The mantle peridotites in Jining, Hannuoba, Yangyuan, and Datong (see Figures 1b and 6) are dominated by fertile low-Mg# (Fo < 92) compositions, which are characteristic of younger ocean-like lithospheric mantle (Tang et al., , 2013 and similar to that in the eastern NCC. However, the existence of highly refractory harzburgites with high-Mg# (Fo > 92) in peridotites from Fanshi and Hebi was documented by petrology and geochemistry studies (Tang et al., , 2013Xu et al., 2008;Zhang et al., 2012). These results suggest the lithospheric mantle was less affected in the interior of the TNCO, contrasting with the considerably modified lithospheric mantle in the northern margin. All of the above geophysical and geochemical observations in the interior of the TNCO suggest structural features similar to those of the stable cratonic lithosphere in the western NCC, that is, the Ordos Block, but distinctly different from those of the significant destroyed lithosphere in the eastern NCC. It is rational to conclude that the thinned lithosphere in the central part of the TNCO is probably cold, dry, and rigid, representing a cratonic remnant that has not been totally destroyed during the craton reactivation.
In summary, the different characteristics in the lithospheric mantle in N-S inferred by geophysics and geochemistry are consistent with the crustal structural heterogeneity and also correspond to the concordant N-S variations in the Moho ( Figure 6) and LAB depths (Figure 7b) to the west of NSGL. This result suggests that the Phanerozoic craton rejuvenation may have affected the central to western parts of the NCC but spatially unevenly on both the crust and lithospheric mantle. The whole lithosphere possibly underwent high-degree lithospheric modifications and the related dynamic processes may still be active beneath the northern margin of the NCC, whereas the lithosphere was less affected and the cratonic remnant has probably been partially preserved in the current crust and lithospheric mantle beneath the interior of the TNCO.

N-S Difference of the Moho Depths to the East of the NSGL
In the eastern NCC, the relatively thick crust (∼32-40 km) beneath the Yan Mountains in the north contrasts with the dramatically thinned crust (∼26-32 km) beneath the Bohai Bay Basin in the south (Figures 6  and 8d). The strong and concentrated PpPs signals of the Moho from stations located at the Yan Mountains (e.g., station A303 in Figure S1) also indicate a fairly sharp Moho beneath this region, as opposed to the gradual crust-mantle transitions that characterize the Bohai Bay Basin to its south . The E-W oriented fold and thrust structures in the Yan Mountains (Davis et al., 2001) differ from the predominant NE-SW trending widespread rifting and extensional structures in the Bohai Bay Basin (Ren et al., 2002), which is covered by very thick sediments (e.g., at >200 km distance along profile P17, fourth panel in Figure 8d). Furthermore, comparing the overall slow structures (Tang & Chen, 2008) and substantially thinning lithosphere (<100 km, Figure 7b) beneath the Bohai Bay Basin in the interior of the craton, distinct high-velocity structures in the crust and uppermost mantle have been identified with surface wave tomography (Tang & Chen, 2008) beneath the Yan Mountains, along with the relatively deep LAB in this region (>100-120 km, Figure 7b and that further north in Chen et al., 2008). These distinct structural features are also consistent with reported low heat flow values (∼30-45 mW/m 2 ; Hu et al., 2000) and negative gravity anomalies in the region (at <200 km distance, Figure 8d), different from those in the south (heat flow is ∼53 mW/m 2 ; Hu et al., 2000); the gravity anomaly is close to zero at >200 km distance in Figure 8d). All of the above observations suggest that the Yan Mountains, as the northern boundary of the eastern NCC, have possibly experienced distinctly different deformational history within the lithosphere during the Mesozoic-Cenozoic from that in the Bohai Bay Basin. In particular, the multiple phases of N-S contractions in Mesozoic probably have brought profound structural complexities in the Yan Mountains, resulting in a different response with the Bohai Bay Basin in the subsequent extensional environment since the Late Mesozoic. Consequently, structural heterogeneities are formed from the surface to the deep lithosphere beneath the Yan Mountains relative to the adjacent Bohai Bay Basin.
The Moho beneath the Bohai Bay Basin further shallows from NE (∼32 km) toward to SW (∼26 km), while the lithosphere presents a southward thickening, and the Bouguer gravity anomalies display a moderate decrease along profile P17 (Figure 8d). Within the Bohai Bay Basin, neither the change in crustal structure nor the Moho depth undulations can be responsible for observed long-wavelength gravity decrease along the profile (at >200 km distance, Figure 8d). Thus, we attribute the decrease to the deeper lithospheric structure beneath the Moho, which requires the underlying lithospheric mantle to be less dense and more buoyant in the interior of the Bohai Bay Basin.
Such an increase to the south in the buoyancy of the lithospheric mantle may indicate either compositional heterogeneity or thermal-induced density variation or both in the lithospheric mantle beneath the Bohai Bay Basin. Considering the dominant fertile lithospheric mantle that has been well demonstrated in the eastern NCC, refractory compositions constrained by geochemical studies (e.g., Wu et al., 2008;Xu et al., 2008; may indicate the existence of local chemical depletion in the thinned lithospheric mantle beneath the eastern NCC. The local chemical depletion in composition may be either recycled cratonic lithospheric materials that accreted to the lithosphere later (e.g., Wang et al., 2016) or in situ remnants. This effect probably contributes to the buoyancy of the lithospheric mantle in the basin interior to some degree, which is one possibility to explain the correlations in seismic and gravity data in Figure 8d. Alternatively, the observed N-S variations may be associated with thermally induced density variations in the lithosphere beneath the Bohai Bay Basin. Consensus is that the region of interest within the Bohai Bay Basin has been cooling during its thermal subsidence process after an intensive rifting in the Paleocene-early-Eocene (Allen et al., 1997;He & Wang, 2003). However, the cooling rate of the lithosphere beneath the basin might be spatially variable, particularly between the margin areas adjacent to the Yan Mountains in the north and the interior of the Bohai Bay Basin in the south. The present lithosphere beneath the Yan Mountains, the northern boundary of the eastern NCC characterized by low heat flow (Hu et al., 2000) and high crustal to upper mantle velocities (Tang & Chen, 2008), might have been colder than other areas in the eastern NCC and caused more efficient cooling in the northern margin of the Bohai Bay Basin. However, the cooling beneath the cratonic interior is less efficient due to its greater distance to the surrounding mountains. Such uneven cooling in the basin may lead to thermally induced density variations, which accounts for the increased buoyancy southward in the lithospheric mantle beneath the Bohai Bay Basin. Regarding the above analyses, we propose that the Cenozoic lithospheric mantle beneath the Bohai Bay Basin may be laterally heteroge-neous chemically and/or thermally even though having been extensively thinned. A better understanding of the relative contribution from these two factors requires further related evidence.

Dynamic Causes of the Lithospheric Destruction in the Northeastern NCC
The primary structural changes in the E-W (NW-SE) direction appear concordant from the surface to the bottom of the lithosphere (Figures 6, 7b, 8a, and 8b). We attribute the general E-W structural variations to the Paleo-Pacific plate subduction under eastern Asia since the Mesozoic that caused significant thinning and fundamental destruction on both the crust and lithospheric mantle in the eastern NCC, in agreement with previous studies (e.g., Chen, 2010;Zhu et al., 2011, and references therein). From the north (NE) to the south (SW), however, the Moho depth distribution also shows substantial variations, and its correlations with other observations appear to be more complex (Figures 6, 7b, 8c, and 8d). These N-S directed structural changes cannot be fully explained by the Mesozoic Pacific subduction.
The N-S directed structural variations in the eastern NCC may reflect the imprints of earlier multiple N-S compressional deformations along the northern margin of the NCC, associated with the amalgamation of the NCC with the Siberian craton that induced the formation of the Yan Mountains (Davis et al., 2001;Meng, 2003). The distinctly different structural features between the boundary areas (i.e., Yan Mountains) and the craton interior (i.e., Bohai Bay Basin) in the eastern NCC also indicate their variable responses to the subsequent lithospheric extension in Late Mesozoic to Cenozoic (Ren et al., 2002).
The N-S oriented structural changes to the west of the NSGL may be associated with the combined effects of multistage tectonic events over the long term from the formation of the craton to the present day. The TNCO and the EW-trending Khondalite Belt close to the northern boundary of the western NCC probably have been mechanically weakened relative to adjacent cratonic nucleus beneath the Ordos since their formation in the Paleoproterozoic (Zhao et al., 2005). These preexisting structural belts may have been repeatedly reactivated by successive tectonic events during the long-term evolution of the NCC. Thus, the underlying lithosphere may have been weakened and deformed to some extent during that time. Recently, the India-Eurasia collision that drove the differential counterclockwise rotations of the Ordos with respect to its surrounding regions (Zhang et al., 1998(Zhang et al., , 2003 has probably further contributed to deformation of the weakened lithosphere under the preexisting structural belts and thus intensified the structural contrasts between the margin areas and interior of the NCC. In particular, the inferred active underplating and upwelling processes around the ZJDY area, along with the extensional rifting and strong seismicity (Figure 1b), may represent the last stage of Cenozoic lithospheric modifications involving both crust and lithospheric mantle in that region. The ZJDY area is close to the ancient collisional Khondalite belts and the northern boundary of the NCC, where the lithospheric structure may have been inherently weaker. Spatially, this region is mainly confined to the west of the NSGL in the northern margin of the NCC, a region that has been less influenced by the Mesozoic Paleo-Pacific Plate subduction processes. Furthermore, the India-Eurasia collision has been proposed to exert a major far-field influence on the lithospheric tectonics of the central and western NCC during the Cenozoic (Chen, 2010;Menzies et al., 1993;Xu, 2007), whereas the influence of Pacific subduction on the NCC and the eastern China continent may have relatively weakened through the Cenozoic (Chen, 2010). Given the spatiotemporal concordance between the Cenozoic lithospheric modification with the tectonic processes induced by the India-Eurasia collision, we propose that the Cenozoic lithospheric modifications at the boundaries of the Ordos Block may have initiated in the northern marginal areas and may be predominantly influenced by the India-Eurasia collision, less intensive and much later than the Mesozoic lithospheric reactivation and destruction in the eastern NCC. This inference also corroborates the diachronous lithospheric thinning scenario suggested by Xu (2007), based primarily on petrological and geochemical data from mantle xenoliths in the northern part of the Shaanxi-Shanxi rift system. However, the lithosphere beneath the central part of the TNCO is less affected. Consequently, N-S structural heterogeneities in the crust and lithospheric mantle are formed to the west of the NSGL.
In summary, a series of tectonic events other than the Paleo-Pacific subduction have prominently affected the crust and lithospheric mantle in the northeastern NCC. Unlike the oceanic lithosphere, of which the first-order structural features, such as crustal or lithospheric thickness, can be directly attributed to the sea-floor age through plate tectonics in most cases (e.g., Kawakatsu et al., 2009), the continental lithosphere exhibits much greater structural diversity typically related to the long-term evolution and history of the Earth (e.g., Begg et al., 2009;Chen, 2010;Cooper et al., 2017). Our study on the Moho depth variations in the northeastern NCC provides very important insights into the lithospheric structural heterogeneities and variable responses of continental lithosphere to complex and long-term tectonic evolution.

Moho Depth Variation and Regional Seismicity
The study region is not only structurally complex but also seismically very active (Figures 1b and 6). Projection of the seismicity in the Moho depth map shows that most earthquakes are concentrated in the regions with very large variations of Moho depths ( Figure 6). This aspect is particularly well observed in the major earthquake belts, that is, nearly E-W trending Zhangbo (Zhangjiakou-Bohai) earthquake belt, the NE-to-SW directed Shanxi earthquake belt, and the earthquake belt within the Bohai Bay Basin (Figure 6). Such a correlation between earthquake distributions and the Moho depth variations in northeastern NCC is reported for the first time and is intriguing because the seismicity and tectonics at the shallow depths may be directly related to dynamic processes in the deep lithosphere in the NCC. Further analyses are currently hindered because detailed information about the focal depths of earthquakes and the rheological properties of the crust and upper mantle is not available in this region. However, this correlation is indicative that more studies on the Moho depth, earthquake distribution and mechanisms are needed to better understand the lithospheric destruction and seismic hazard in the NCC.

Conclusions
In this study, we constructed a high-resolution Moho depth map in the northeastern NCC by employing a 2-D wave equation-based migration method of P RFs from a dense broadband seismic array. The variations in the Moho depth are heterogeneous in both E-W and N-S directions and are closely related to the long-term tectonic evolution of the NCC. In combination with other geophysical observations, petrographic and geochemical data, gravity data and regional tectonics, we identified the following tectonic implications: 1. The N-S variations in the Moho depth on the western side of the NSGL may be associated with intense crust-mantle interactions related to magmatic underplating and mantle upwelling processes in the northern margin of the central to western NCC and the preservation of cratonic crust and lithospheric mantle in the interior of the TNCO. This finding implies uneven modifications of the lithosphere to the west of the NSGL. 2. The N-S Moho depth variations beneath the Yan Mountains and the Bohai Bay Basin may reflect the imprints of earlier N-S compressional deformation along the northern margin of the NCC associated with its amalgamation with the Siberian craton, in addition to different responses between the boundary area and craton interior to the subsequent lithospheric extension since the Late Mesozoic. 3. Observations on the N-S variations in Moho depth, LAB depth, and gravity from the margin to interior of the Bohai Bay Basin can be reconciled with the presence of a more buoyant lithospheric mantle in the interior of the basin, which may indicate that the lithospheric mantle beneath the Bohai Bay Basin is laterally heterogeneous chemically and/or thermally. 4. The dominant E-W structural variations may have been associated with Paleo-Pacific plate subduction under the eastern Asia since the Mesozoic time. However, the newly observed small-scale N-S structural changes may have resulted from the structural heterogeneity of the cratonic lithosphere inherited since the formation of the NCC in the Paleoproterozoic, or spatially uneven effects on the cratonic lithosphere of subsequent thermotectonic events during the long-term evolution of the craton, or both.