S wave Velocity Structure of the Crust and Upper Mantle Beneath Shanxi Rift, Central North China Craton and its Tectonic Implications

The Shanxi rift is located in the central part of the North China Craton (NCC). With strong tectonic deformation and intense seismic activity, its crust‐mantle deformation and deep structure have always attracted wide attention. Using teleseismic events observed in a dense network of 610 temporary and 127 permanent stations in the central NCC, we obtained the crust‐mantle S wave velocity structures by the joint inversion method of receiver functions and surface wave dispersion data. Our results show that the crust thickens in the northern part of the Shanxi rift (41 km) and thins in the southern part (35 km). The Taiyuan and Linfen basins, located in the central part, have high‐velocity zones in the lower crust and upper mantle; beneath the Yuncheng basin in the south, there are low‐velocity zones in the lower crust and uppermost mantle (30–80 km); the Datong basin, located in the northern part of the rift, exhibits a wide range of low‐velocity anomalies in the middle to lower crust and upper mantle. We speculate that the destruction of the NCC and associated lithospheric thinning had a significant impact on the southern part of the rift, but is still in its early stage in the central part, which retains most of the craton features. An upwelling of the asthenospheric magma occurred beneath the western part of the Datong basin. The horizontal deflection of the asthenospheric flow causes low‐velocity anomalies in its surrounding area, which is consistent with LAB topography.

rifting. Detailed velocity structures could provide important clues to this issue and may provide additional support for one of the competing tectonic models.
To probe the crust-mantle velocity structure of the NCC, numerous tomographic methods have been applied, such as surface wave tomography (Huang et al., 2009;Shen et al., 2016;Wang et al., 2012), ambient noise tomography (Fang et al., 2010), body wave tomography (Zhao et al., 2009), and receiver function inversion (Tang et al., 2010). The consensus of these studies is that there are extensive low-velocity zones in the crust and upper mantle beneath the northern part of the central NCC. However, the body wave tomography is unable to resolve detailed crustal structures, and the study using the receiver function inversion only gave the migration images for two profiles crossing the Shanxi rift. The models of surface wave tomography can achieve a horizontal grid resolution of 0.5° × 0.5°, but the resolution of deep vertical structures using longer periods of surface waves is inadequate. The resolution of these results is therefore insufficient to detect crust-mantle velocity characteristics in the Shanxi rift.
Since the P wave receiver function is sensitive to velocity discontinuities, by making full use of their complementary advantages, the joint inversion of the receiver functions and surface wave dispersion data can effectively suppress the nonuniqueness of the inversion results and obtain a more reliable S wave velocity structure (Julià et al., 2000). Recently, this method was applied to the Shanxi rift to study the S wave velocity CAI ET AL. structures . The results showed notably different features between the northern and southern part of the Shanxi rift. However, the permanent stations they used have coarsely distributed, so the ambiguity in the interpretation of the deformation and evolution of the Shanxi rift mentioned above could not be fully resolved.
In this study, we will use observations from a total of 737 permanent and temporary seismic stations (Figure 1a), which is the densest network in the study area by far. The joint inversion method of receiver functions and surface wave dispersion data will be applied to obtain the S wave velocity structure of the crust and upper mantle. With these results, we expect to shed light on the crust-mantle deformation and the current evolution state of the Shanxi rift.

