Crustal Anisotropy Beneath the Trans-North China Orogen and its Adjacent Areas From Receiver Functions

As an important segment of the North China Craton, the Trans-North China Orogen (TNCO) has experienced strong tectonic deformation and magmatic activities since the Cenozoic and is characterized by significant seismicity. To understand the mechanism of the crustal deformation and seismic hazards, we determined the crustal thickness (H), Vp/Vs ratio (κ) and crustal anisotropy (the fast polarization direction φ and splitting time τ) beneath the TNCO and its adjacent areas by analyzing receiver function data recorded by a dense seismic array. The (H, κ) and (φ, τ) at a total of 309 stations were measured, respectively. The Moho depth varies from ∼30 km beneath the western margin of the Bohai bay basin to the maximum value of ∼48 km beneath the northern Lüliang Mountain, which shows the positive and negative correlations with the elevation and the Bouguer anomaly. The average φ is roughly parallel to the strikes of the faults, grabens and Mountains in this study area, whereas a rotating distribution is shown around the Datong-Hannuoba volcanic regions. Based on the φ measured from the Moho Ps and SKS/SKKS phases, we propose that the crustal deformation and seismic hazards beneath the TNCO could be due to the counterclockwise rotation of the Ordos block driven by the far-field effects of the India-Eurasian collision.


INTRODUCTION
The North China Craton (NCC) as the largest and oldest known Archean craton in China is located in the eastern margin of the Eurasian plate (Zhao et al., 2001). The NCC consists of the western NCC (WNCC) and the eastern NCC (ENCC), which are separated by a Paleoproterozoic orogen, the Trans-North China Orogen (TNCO) (Figure 1). The WNCC, mainly including the Ordos Block (OB) and its surrounding Grabens, the Hetao Graben (HTG) and the Weihe Graben (WHG), is dominated by the stable and thick lithosphere, while the ENCC underwent significant reactivation and destruction during Mesozoic and Cenozoic (Chen and Ai, 2009;Zhu et al., 2011). The TNCO, a composite unit composed of the Lüliang Mountain (LLM), the Taihang Mountain (THM) and the Fenhe Graben (FHG), was formed by the collision of the ENCC and WNCC in the Late Paleoproterozoic (∼1.8 Ga), resulting in the final amalgamation of the NCC (Zhao et al., 2001). As a transition zone of the surface topography, lithospheric thickness (Chen and Ai, 2009) and gravity anomalies (Deng et al., 2014) between the ENCC with lithospheric thinning and the WNCC remaining stable cratonic roots, the TNCO has been accompanied by complicated tectonic deformation and magmatic activity during the episodic geological evolution (Ren et al., 2002;Xu et al., 2004;Zhu et al., 2011). The relative movement of the TNCO estimated from GPS is at a rate of ∼5-11 mm/year (Shen et al., 2000;Zhao et al., 2015), which is characterized by the distribution of strong earthquakes (Xu and Ma, 1992;Gao et al., 2020) (Figure 1). The mechanism of crustal deformation is important in understanding the tectonic evolution and seismic hazards in this study area. However, the details on crustal deformation beneath the TNCO and its adjacent areas are still being discussed and debated recently (Tian and Zhao, 2013;Zhang et al., 2016;Yang et al., 2018;Schellart et al., 2019;Zheng et al., 2019;Cai et al., 2021;Chang et al., 2021). Further studies on crustal structure and anisotropy are essential for constraining the crustal deformation mechanism.
Crustal deformation can produce anisotropy on the wavelength scale of seismic waves (Nicolas and Christensen, 1987;Mainprice and Nicolas, 1989). Therefore, measurements of crustal anisotropy can also help us understand the mechanism of crustal deformation and the process of tectonic evolution. Generally, the upper crustal anisotropy is attributed to the stressinduced alignment of cracks (Crampin and Peacock, 2008), while seismic anisotropy in the mid-to-lower crust and upper mantle is more likely to be caused by strain-induced lattice preferred orientation of minerals (Meissner et al., 2002). S waves from local earthquakes located in the upper crust have been used to estimate the upper crustal anisotropy beneath the ENCC (Gao et al., 2011). In the upper mantle, seismic anisotropy is commonly determined by the splitting of the SKS/SKKS phases (Zhao et al., 2008;Chang et al., 2017;. There is a gap, obviously, between seismic anisotropy derived from local S waves and teleseismic SKS/SKKS phases to measure the mid-to-lower crustal anisotropy. Sherrington et al. (2004) pointed out that seismic anisotropy inside the mid-to-lower crust cannot be ignored when we look at the crustal deformation. The Ps phase from receiver functions, a P-to-S converted wave at the Moho, is an ideal phase to estimate seismic anisotropy within the whole crust.
We noticed that Liu and Niu (2012) developed a comprehensive technique that employed the radial (R) and transverse (T) receiver functions to compute the crustal anisotropy and Sun et al. (2012) further added a harmonic analysis to enhance the robustness of the measurements. Based on these methods, crustal anisotropy beneath the southeastern and northeastern Tibetan plateau was calculated to investigate the mechanism of the crustal deformation and thickening (Sun et al., 2012;Wang et al., 2016;Xu et al., 2018). Yang et al. (2018) and Zheng et al. (2019) also measured the splitting of the Moho Ps phases from receiver functions and observed the spatial distribution of crustal anisotropy beneath the northeastern TNCO, respectively. Results mentioned above confirmed that crustal anisotropy has been successfully estimated from the splitting of the Moho Ps phases and performed to study the crustal deformation.
To further improve the understanding of the crustal structure and deformation beneath the TNCO and its adjacent areas, we determined the Moho depth, Vp/Vs ratio and crustal anisotropy using receiver functions extracted from the teleseismic data recorded by a temporary dense array of broadband seismographs in this study. We further discussed the potential implications of the crustal structure evolution and deformation mechanism by comparing other published models (Deng et al., 2014;Shen et al., 2016;Wang et al., 2020;Chang et al., 2021;Huang et al., 2021).

Data
The ChinArray project plans to roll over mainland China with a transportable array consisting of more than 1,000 broadband seismographs with a station spacing of ∼30-40 km in order to understand the structure of Earth's interior (ChinArray-Himalaya,  (Kreemer et al., 2014).
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 753612 2011). The third phase deployment of the ChinArray project covers an area of 8°× 10°centered at the TNCO ( Figure 1). Each seismograph was equipped with a Guralp CMG-3ESP or CMG-3ESPC seismometer and a Reftek 130 data logger. In this study, we used a total of 309 stations with high quality three component waveform recordings from April 2016 to January 2019 and then selected 440 earthquakes (Supplementary Figure S1) with a magnitude greater than 5.0 and an epicentral distance between 30°and 90°to calculate the receiver functions. As shown in Supplementary Figure S1, most of the teleseismic events come from inside the western Pacific subduction zone and the Java trench. However, with the other earthquakes, the overall coverages in backazimuth and distance are reasonably good.