Data
In this study, we used observations from 336 temporary stations of Phase-III ChinArray project (808 stations in total), with an average station distance of ∼30 km (ChinArray-Himalaya, 2011). The additional data we used here are from 127 permanent stations (SEISDMC, https://doi.org/10.11998/SeisDmc/SN), and 274 temporary stations available at the IRIS Data Management Center ( Figure 1a). We selected teleseismic events with epicentral distances between 30° and 90° and magnitudes larger than 5.5. In order to obtain reliable receiver functions, we calculated the P wave signal-to-noise ratio (SNR) in the frequency domain using the average amplitude ratio of the signal and noise in the 30 and 40 s before and after the P wave, and selected the waveforms with SNR exceeding 4.0, which resulted in a catalog of 3,863 teleseismic events ( Figure 1c).

Receiver Functions
We extracted the P wave receiver functions by using the maximum entropy deconvolution method (Tselentis, 1990;Wu et al., 2007), which uses maximum entropy as a criterion for determining the autocorrelation and cross-correlation functions, thus efficiently avoiding subjective assumptions outside of the time window and improving the resolution of the receiver functions. The original seismic recordings were processed using a bandpass filter of 0.02-1 Hz. Subsequently, the horizontal components were rotated into a radial and transverse coordinate system, and the receiver functions were obtained by deconvolving the vertical from the radial component. A Gaussian filter with an α-coefficient of 2.5 and a volume factor of 0.001 were applied to the deconvolution method. The receiver functions with clear P wave phases were selected for further study.

H-κ Stacking
The Moho depth and crustal V p /V s ratio was estimated using the H-κ stacking method (Zhu & Kanamori, 2000) with the primary converted phase at the Moho interface Ps, and the multiple reflected phases PpPs, PpSs/PsPs. At the interface with depth H, the time differences between the Ps, its multiple-wave PpPs, PpSs/PsPs phase and the direct P wave are: Here, p is the ray parameter of the P wave.   / P S V V . P V and S V are the P and S wave velocities of the interface, respectively.
We then use the amplitudes of the Ps, PpPs and PpSs/PsPs to construct the objective function S(H,κ): where   F t is the amplitude of the receiver functions; 1 t , 2 t , and 3 t are the calculated arrival time of the Ps, PpPs, and PpSs/PsPs with the interface depth H and the crustal V p /V s ratio κ, respectively;  1,  2, and  3 are the weighting factors, which are set to 0.6, 0.3, and 0.1. When H and κ are close to the real crustal thickness and V p /V s ratio, the objective function will reach its maximum. The variance is estimated from the flatness of the objective function at the maximum value:   be caused by multiple reflections from the sedimentary layers. Despite that, the crustal thickness and Vp/Vs ratio in these cases can still be obtained correctly, but with a greater uncertainty range.

Surface Wave Dispersion
The surface wave dispersion curves used in our study are obtained from Shen et al. (2016). They used ambient noise data from more than 2000 seismic stations throughout China to produce Rayleigh wave phase velocity maps with a period range of 8-50 s and a spatial resolution of about 0.5° × 0.5° (Shen et al., 2016). We used phase velocities with a period range of 10-50 s ( Figure 3) and constructed the phase velocity curves for each station by interpolating adjacent grid points.

Joint Inversion
We used the linear joint inversion method of Julià et al. (2000) to obtain the S wave velocity structure of the crust and upper mantle in the study area. To set up the initial velocity model used for the joint inversion, the crustal thickness and Vp/Vs ratio were from the H-κ stacking in this paper. The sedimentary layer was set to 5 km according to the artificial sounding profile (Jia et al., 2014). The search range for the crustal thickness was limited to ±8 km of the initial crustal thickness. The initial S wave velocities were set to 2.5 km/s, 3.5 km/s, and 4.5 km/s for the sedimentary layer, crust, and upper mantle, respectively. The mantle velocity structure is more stable compared to the crust, and the average mantle Vp/Vs ratio is about 1.80 according to the global average velocity model (PREM, Ak135), so we fixed the Vp/Vs ratio of the mantle to 1.80 in the joint inversion. The density was estimated using Birch's law ρ =  0.32 0.77 p V (Birch, 1961a(Birch, , 1961b. Since the receiver functions and surface wave dispersion curves have their own physical units and number of data points, Julià et al. (2000) defined a joint method to predict the inversion residuals. The maximum number of iterations was set to 120, and the inversion process terminated when the residuals reached a minimum. Here we chose the inversion results with the joint residuals less than 0.1, which meant the similarity of the data to the theoretical waveforms was larger than 90%.
We compared the joint inversion results of the receiver functions and the surface wave dispersion curves with different weighting factors ( Figure 4). When the weight of the receiver functions is larger than 0.7, the surface wave dispersion curves are not well fitted. If the weight is too small (0.2), which means the weight of the surface wave dispersion curve is large, the velocity structures become very smooth and does not explain the receiver functions well. In the range of 0.4-0.6, the obtained velocity structures are almost similar to each other. In order to make better use of the characteristics of both receiver functions and surface wave dispersion curve, the two data sets are weighted equally (weight = 0.5).
Since each station usually has a lot of receiver functions, the inversion process will consume a large amount of computation time if we fit all of them individually. Moreover, heterogeneities near the stations can affect the obtained receiver functions at different epicentral distances, thus decreasing their similarities. To reduce the effect of near-station structure, it is advisable to select several receiver functions with different epicentral distances. In this study, we used stacked receiver functions with different average slowness (epicentral distance) for the inversion. Specifically, we counted the number of receiver functions in each slowness interval of 0.002 s/km from 0.04 to 0.08 s/km and selected three or four intervals with the highest number of the receiver functions to be stacked separately (Figure 4a).
It is known that the receiver functions penetrating a fault is affected by the structure of the fault zone and exhibit different waveforms from the receiver functions parallel to the fault. Therefore, we use receiver functions parallel to the faults to avoid fault effects and influences from other tectonic structures. Figures 5  and 6 show two examples of the joint inversion results at station CN.HMA and S0.B037, respectively, by using the receiver functions with back azimuths in the range of 0-180°, 0-90°, and 90-180°. CN.HMA is located in the southern part of the Shanxi rift and S0.B037 is located in the Ordos block. Both stations are circled in Figure 1a. Station CN.HMA is surrounded by several NE-SW striking faults (Figure 1a). Figure 5 shows that the receiver functions in the range of 0-90° and 90-180° exhibit clear differences. In the 0-90° range, the receiver functions and surface wave dispersion curves fit better and the corresponding residual is the smallest, thus we use the inversion results in the range of 0-90°. At station S0.B037, where there is no significant fault around (Figure 1a), the residual in the range of 0-180° is the smallest (Figure 6) and the results are used in the subsequent analysis. Following this approach, we divided the receiver functions of all stations into three groups of 0-180°, 0-90°, and 90-180° according to their back azimuths to stack, and fitted the stacked receiver functions separately in the three groups. Subsequently, we selected the inversion results with the smallest fitting residuals among the three groups for further analysis.
CAI ET AL.
10.1029/2020TC006239 6 of 16 We obtained the crustal thickness by calculating the vertical gradient of the velocity curve at each station ( Figure 7). Specifically, the shallowest depth at which the S wave velocity is 4.2 km/s was searched for between the range of 0 and 80 km and used as a reference depth. The depth with the largest vertical gradient within ±8 km of the reference depth was searched for as the crustal thickness and then confirmed manually.

Crustal Thickness
Based on the H-κ stacking method, we obtained the distribution of crustal thickness and Vp/Vs ratio in the study area (Figure 8). The results show that the North China basin is underlain by a thin crust of about 32 km, and it increases to nearly 40 km in the Taihang mountain. In the Shanxi rift, the crust gradually thickens from south to north. Specifically, the crustal is approximately 35 km thick in the Yuncheng basin and thicker in the Linfen basin and Taiyuan basin at about 38 km. In the western part of the Datong basin, the crust is thicker than 41 km, and it decreases to about 38 km in the eastern part. The Ordos block is underlain by a thick crust between 42 and 44 km, which is thicker in the south and north and thinner in the center. The results are mostly consistent with previous H-κ stacking results using data only from permanent stations (Liu et al., 2011;Ge et al., 2011), gravity inversion (Zhang et al., 2011) and deep seismic sounding profiles (Figures 9b and 9c, Jia et al., 2014). Figure 8b shows the distribution of Vp/Vs ratio. The average Vp/Vs ratio is relatively low in the central part of the Ordos block, but higher in its southern and northern parts. The Vp/Vs ratio in the Taihang mountain is about 1.73. The Shanxi rift shows a high Vp/Vs ratio in the south and north, and relatively lower in the central part, which is generally in agreement with the results of the previous study (Ge et al., 2011 Figure 9a shows the crustal thickness beneath the study area by using the joint inversion of receiver functions and surface wave dispersion curve. Comparing the crustal thickness obtained by the H-κ stacking method and joint inversion method (Figures 8a and 9a, 9b), we can find that the overall features of the two results are similar. Both results are characterized by a thicker crustal thickness in the western region and the thickness of the crust increases from south to north along the Shanxi rift. The difference is that the crustal thickness in the northeastern Ordos block and Yinshan mountain from the joint inversion is a little thicker than that from the H-κ stacking results (Figure 9b). Comparison with the results of a deep seismic sounding profile in the study area (Jia et al., 2014) shows that the crustal thickness obtained by the two methods do not differ significantly from the results from the artificial seismic profile analysis (less than 2 km, Figures 9c  and 9d), indicating that the crustal thickness obtained in this paper is reliable. Compared with previous results using only permanent stations (Ge et al., 2011;Liu et al., 2011), the application of temporary stations has a better constraint on the crustal thickness in the basins and provides a more reliable reference for studying the crust-mantle structure in central North China Craton.

S wave Velocity Structures
The distribution of the S wave velocity structures at different depths is shown in Figure 10. The structural heterogeneities are visible and will be discussed in detail below.  thinner sedimentary layer of about 0-2 km (Jia & Zhang, 2005). Therefore, the velocity difference in the shallow layers is consistent with changes in sedimentary thickness.
At the depth of 13 km, the Shanxi rift, North China basin and Hetao basin show relatively high-velocity anomalies. The central part of the Ordos block is characterized by high velocities, while its northern and southern parts show relatively low velocities. This is consistent with its thin crustal thickness in the center and thicker crust in the north and south, indicating that its upper crust may have a similar structure to that of the Moho. In the northeastern part of the Ordos block, the upper crust presents high-velocity anomalies. As shown in the velocity profile of Ai et al. (2019), the high velocity at this depth is also confirmed by the presence of uplift of the upper and middle crust in this region.
At the depth of 26 km, the central part of the Shanxi rift exhibits high velocities. The Datong basin and Yinshan mountain have large-scale low-velocity anomalies, which is similar to the velocity structure of Ai et al. (2019) in the lower crust. Furthermore, Tang, Chen, et al. (2013) combined ambient noise with surface waves and showed inversion results of similar characteristics at depths of 20km-Moho. One improvement, however, is that our results show low-velocity anomalies around the Yuncheng basin, which can clearly reveals its exceptional structure.   also reflects the deformation of the lower crust is to some extent causing large areas of crust thickening and topography uplift.
At the depth of 70 km, the Taiyuan basin and Linfen basin, located in the central part of the Shanxi rift, as well as the Ordos block, present high velocities, showing the features of their upper mantle, but the velocity is relatively low in the south of the rift. Large areas of low-velocity anomalies span from the northeast of the Ordos block to the Datong basin, and the North China basin. These low-velocity features, which have also been revealed in previous studies Tang, Chen, et al., 2013), suggest the possibility of thermal upwelling of the upper mantle.
At the depth of 105 km, the Ordos block still appears as a rigid block with high velocities, and the southern part of the Shanxi rift presents a high-velocity feature comparable to that of the central part. The Datong basin and North China basin still exhibit extensive low-velocity features, which may reflect the thermal effects from the deep mantle structure.

The Effect of Destruction and Lithospheric Thinning in the Central and Southern Parts
As can be seen in Figures 10 and 11, the crust and upper mantle of the Ordos block present high velocities, showing a stable craton basin. In the central part of the Shanxi rift, the lower crust and upper mantle display high velocities comparable to those of the Ordos block (Figures 11b and 11d), and the high velocities extend eastward to 112.5°E. From the imaging results of Tang, Chen, et al. (2013), cratons are ancient cratons with cold and thick lithospheric mantles. This area still has the characteristics of the Archean craton, as indicated by the distribution of the high-velocity Archean rocks. Different from the northern and southern part of the rift, it has a thick lithosphere of 110-120 km (Tang, Chen, et al., 2013;Zhang et al., 2019) and low temperature in the crust and upper mantle, only slightly warmer than that in the stable Ordos block . We speculate that the current destruction and lithospheric thinning in the central part is not very clear or in an early stage, so the rocks have not been weakened and retain most of the craton features, which supports the second geological model introduced above, that the rifting CAI ET AL.
10.1029/2020TC006239 11 of 16 process of the Shanxi rift developed from both its southern and northern extremities toward the central part (Xu et al., 1993).
In the southern part of the Shanxi rift, velocities in the lower crust are low but relatively high below 80 km (Figures 10 and 11c). However, the images of Tang, Chen, et al. (2013) showed that the southern part has widespread low-velocity anomalies down to depth 190 km. This may be due to the insufficient resolutions of their results to illustrate the thin high-velocity layer. The thickness of the lithosphere is only about 85 km (Tang, Chen, et al., 2013), showing a thinning lithosphere. The temperature from the crust to the upper mantle is also higher in the south than that in the central part . The low velocities in the crust-mantle with thin lithosphere and high temperature indicate that the destruction has a significant CAI ET AL.
10.1029/2020TC006239 12 of 16 effect on the southern part of the rift. The thin high-velocity layer at 80-105 km may be a remnant craton structure, which has not been completely destructed yet. The high-Mg xenoliths at Hebi county, located in the southern Shanxi rift, also attested to the presence of the Archean lithospheric remains at relatively shallow depths (Zheng et al., 2001).

Asthenospheric Magma Flow in the Northern Rift
The middle and lower crust to the upper mantle of the Datong basin exhibits a wide range of low velocities ( Figure 10). The profile AA′ and DD′ (Figures 11a and 11d) also show that there is a low-velocity layer in the middle and lower crust and a bottom-up low-velocity zone in the upper mantle. Although the results obtained by Ai et al. (2018) using permanent stations also show a low-velocity zone in the upper mantle, they did not clearly reveal the low-velocity layer in the middle-lower crust. The thermal structure research shows that the Datong volcano area has a high temperature from the surface down to 300 km depth . In addition, there are high-conductivity layers spanning from the middle and lower crust to the upper mantle (Zhang et al., 2016). These characteristics indicate the existence of the hot materials with low velocity and high temperature beneath the Datong basin. It was found that Cenozoic basalts are widely distributed in the Datong basin (Zhang et al., 2003). Chemical elements analysis of the rocks proved that these basalts mainly came from the deep asthenosphere, indicating that there may be an upwelling channel in the asthenosphere (Xu et al., 2005). The teleseismic tomography results show that the low-velocity anomaly exists down to 400 km depth (Huang & Zhao, 2006;Lei et al., 2012), which suggests that the asthenospheric upwelling is possibly caused by the dehydration of the stagnant Pacific slab (Huang & Zhao, 2006;Li & van der Hilst, 2010). However, the high-velocity anomaly representing the stagnant Pacific slab is far away from the Datong basin, so it is argued that the upwelling may represent a mantle plume from the lower mantle (Lei, 2012). Although no unified conclusion on the deep source has been developed, it's reasonable to assume that the lithospheric destruction in the Datong basin is related to the asthenospheric upwelling. Under the far-field effect from the continuous NE compression of the Tibetan Plateau and slab rollback related to the subduction of the Pacific plate (Ren et al., 2002;Xu & Ma, 1992), the Shanxi rift is in a maximum NW-SE extensional stress field (Li & Van Der Hilst, 2010;Qu et al., 2014). The SKS splitting results illustrate that the fast polarization axis of the SKS is predominantly NW-SE in the Datong basin (Huang et al., 2011;Zhao et al., 2011), which is explained by the horizontal deflection of the asthenospheric flow as it rises to the LAB . The lithosphere is less than 110 km in the Datong volcano area but thicker in the surrounding regions (Tang, Chen, et al., 2013;Zhang et al., 2019). Combining the stress system with the anisotropy and velocity structures, we could further extrapolate that when the asthenospheric flow rises to the LAB, it is deflected horizontally and centrifugally, mainly in the same direction as the NW-SE extension ( Figure 12). We also found that a relatively small velocity gradient zone appears at the crust-mantle boundary beneath the Datong basin (Figures 11a), indicating that the differentiation of crust-mantle material is not complete and a clear and sharp Moho interface has not yet been formed.

Conclusions
We applied the joint inversion of receiver functions and surface wave dispersion curve to investigate the crustal thickness and S wave velocity structure of the crust and upper mantle beneath the study area (34°N-42°N, 108°E−115°E). With the help of the data from the densest network of temporary and permanent stations in the central NCC to date, the results have a horizontal resolution of approximately 30 km, comparable to the density of the station distribution. Combined with the geological and other geophysical data, the deformation mechanism of the Shanxi rift is explored, providing a new reference for understanding the evolution process of the Shanxi rift.
1. In the central part of the Shanxi rift, the lower crust and upper mantle show high velocities, comparable to that of the Ordos block. Its lithosphere does not show clear thinning and intense deformation, so we speculate that the destruction in the central part is at an early stage and it retains most of its craton features 2. The lower crust and upper mantle in the south of the Shanxi rift generally shows low velocities, but with a thin high-velocity layer below 80 km. It has a thin lithosphere and high temperature, implying that the destruction has a significant but not complete effect on the southern part of the rift and that the thin high-velocity layer maybe a relic of the Archean lithosphere 3. Datong basin and its surrounding regions exhibit large-scale low-velocity anomalies in the middle-lower crust and upper mantle, which may be related to the upwelling erosion of the asthenospheric hot materials. The upwelling center may be located in the west of the Datong basin. Under the resistance of the LAB, the asthenospheric flow may deflect horizontally, resulting in the low velocities in the surrounding regions of the Datong volcano area