Methods
The true orientation of the two horizontal components is important when we use the three-component recordings of the teleseismic events to generate receiver functions (Zeng, et al., 2020). Therefore we first employed the P wave particle motions of teleseismic events to estimate the true sensor orientation for each station with the method proposed by Niu and Li (2011) before rotating the two horizontal components into the radial (R) and transverse (T) directions based on the great arc raypaths connecting the events and stations.
Following the previous studies (Vinnik, 1977;Niu and Kawakatsu, 1998;Niu et al., 2007), we further projected R and T components to the principal directions (longitudinal, P, and in-plane transverse, SV) estimated from the covariance matrix in order to minimize the P wave energy in the receiver function. We then used the "water-level" deconvolution technique (Clayton and Wiggins, 1976;Ammon, 1991) to compute receiver functions in the frequency domain. The "water level" and the corner frequency of the Gaussian low pass filter were set to be 0.01 and 1.5 Hz, respectively. After careful inspection, a total of 33,772 receiver functions from the 309 stations were selected with a station average of ∼109 receiver functions.
The crustal thickness and average Vp/Vs ratio were estimated at each station by the two-step analysis following Niu et al. (2007). We first computed the initial crustal thickness using the N th -root stacking technique (Muirhead, 1968;Kawakatsu and Niu, 1994) with the Moho Ps phase alone. By searching for the range of 20-70 km, we defined the depth with the maximum amplitude as the initial crustal thickness. Figure 2A shows the depth stacking at the station 15,819 and a clear P to S conversion peak at ∼48 km. Then a refined H-κ analysis (Zhu and Kanamori, 2000) was used to estimate the final crustal thickness, H, and Vp/Vs ratio, κ: Here K denotes the total number of receiver functions at a given station and r i (t) is the amplitude of the i th receiver function at the relative arrival times of the 0p1s (t 1 ), 2p1s (t 2 ) and 1p2s (t 3 ), following the phase notation of npms introduced by Niu et al. (2007), with respect to the direct P wave. w 1 , w 2 , and w 3 are the weights of the three time windows, which are assigned to 0.5, 0.25, and 0.25, respectively. c (κ) is a coherence index of the three phases, which is introduced by Niu et al. (2007) to reduce the trade off between H and κ. We searched for H within ±20 km of the initial depth derived from the depth stacking and κ in the FIGURE 3 | Results of crustal anisotropy from the joint analysis with multi-component receiver function data at station 15,819. (A-C) show three individual methods to estimate seismic anisotropy: 1) radial energy maximization with a cosine moveout correction; 2) radial correlation coefficient maximization; 3) transverse energy minimization. (d) is the joint solution. (φ, τ) are searched in the range of 0-360°and 0-1.5 s with an increment of 1°and 0.02 s, respectively. White plus marks the measured (φ, τ) where the objective functions reach to maximum. (E,F) show the calculated SNRs of stacked receiver function data as a function of the square root of the subsample numbers, S 1/2 . Open and filled symbols in (E,F) indicate SNR calculated from stacks of receiver functions before and after the removal of seismic anisotropy measured by the joint receiver function. More specifically, opened squares shown in (E) indicate stacks from T receiver functions without anisotropy correction but with polarity correction, while opened circles are from T receiver function with no corrections of anisotropy and polarity. Filled squares and circles represent stacks after correction of anisotropy, and with and without polarity correction, respectively.
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 753612 range of 1.5-2.0. H and κ were determined when the summed amplitude, s (H,κ), reached its maximum. The result of H-κ analysis at the station 15,819 is shown in Figure 2B, which gives an estimate of H 46.1 km and κ 1.744. For the comprehensive analysis of crustal anisotropy , we first calculated the moveout and made corrections for each station using the estimates of the above H-κ analysis so that all the cluster receiver functions were expected to have a Moho Ps arrival time equivalent to the one with an epicentral distance of 60°and a source depth of 0 km. The R and T receiver functions were further normalized with the peak amplitude of the P wave. Then, we analyzed the R and T receiver functions as a function of back azimuth to study the systematic variations in the peak Ps arrival time of the R receiver functions and polarity changes of the T receiver functions. The normalized maximum amplitude (A n,max ), maximum energy (E n,max ), and minimum total residual (R n,min ) between each receiver function and the stacked receiver function were calculated and summed along a harmonic moveout with a degree varying from 1 to 8 by the method proposed by Sun et al. (2012). In Figures 2D,E, we show the R receiver functions plotted as a function of back azimuth and the result of harmonic analysis at the station 15,819. Then, we employed three individual objective functions (IFOs) and one joint objective function (JOF), and performed a statistical analysis on the robustness to estimate the crustal anisotropy beneath each station. A more detailed description on the IFOs, JOF and the statistical analysis can be found at Liu and Niu (2012).
Here, the Figure 3

RESULTS
Based on the H-κ analysis and the crustal anisotropy estimation, we obtained 309 measurements of Moho depth and Vp/Vs ratio (Supplementary Figure S2), and 62 robust parameters of the crustal anisotropy ( Figure 5). All estimated results are listed in Supplementary Table S1, which is organized by grouping stations in the following tectonic units: the Central Asia Orogenic Belt (CAOB), the Ordos Block (OB), the Lüliang (LLM) and Taihang (THM) Mountains, the Bohai Bay Basin (BHBB), the Hetao Graben (HTG), the Weihe-Fenhe Graben (WFG) and the Datong-Hannuoba (DHV) Volcano. We also present the average crustal thickness, Vp/Vs ratio, and crustal anisotropy, together with other parameters, in Supplementary Table S2 for the further comparison.
The estimated (H, κ) at 309 stations were interpolated into 1,419 meshed grids of 0.25°× 0.25°to show the lateral variations of Moho depth ( Figure 4A) and Vp/Vs ratio ( Figure 4B) beneath the study area, respectively. The Moho depth varies from ∼30 to ∼48 km ( Figure 4A). Overall, the crust across the study area gradually tends to be thick from the southeast to northwest. Specifically, the thinnest crust, less than 32 km, is observed beneath the western margin of the BHBB while the thickest part (∼48 km) is nearby the northern LLM. The eastern HTG appears to have the highest Vp/Vs ratio up to ∼1.90, while the lowest Vp/Vs ratio, generally below 1.65, is observed between the OB and LLM. The average Vp/Vs ratio of the CAOB is relatively lower in this study area (Supplementary Table S2), but most of the study areas seem to be higher than the global average 1.74 (Kennet et al., 1995;Zhu, 2018).
We measured the crustal anisotropy at all the stations in the study area and listed the estimated parameters in the Supplementary Table  S1. However, the crust with weak anisotropy or isotropy beneath the stations does not lead to the corresponding pattern of the harmonic degree (n 2) and the SNR test ( Figure 2E and Figures 3E,F). Among the 309 stations in the study area, 7 stations did not pass the SNR test. We used "−" to indicate the measurements of 42 stations without enough back azimuthal coverage for anisotropy analysis and 7 stations with the Ps phase splitting time greater than the maximum time shift (1.5 s) in the Supplementary Table S1. Combined with the harmonic analysis, a total of 62 reliable measurements with a harmonic order of n 2 are finally shown in the Figure 5 to characterize the crustal anisotropy beneath the TNCO and its adjacent areas. In addition, we take station 15,819 as an example to exhibit the crustal anisotropy beneath the western DHV (Figure 1), which was measured with the fast polarization direction φ 22°and splitting time τ 0.26 s, respectively. The average (φ, τ) of each geological block is also listed in Supplementary Table S2. The average φ measured from most of stations is roughly parallel to the strike of major geological blocks, for example, the strike of the CAOB, LLM, THM, HTG and WFG. The average τ at most of tectonic blocks is approximately 0.3 s and the peak average value (τ 0.45 s) is located at the WFG. Particularly, stations located at the DHV and their surrounding regions (marked by white ellipse in Figure 5) are distributed in rotating fast polarization directions.
We also compared our Moho Ps fast polarization direction and splitting time at the 62 stations with those measured from the SKS/ SKKS phases (blue lines, Chang et al., 2021), the GPS velocity (green arrow, Zhao et al., 2015) and the absolute plate motion (APM, thick black vectors, Kreemer et al., 2014) in the Figure 5. From the average results of these four datasets at each tectonic unite (Supplementary Table S2), we found that the datasets of the CAOB show roughly similar directions, but certain differences at other blocks, the OR, WFG, DHV, LLM and THM. Besides, the Moho Ps splitting times of the tectonic unites are much smaller than those of SKS/SKKS phase. In other words, the SKS/SKKS splitting times are mainly originated inside the mantle beneath the TNCO and its adjacent areas.

DISCUSSION
The distributions of the estimated (H, κ) and (φ, τ) show significant lateral variations, which suggests the crustal Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 753612 structures beneath the TNCO and its adjacent areas are quite complicated. In fact, previous investigations revealed that the TNCO is dominated by the strong tectonic activities and intense seismicity (Xu and Ma, 1992;Liu and Wang, 2012;Li et al., 2015). At least ten earthquakes with magnitude no less than Ms 7.0 in Chinese history occurred beneath the TNCO. Especially, three of devastating earthquakes located in the WFG, i.e., the 1,303 Hongdong Ms 8.0, 1,556 Huaxian Ms 8.3 and 1,695 Linfen Ms 7.8 earthquake, caused serious damages and casualties (Wu and Jia, 1981;Rao et al., 2017). In this study, we utilized the obtained (H, κ) and (φ, τ) to understand the crustal deformation and seismic hazards beneath the TNCO and its adjacent areas. Moho depth and Vp/Vs ratio are key parameters in constraining the structure and bulk average composition of the crust (Zandt and Ammon, 1995;Christensen, 1996). The crustal composition is classified as low Vp/Vs ratio (κ ≤ 1.76), intermediate values (1.76< κ ≤ 1.81) and high values (κ > 1.81) as the relative abundance of quartz (Vp/Vs 1.49) and plagioclase (Vp/Vs 1.87) fluctuates (Zandt and Ammon, 1995;Wang et al., 2010). An increase of plagioclase content and a decrease of quartz content usually give rise to the increase of the Vp/Vs ratio of a rock. Generally, felsic rocks tend to low Vp/Vs ratio as compared to mafic rocks. In this study, the lower Vp/Vs ratios observed at some local regions of the CAOB and central LLM ( Figure 4B) imply that the crust has a more content of felsic rocks. The HTG with the highest Vp/Vs ratio (∼1.90) represents more mafic crustal composition here. In fact, the Vp/Vs ratios measured at most of stations are higher than the global average in this study, which is likely to be intermediate to mafic crustal composition.
In addition, the Moho depth gradually thickens from the southeast to northwest in this study area ( Figure 4A), which is consistent with the previous results from the receiver functions (Li et al., 2014;Cai et al., 2021), surface wave dispersion (Shen et al., 2016) and the Bouguer anomaly (Deng et al., 2014). According to the Airy isostasty and WGM2012 global model (Bonvalot et al., 2012), the larger anomalies of the negative gravity are usually located at the regions with the higher mountains and thicker crust. The Bouguer anomaly is negatively correlated with the elevation and Moho depth (Supplementary Figure S3A-B) in this study, which indicates that our Moho depth is reliable and robust. Here, the Airy isostatic theory was used to examine the variation of the Moho depth beneath the TNCO and its adjacent areas. Assuming the crustal density as empirical value, Wang et al. (2010) and Tugume et al. (2012) suggested that 1 km uplift of the surface topography corresponds to ∼10 km predicted crustal thickening. The positive correlation between the Moho depth and elevations is shown in Supplementary Figure S3C, but its correlation coefficient is only 0.43. The weak correlation implies that a certain amount of mantle mafic materials compensate for the crust to keep the Airy equilibrium (Ji et al., 2009;Tugume et al., 2012). Furthermore, the Vp/Vs ratio increases roughly with the crustal thinning (Supplementary Figure S3D). Based on the widespread Cenozoic basalts and pyroclastic rocks in the TNCO and its adjacent areas (Xu et al., 2004;Qian et al., 2017), the variation and correlation of our datasets argued above support the model (Ji et al., 2009;Hu et al., 2020) that the underplating of mafic magmas from the partial melting upper mantle intruded into the deep crust to balance the crustal buoyancy. The higher Vp/Vs ratio demonstrates that the crustal composition was compensated by the underplated mafic intrusion when the local crust was gradually thinning during the Cenozoic tectonic extension.
We further argue that the observed crustal anisotropy, associating with other parameters (e.g., the GPS velocity, the APM direction and the φ of the SKS/SKKS phase), reflects the Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 753612 6 complex crustal deformation beneath the TNCO and its adjacent areas. In the north of the study area, the CAOB as the largest Phanerozoic accretionary orogenic belt has grown significantly in the crust and developed a series of E-W strike slip faults since the Paleo-Asian Ocean closed (Xiao et al., 2003;Windley et al., 2007). The φ measured from Moho Ps and SKS/SKKS phases, as well as the directions of the GPS and APM, are nearly E-W beneath the CAOB ( Figure 5 and Supplementary Table S2), which suggests that the whole lithosphere deforms coherently along depth with a compressional direction of N-S. Around the Ordos block, the average φ of the Moho Ps phase trends NE in the southeast and rotates to E-W trending in the north ( Figure 5 and Supplementary Table S2). This counterclockwise rotation pattern was also observed by the φ of the SKS/ SKKS phase and speculated to be derived from the corresponding asthenospheric flow caused by the northeastward growth and expansion of the Tibetan plateau (Chang et al., 2017(Chang et al., , 2021. The ongoing collision of the Indian and Eurasian plates since ∼50 Ma has driven the asthenospheric flow (Chang et al., 2021;Wang et al., 2020) and counterclockwise rotation of the Ordos block (Zhang et al., 1998), which may contribute to generating the leftlateral shear stress in the lithosphere and observed seismic anisotropy around the Ordos block. Another noteworthy feature revealed by this study is the rotating distribution of the fast Moho Ps phase directions around the DHV regions in Figure 5. The mode of the φ around the DHV regions is also consistent with the recently studies of SKS/SKKS splitting (Chang et al., 2021) and surface wave dispersion (Huang et al., 2021). At some active volcanic areas (e.g., the Mount Ruapehu volcano, New Zealand), the anisotropic orientation agreed with our observations and was further inferred to be due to the stress-aligned microcracks caused by the magmatic eruption (Gerst and Savage, 2004;Johnson et al., 2011). The prominent low velocity anomalies from both the traveltime and full waveform tomography have confirmed the existence of the mantle upwelling beneath the DHV regions related to the stagnancy and dehydration of the Pacific slab in the mantle transition zone (Lei, 2012;Tao et al., 2018). In addition to the anisotropies from P wave tomography (Tian and Zhao, 2013) and SKS/SKKS splitting (Chang et al., 2021), we propose that the fossil crustal anisotropy with the rotating φ around the DHV regions is resulted from the asthenosphere upwelling induced by the Pacific plate subduction.
We noticed that the historical destructive earthquakes occurred in the extensional WFG ( Figure 1) where the observed average τ is the largest, 0.45 s, in this study area (Supplementary Table S2). GPS measurements, still extending with a rate of ∼4 mm/year (Shen et al., 2000), and large strike-slip rate of master-faults, 5.68-7 mm/year (Xu and Ma, 1992), further indicated the strong crustal deformation inside the WFG. By analyzing the causes of crustal anisotropy and relationship between Moho depth and Vp/Vs ratio above, we speculate that the counterclockwise rotation of the Ordos block driven by the far-field effects of the India-Eurasian collision facilitated the left-lateral shear stress and extensional crustal deformation of the WFG and in turn generated the high seismic hazards.

CONCLUSION
In this study, we investigated the Moho depth, Vp/Vs ratio and crustal anisotropy beneath the TNCO and its adjacent areas using receiver function data with the H-κ stacking method and the joint inversion scheme. The variations of the Moho depth and Vp/Vs ratio might reveal that the underplated mafic intrusion compensated for the crustal composition and Airy equilibrium. Based on the φ of the Moho Ps and SKS/SKKS phases, we proposed that the CAOB showed a coupling deformation between the crust and mantle, under the N-S compressive stress, with the closure of the Paleo-Asian Ocean. The asthenosphere upwelling induced by the Pacific plate subduction caused the stressaligned microcracks and the rotating distribution of the crustal anisotropic azimuths around the DHV regions. The observed average φ of the Moho Ps phase and the high seismic hazards are attributed to the counterclockwise rotation of the Ordos block driven by the far-field effects of the India-Eurasian collision.
FIGURE 5 | Comparison of crustal anisotropy measured from the Moho Ps (the red bar lines) with measurements of the SKS/SKKS splitting (the blue bar lines, Chang et al., 2021), the GPS velocity (the green arrows, Zhao et al., 2015) and the absolute plate motion direction of the GSRM v2.1 model (the black arrows, Kreemer et al., 2014). White solid ellipse outlines the rotating distribution of the fast polarization directions measured from the Moho Ps. The abbreviations are the same as those defined in the caption of Figure 1.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
XX estimated the receiver functions, calculated the Moho depth, Vp/Vs ratio and crustal anisotropy and wrote the manuscript. ZD provided the raw data and guided the project. LL provided many useful instruction and suggestions. FL offered the corresponding codes and proposed the main ideas of this study. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
The waveform