Seismic Evidence on Different Rifting Mechanisms in Southern and Northern Segments of the Fenhe‐Weihe Rift Zone

The Fenhe‐Weihe Rift (FWR) system is one of the most active rift zones in east China, and plays an important role in the reactivation of the North China Craton. In this work, a high‐resolution 3‐D lithospheric S wave velocity model of the FWR is built by joint inversion of receiver functions and ambient noise. Strong correlation between topography and tectonic features is observed at shallow depths. At great depths, strong heterogeneities are observed along the FWR: (1) high‐velocity anomalies emerge beneath the Weihe and Taiyuan basins in the lower crust; (2) two high‐velocity anomalies are imaged in the uppermost mantle beneath the south to central FWR; and (3) two strong large‐scale, low‐velocity zones are observed beneath the north Trans‐North China Oregon: one is at ~60 km and the other one is much larger with a deeper source. Based on the petrologic properties, geochemical analyses, and seismic images obtained in this work, we propose that the rifting mechanisms between the south and north FWR are different. In the south, the rifting could be mainly caused by the passive effects from the uplift of the Tibetan Plateau, while in the north the rifting is due to the combined effects of the asthenosphere upwelling and the counterclockwise rotation of the Ordos block. The uppermost mantle high‐velocity anomalies underneath the south to central part of the FWR could be the remnants of strong ancient cratonic root, and act as mechanical barriers to block the rifting processes between the south and north FWR.


Introduction
The Cenozoic Fenhe-Weihe Rift (FWR) consists of a series of asymmetric basins located at the east boundary of the Ordos block in North China, including the Datong, Xingding, Taiyuan, Linfen, Yuncheng, and Weihe basins from north to south (Figure 1). Geologic investigations showed that the FWR and its neighboring Shanxi highlands split the North China Craton (NCC) into three tectonic parts: (a) the North China Plain in the eastern NCC, (b) the Ordos Block in the western NCC, and (c) the Trans-North China Oregon (TNCO) that is composed by the FWR and the Shanxi highlands (e.g., Chen et al., 2009Chen et al., , 2014; Y. C. . The western NCC has been a tectonically stable block ever since the Precambrian amalgamation of eastern NCC and western NCC (Zhao et al., 2005). On the contrary, the eastern NCC has experienced substantial and complex tectonic deformations during the Phanerozoic era, resulting in a rather thin and hot lithosphere beneath the eastern NCC (e.g., Menzies et al., 2007;Zheng et al., 2007). As a transition zone between the eastern NCC and the western NCC, the FWR accommodates dramatically differential movements between the eastern NCC and western NCC and plays an important role in the rejuvenation of the NCC. GPS measurements revealed that the north part of the rift system is still under extension near N105°E at a rate of 4 ± 2 mm/year (Shen et al., 2000), which leads to highly active seismicity in the FWR. Therefore, understanding the rifting mechanism of the FWR is important in studying the evolution of the NCC and evaluating the seismic hazards in Shanxi Province.
To date, the rifting mechanism of the FWR is still in debate. Some researchers proposed that the evolution of the FWR is governed by the regional stress fields caused by the collision of the Indian and Eurasian plates (Xu & Ma, 1992), which drove the counterclockwise rotation of the Ordos block (Tapponnier & Molnar, 1976). However, increasing work found there are many different tectonic properties in the FWR from north to south. Bouguer gravity anomalies and seismic sounding lines revealed a 2-to 3-km uplift of Moho discontinuity beneath the central FWR and 5 to 6 km beneath the southwestern FWR (Xu & Ma, 1992). Moreover, most earthquakes in the south FWR occurred at 5-to 25-km depth, while the focal depths of the earthquakes range from 5 to 15 km beneath the Taiyuan basin (Song et al., 2012; Figure 1). This seismicity pattern demonstrates that the focal depths of the earthquakes in the FWR get shallower from south to north. Meanwhile, magnetotelluric images exhibited a significant high-conductivity zone in the southern part of the FWR, beneath the eastern margin of the Linfen Basin. This high-conductivity zone is deeply rooted and may facilitate the evolution of the rift zone under the far-field effects of the India-Eurasian Collision (Yin et al., 2017). A high-conductivity zone is also imaged beneath the Datong volcanic zone in the northern FWR (Zhang et al., 2016) and is shallower than the high-conductivity body beneath the Linfen Basin. These two separated high-conductivity zones and the unevenly distributed seismicity indicate significant lateral rheological heterogeneities along the FWR. The seismicity pattern along the Fenhe-Weihe Rift in right-side subplot is modified from Song et al. (2012). In the lower right corner, the location of the study region is plotted. Stations WAT, BOD, and ZOZT presented in Figure 11 as examples are marked with yellow, green, and red squares, respectively.
The difference of geophysical and geological characteristics between the north and south segments of the FWR implies that the north and south FWR might be experiencing different rifting processes. To fully understand the rifting mechanism, a high-resolution 3-D shear wave velocity structure of crust and uppermost mantle of the FWR is desired. Recently, several tomographic studies have been conducted in the FWR and its surroundings (e.g., Ge et al., 2011;Huang et al., 2009;Jiang et al., 2013;Y. C. Tang et al., 2013;Zheng et al., 2011). For example,  build a 3-D thermochemical model of the NCC by a multiobservable probabilistic inversion method and demonstrate the extent and possible causes of cratonic lithospheric modification. Jiang et al. (2013) present a lithospheric model of the central and western NCC from teleseismic Rayleigh wave tomography and conclude that the recent Cenozoic lithospheric reworking of the TNCO is driven by the combined effects of Pacific slab subduction and India-Eurasia collision. Beside these 3-D lithospheric models, high-resolution transects across the TNCO are also obtained from receiver functions (e.g., Chen, 2009;Chen et al., 2009;Xu et al., 2014). Specifically, Chen (2009) reveals that the eastern and central NCCs are under different lithospheric tectonics from both S and P receiver function migration. However, studies mentioned above focused mainly on either large-scale lithospheric structure of the NCC or the Moho and lithospheric discontinuity variations beneath certain seismic lines across the FWR. Therefore, the resolutions of these models are relatively coarse in the FWR system and cannot resolve the detailed features in the interior of the FWR. Shen et al. (2013a) developed a Bayesian Monte Carlo approach for jointly inverting surface wave dispersions and receiver functions for shear wave velocities. This method is particularly suitable to study the crustal and uppermost mantle structures because it employs the short-to medium-period ambient noise dispersions and P to S converting phase of receiver functions. Moreover, Monte Carlo sampling in the model space can provide us the uncertainty of a certain model parameter, which is of great importance to estimate the confidence level of the results. In this study, by taking advantage of this approach and the dense broadband seismic array deployed by China Earthquake Administration (CEA; Zheng et al., 2010), we build a high-resolution 3-D shear wave velocity model of the FWR. The S wave velocity model and associated uncertainties can further improve our understanding on the possible rifting mechanism of the FWR to a great extent.

Ambient Noise Tomography
We collect one-year continuous seismic waveforms of vertical component recorded at 212 stations from CEArray (SEISDMC; doi:10.11998/SeisDmc/SN) in 2015. All seismic stations used in ambient noise tomography are located in the NCC, especially around the FWR (Figure 1).
Following the methodology of Bensen et al. (2007) and Lin et al. (2008), empirical Green's functions (EGF) are extracted from cross correlation of each station pair. First, seismic records are decimated to 1 Hz and are cut on daily basis. Afterward, we remove the instrumental response and filter the data at 3-150 s. Subsequently, a weighted running average method is used to normalize the time series in time domain to reduce the disturbances of earthquake signals and instrumental irregularities. Finally, a spectral whitening is applied to broaden the bandwidth of the ambient noise and smooth the imbalance spectrum that may be caused by persistent localized microseisms (Bensen et al., 2007;Zeng et al, 2010;Zheng et al., 2011]. After computing cross-correlation functions for all possible station pairs, we stack the daily cross-correlation functions for a whole year. Figure 2a shows the final stacked cross-correlation functions between station pairs in the Inner Mongolia Province. Persistent localized sources of microseisms may produce disturbances on EGFs, which often arise as precursory signals in addition to the Rayleigh wave signals . Zeng and Ni (2010) identified a persistent microseismic source from the Aso volcano in Kyshu Island. The energy radiated from the Aso volcano can be clearly identified in our EGFs and caused serious disturbance (Figure 2b). Before picking out dispersion curves for each station pair, the disturbance from the Kyushu microseism must be removed. Following Zheng et al.'s strategy (2011), we minimize the spurious signals on EGFs at periods less than 18 s. If the arrival time of the disturbance signal on EGF is close to that of the surface wave signal, we keep only the opposite lag of the EGF that is free from the disturbance. On the other hand, if the two signals can be clearly distinguished from each other, we just remove the Aso signal off from the EGF and retain the symmetric component to improve the signal-to-noise ratio (SNR). For periods longer than 18 s, we retain the symmetric component of EGF since the disturbance of Aso signal is weak at longer periods.
We then measure dispersion curves from the EGFs for both group and phase velocities simultaneously using a frequency-time analysis method . We track phase velocity dispersion curves from long periods to short periods constrained by a reference global phase velocity model (Ekström, 2011) to unwrap the intrinsic 2π phase ambiguity of the phase measurements. Two criteria are adopted to ensure the reliability of the dispersion measurements: (1) SNR of the EGFs should be greater than 10 and (2) interstation distance should be longer than one wavelength. We relax the interstation distance limitation from three to one wavelength to include more measurements at long periods (Luo et al., 2015). Following Barmin et al. (2001), we construct dispersion maps for both group and phase velocities on a 0.25°× 0.25°grid at periods from 6 to 45 s. The spatial resolution of the velocity maps is defined as twice the standard deviation of a 2-D Gaussian fit to the resolution topography at each grid (Levshin et al., 2005). We perform the surface wave inversion twice. In the first inversion, the dispersion measurements with misfits larger than 6 and 3 s are discarded for group and phase velocity inversions, respectively. In the second inversion, the remaining dispersion measurements are inverted to obtain the final dispersion maps. The smoothing and damping in surface wave tomography here are chosen empirically to ensure the consistency between different dispersion maps. Generally, α, σ, and β in Barmin's method are set to be around 600, 60, and 30 in this work.

Receiver Functions
Three component teleseismic records are collected from Ge et al. (2011). The stations and teleseismic event distributions are shown in Figure 1. In extracting the receiver functions, a low-pass Gaussian filter, G(w) = exp (−w 2 /(4a 2 )), is applied to clean up the high-frequency noise. In this work, we set a to be 2.5. With a being 2.5, G(1.2) = 0.1, for example, the frequency (Hz) at which G = 0.1 is 1.2 Hz. The raw azimuthal dependent receiver functions, which are calculated with a time domain iterative deconvolution technique (Ammon, 1991;Ligorría & Ammon, 1999), are normalized to a reference slowness of 0.06 s/deg. Because of the different time and amplitude dependence on incidence angle for reverberated phases, we retain only the receiver functions within 9-s window after the direct P signal. Following Shen et al. (2013a), three criteria are used to pick out reliable receiver functions. First, we retain receiver functions whose variance reductions in the iterative deconvolution are greater than 80%. Second, we discard receiver functions with either too large (greater than 1) or too small (smaller than 0) amplitudes at zero time. Finally, receiver functions are selected with misfits less than 0.08. An example of the quality control process for station HEYON located in Hebei Province is shown in Figure 3. All 87 untreated receiver functions are arranged by azimuth in Figure 3a. After applying the first two criteria, 69 residual receiver functions are retained, as shown in Figure 3b. Actually, apparent disagreements among the original receiver functions are wiped out. In the third step, 67 azimuthally smooth receiver functions are selected for the final harmonic stripping, as presented in Figure 3c. To estimate the receiver functions for a horizontally layered isotropic equivalent medium, "harmonic stripping" is applied to extract azimuthally independent components from all the quality controlled receiver functions (Shen et al., 2013a). The results of harmonic stripping for station HEYON are shown in Figure 4. The three harmonic components are shown in Figures 4b and 4c, and Figure 4e displays the synthetic receiver functions H(θ, t) which are dependent on incident azimuths.

Joint Inversions
To bridge the Moho, a Bayesian Monte Carlo approach of joint inversion of dispersion curve and receiver function (Shen et al., 2013a) is applied. Following Shen et al. (2013a), we parameterize the 3-D shear wave velocity model as a combination of a set of 1-D depth-dependent models. For each grid, the 1-D model is composed by a sedimentary layer, a crustal layer, and an uppermost mantle layer. The sedimentary layer here is parameterized with linear gradient velocities defined by the layer thickness and velocities at the top and bottom of the layer. The crustal and mantle layers are characterized by four and five cubic Bsplines, respectively. The crust thickness is also a free parameter during the inversion but the mantle layer thickness is updated automatically as the total thickness of the 1-D model is fixed to 120 km. Totally, 13 parameters are used for model parameterization here.
We define the model space as a reference model subject to allowed perturbations and constraints. The reference model is taken from Shapiro and Ritzwoller (2002) but with sedimentary layer replaced with that taken from CRUST 1.0 (https://igppweb.ucsd.edu/~gabi/rem.html). Besides, the sedimentary thickness of the reference model is set to 1 km uniformly. In this work, the allowed constraints and perturbations from the reference model are the same as that defined by Shen et al. (2013b). But at stations that the receiver functions are not fitted well during inversions, the range of the sedimentary thickness perturbation is enlarged to be ±500%.We set V p /V s ratio to 2.0 in the sedimentary layer, and 1.75 in the crustal and mantle layers empirically. Density in the sedimentary and crustal layers is based on the scaling relationship from Christensen and Mooney (1995) and Brocher (2005). In the mantle, density is determined from the scaling relationship from Karato and Wu (1993). Attenuation factor Q is taken from PREM model (Dziewonski & Anderson, 1981) to make physical dispersion correction. It is reported that the biggest deviation of Q from PREM results in 1% shear wave velocity variations within the 1σ uncertainties (Shen et al., 2013b). Thus, Q from PREM in the mass can be used in the present inversion.
The joint inversion misfit (S joint ) for a given station is Here C i (m) and U j (m) are the calculated phase and group velocity from model m at periods i and j, C obs i , and U obs j are the corresponding phase and group velocities from ambient noise tomography. The time series A 0 (t k ) is the azimuthal independent receiver function from harmonic stripping, and R k (m) is the calculated receiver function. The parameters α, β, and s are the associated uncertainties. A factor κ is used to scale the weights between dispersion curves and receiver functions. Empirically, κ ranges from 1 to 40 (Shen et al., 2013a). In this work, both dispersion curves and receiver functions are equally reliable; thus, we set κ to 2.5 in the joint inversions following Shen et al. (2013a). For nodes without receiver function data, the misfit function consists of only S ph and S gp .
Data uncertainties are the direct input for shear wave velocity inversion. In this study, the harmonic stripping estimates the uncertainties of the azimuthal independent receiver functions, while the method of Barmin et al. (2001) for surface wave tomography produces no uncertainties but gives the spatial resolution of the dispersion maps. In this situation, we estimate the dispersion uncertainties empirically following  by scaling resolution to uncertainty via the following relationship: Here σ(r) is the position-dependent uncertainty that we need, R(r) is the resolution derived from Barmin's method, and σ 0 is the reference uncertainty at reference resolution R 0 . We set R 0 to 50 km and σ 0 as the mean value of the uncertainties in east China from a reference model , in which the uncertainties are scaled from eikonal tomography and resolution.
In this study, the joint inversion is performed for each seismic station, and surface wave dispersion inversion alone is performed for each node on the 0.25°× 0.25°grid on which we perform ambient noise tomography.
To construct an integrated model in the vicinity of FWR, a simple interpolation procedure is employed for integration of the results. At each grid, the V s value for a given depth is a weighted average of the values derived from the nearby joint inversions and the surface wave inversion alone at the current grid. The weights for values at nearby seismic stations are defined as where d i represents the distance from the current grid to seismic station i and σ i is the related uncertainty from the Monte Carlo inversion. The weight for the grid where we perform surface wave inversion alone is fixed to be a constant 0.01, a so small value to fully assimilate joint inversion results in our 3-D model. The threshold for d i is chosen to be 0.25 to ensure that there is no more smoothing effects on our results. The 3-D model is eventually constructed on a 0.25°× 0.25°spatial grid by performing this interpolation approach.

Results
3.1. The 2-D Surface Wave Tomography 3.1.1. Assessment on the 2-D Surface Wave Tomography Resolution maps at periods from 6 to 30 s for phase velocities are presented in Figure 5. The resolution of ambient noise tomography is defined as twice the standard deviation of a 2-D Gaussian fit to the resolution topography at each inversion node (Levshin et al., 2005). The resolution maps display that almost the whole study region has a resolution higher than 100 km. Especially, along the FWR, the resolution is higher than 40 km, indicating that our results are capable of resolving small-scale structures of the FWR.
We also present the probability distributions of misfits to observations from the final estimated Rayleigh wave phase dispersion maps in Figure 6. Means of the distributions are nearly zero, indicating that there is no systematical bias in the ambient noise tomography. The best fit occurs at period of 20 s, in which the standard deviation of the travel time misfits reaches to a least value of 0.59 s. But the quality of the inversion degrades at shorter periods or longer periods. Especially at 6 s, the travel time misfits on average are 1.19 s.

Journal of Geophysical Research: Solid Earth
This degradation at short periods may be due to the less reliability of the extracted empirical Green's functions at short periods as the interstation distance is on average 50 km in the study region. At longer periods, the degradation could be ascribed to the increased random errors in phase velocity measurements from signal with longer wavelength. Overall, our measurements from empirical Green's functions are coherence with one another and data misfits from the ambient noise tomography are satisfactory. Figures 7 and 8 show phase and group velocity perturbations at periods from 6 to 30 s, respectively. In general, group velocities are more sensitive to shallower structures than phase velocities (e.g., Lin et al., 2007); thus, the inversion involving group velocities can better resolve the shallower crustal structures. As presented in Figure 7, phase velocities at periods from 6 and 10 s correlate well with the surface topography. Typically, high velocities are found at the Taihang and Lvliang uplifts, while low velocities at the rifted basins are observed, especially at the Huabei basin. An interesting velocity variation pattern with velocities increasing from south to north in the FWR is observed. At 14-s period, the Cenozoic FWR is imaged as an S-shaped low-velocity belt surrounded by high-velocity anomalies in the Lvliang uplift to the west and the Taihang uplift to the east. At longer periods, the velocity maps are affected mainly by the deeper crust and the uppermost mantle. In overall, phase velocities at periods from 18 to 30 s show a similar pattern: the phase velocities in the southern FWR are notably higher than that in the northern FWR. Meanwhile, the group velocity maps at periods from 18 to 30 s exhibit stronger lateral variations. Especially at 30-s period, the group velocities in the Ordos block are much lower than that in the Huabei basin, indicating that the Moho beneath the Ordos block is significantly deeper than that beneath the Huabei basin.  Figure 9 as an example. The surface waves resolve well the velocity in the middle crust (Figure 9a), but are less sensitive to the velocity discontinuity (Figures 9c and 9e). By introducing receiver functions in the Monte Carlo inversion, the crustal thickness and the V sv jump across the Moho are better resolved than that from surface wave inversion alone (Figures 9d and 9f). The 1-D V s models at station HEYON obtained from both joint inversion and surface wave inversion alone are presented in Figure 10. The width of uncertainty corridor of velocity profile based on the joint inversion is obviously thinner than that obtained from surface inversion alone. Overall, involving the receiver functions in the inversion significantly reduces the trade-off between velocities in the lower crust and uppermost mantle. Figure 11 exhibits three depth-dependent models located in Heibei, Shanxi, and Shaanxi Provinces, respectively, from joint inversions along with the observed and synthetic receiver functions and dispersion curves. Generally, the phase velocities are fitted well and have a reduced chi-square value less than unity (χ 2 ph <1).

Journal of Geophysical Research: Solid Earth
The group velocities are not fitted well at periods shorter than 15 s, but it is overall acceptable considering that the group velocities measured at short periods are easier to be biased (e.g., Zheng et al., 2011). In the receiver functions, the Moho conversion phases are also fitted well, and the reduced chi-square value for receiver functions are around unity. Overall, fit to the data is satisfactory in the joint inversion.
One standard deviation of the posterior distribution output from Monte Carlo sampling defines the model uncertainty. Shear wave velocity uncertainty maps are presented in Figures 12 and 13. The uncertainty maps help evaluate the reliability of the features resolved from the Monte Carlo inversion. Due to the sparse distribution of seismic stations, the joint inversion strategy cannot constrain the Moho depth satisfactorily in the entire study region. Meanwhile, the uncertainties at great depths are also high because the relative short-period dispersions (<40 s) cannot constrain the deeper structures well. Generally, velocities of the crust are resolved well with uncertainties less than 0.04 km/s except regions near the Moho and somewhere near the surface. Similarly, velocities in mantle are with uncertainties less than 0.1 km/s apart from regions below the Moho immediately. Fortunately, the features mentioned in the following paragraphs are all in regions with small uncertainties, thus making it possible to analyze the small-scale features observed in the FWR.

Horizontal Shear Wave Velocity Maps
In the upper crust as presented in Figure 14a, the Taihang and Lvliang mountains are imaged as strong highvelocity zones. By contrast, the Hetao and Weihe rifts, together with the Huabei basin, appear as significant low-velocity zones (LVZ). These LVZs are consistent with the thick sediments in these basins, as is reported that the early Tertiary extension of the NCC first occurred in the Huabei, Hetao, and Weihe basins (Zhang et al., 1998). Similar to the results presented in the dispersion maps, the S wave velocities in the Taiyuan and Datong basins also exhibit low values, but are much higher than that in the Weihe basin at the southern FWR. It is well known that the Cenozoic FWR extended from south to north, resulting in a cluster of Sshaped trans-tensional basins (Zhang et al., 1998). As a result, the southern FWR may have experienced a longer and more complex rift process in history, having a weaker crust compared to the northern FWR. In the middle crust as presented in Figures 14b-14d, S wave velocities of the western Ordos are slower than that of its eastern counterpart. Geological investigations showed that the Ordos basin was developed during the middle-late Triassic and early Cretaceous, and has gone through uplifting in its eastern part and subsiding in its western part during the Mesozoic and early Cenozoic (Liu et al., 2006). Consequently, the western Ordos basin possesses a thicker sedimentary layer than its eastern counterpart, causing the lower velocities in the western part. At 20-km depth, the FWR appears as an entire LVZ belt between the Taihang and Lvliang uplifts, except for the Shilingguan and Lingshi blocks which are located to the north and south of the Taiyuan basin, respectively. This velocity pattern indicates that the regions with the smallest extension in FWR are located at Shilingguan and Lingshi blocks, which is consistent with the tectonic features observed on the surface. Another intriguing phenomenon in the upper to middle crust is that the LVZ beneath the Taiyuan basin migrates eastward with increasing depths (Figures 14a-14c). It may indicate that the weakest part of the Taiyuan basin is located in its eastern part. In the lower crust as presented in Figures 14e and 14f, velocities beneath the Taiyuan and Weihe basins are higher than the rest of FWR. The corresponding model uncertainties at these depths are low (<0.05 km/s), guaranteeing the reliability of this feature. At 40-km depth, the southern FWR is dominated by high velocities (>4.2 km/s), while the Ordos block and the northern FWR exhibit significantly low velocities (<3.9 km/ s). Considering that 40-km depth is close to the average Moho depth in this region, the Moho beneath the southern FWR may be shallower than 40 km. From the lower crust to the uppermost mantle (Figures 14g-14i), the Datong volcanic zone appears as strong low-velocity anomalies, although the center of lowest velocities varies as depth changes. As expected, the Ordos block is imaged as a rather highvelocity zone in the upmost mantle, reflecting its stable cratonic characteristics.
The lowermost crust structure, V s jump across Moho discontinuity, crust thickness, and the associated uncertainties are presented in Figure 15. High-velocity anomalies can be clearly identified beneath the Weihe and Taiyuan basins (Figure 15a). Correspondingly, the Weihe and Taiyuan basins display small-

10.1029/2018JB016476
Journal of Geophysical Research: Solid Earth velocity changes across the Moho. As expected, the crustal thickness of the Ordos block is thicker than that of the Huabei basin, suggesting the stable properties of the Ordos block and lithospheric extension of the Huabei basin (Menzies et al., 2007). Meanwhile, the crust of the southern FWR is thinner than that of the northern FWR. In the southern segment, the crustal thickness is 35-38 km, while in the northern part the crust is thicker than 42 km, suggesting nonsynchronous rifting processes between the southern and northern parts of the FWR.

Vertical Shear Wave Velocity Profiles
Cross sections aid to better examine the features of the crust and uppermost mantle. Here five cross sections are presented in Figure 16. Significantly low velocities are observed beneath the Datong volcanic zone and can be tracked down to depths greater than 60 km from cross sections AA′, DD′, and EE′, which coincides with the high-conductivity bodies resolved from magnetotelluric result in this region (Zhang et al., 2016). The low-velocity and high-conductivity anomalies suggest that a hot magma chamber could be rooted beneath the Datong volcanic field, and the center of the chamber is located beneath the northeastern Ordos block as shown in cross section AA′. Furthermore, the hot chamber beneath the Datong volcanic zone is found to be connected to a lesspronounced low-velocity zone in the upmost mantle north to the Taiyuan basin, as shown in cross section DD′. Contrasting to the strong low velocities beneath Datong volcanic zone, discrete high-velocity bodies are observed in the uppermost mantle in the south to central FWR as shown in cross sections CC′ and DD′. The strong lateral heterogeneities along the FWR play important roles in the development of FWR, which will be discussed in the latter section.

Ongoing Asthenosphere Upwelling Beneath the Datong Volcanic Zone
Compared to other intracontinental rifts in the world, volcanic activities in the Cenozoic FWR are rare and mainly restricted to its northern part, and most magmatic events in FWR occurred in the Datong volcanic zone during the Quaternary (Xu et al., 2005;Xu & Ma, 1992). Our results exhibit persistent low velocities in both the lower crust and uppermost mantle beneath the Quaternary Datong volcanic field ( Figure 14). As shown in cross section AA′ in Figure 16, our model resolves a significant LVZ (named L1) beneath the Datong volcanic zone, and it is deeply rooted up to a depth of 60 km. Moreover, it is located to the west of the Datong basin, which is consistent with the results from a recent magnetotelluric study (Zhang et al., 2016). Zhang et al. (2016) interpreted the highconductivity zone as magma source for the Datong volcanoes, and proposed that it was driven by the eastward escaping mantle flow from Tibetan. However, we image an another great LVZ (named L2) with its size larger than L1, beneath the northern TNCO at greater depths in the mantle, and this deeper LVZ is connected to the magnetotelluric conductive zone where L1 is located, as presented in cross sections DD′ and EE′. We propose that this deeper LVZ may play a significant role in the development of the Datong volcanoes and the magmatic source of the Datong volcanoes may be related to this deeper LVZ. The Datong volcanic zones may have a deeper and larger magma source besides the L1. Limited to the model resolution, we cannot image the root of L2, but it should be at least deeper than 120 km. Geochemistry analysis of volcanic rocks revealed that the Datong alkali lava suites have lower Al 2 O 3 and CaO, higher SiO 2 , and heavy rare Earth element contents and Na/Ti ratios, which are quite different from the Hannuoba basalts, a typical example of Miocene basalts in eastern China (Xu et al., 2005). The chemical composition of volcanic rocks reflects lower pressure of melt segregation and more depleted nature of the source of Datong alkali basalts, and the forward modeling for Datong basalts suggested that the Datong

10.1029/2018JB016476
Journal of Geophysical Research: Solid Earth volcanic rocks are generated at depth of~70 km, much shallower than that of Hannuoba basalts (Xu et al., 2005). Combined with the seismic evidence that strong low velocities (L1 and L2) appear beneath the Datong volcanic zone in upper mantle, we propose that the temporal thinning of the lithosphere beneath Datong region is due to a possible upwelling of hot asthenospheric materials. The possible mode for the lithospheric thinning and the generating of volcanic magma could be like this: the upwelling asthenosphere heats the uppermost mantle and produces L2 beneath the Datong region, then the hot materials and heat migrate upward to the shallow depth of 60~70 km, about the depth of L1. At the end, the lava separates in L1 and migrates to the crust and generates the relative homogeneous geochemical alkali basalts in the Datong volcanic zone (Xu et al., 2005). Considering the fast rifting rate at the junction between the FWR and Hetao rift (Shen et al., 2000), this hot material upwelling is probably an ongoing process.
Previous P and S wave tomography has imaged a significant low-velocity anomaly beneath the Datong volcanic zone which is thought to originate from at least the transition zone and extend to the base of the TNCO lithosphere (Zhao et al., 2009). Although the mechanism of the likely hot upwelling is still not very clear, its extending to more than 500-km depth indicates that it may be a mantle plume or a deeper convection process due to the possible delamination of the stagnant Pacific slab. It is believed that Cenozoic magmatism in the northern TNCO and modification of the lithosphere are related to the hot upwelling in the TNCO mantle (Jiang et al., 2013). Our model is concordant with previous studies, but resolves more detailed features on the crustal and upper mantle structure beneath northern TNCO.

Limited Mafic Underplating Beneath the TNCO Before the Rifting of the FWR
As presented in Figure 15a, clear high velocities in the lowermost crust are observed beneath the Weihe and Taiyuan basins, and the associated uncertainties presented in Figure 15d imply the robustness of the results. These high velocities have also been observed by receiver function studies with seismometers deployed in the Weihe basin (Xu et al., 2014) and the shear wave tomography in the south FWR (Ding et al., 2000). Considering the Moho uplifting beneath the FWR, especially beneath the Weihe and Taiyuan basins (Xu & Ma, 1992), it is not likely that the crustal eclogitization occurred there, as rocks shallower than~40 km could not be completely converted to eclogite (Ryan, 2001). Actually, geochemical study on mantle peridotite xenoliths indicates that the lithospheric mantle beneath the NCC has gone through repeated mantle metasomatic overprints (Tang et al., 2008;Zhang et al., 2009), and episodic widespread magmatic underplating beneath the NCC is inferred based on a comprehensive synthesis . A more recent shear wave velocity study also proposed that the mafic underplating occurred beneath the central NCC . Combining the facts that no significant low velocities are observed in the mantle beneath Weihe and Taiyuan basins, and the Cenozoic magmatism in the TNCO occurred mainly at early Tertiary (Jiang et al., 2013, and reference therein), we propose that the high-velocity layer observed here should be ascribed to the early mafic underplating before the rifting of the FWR. The underplating magma cooled down during the long period, thus eventually resulting in high velocities therein. As manifested by our lowermost crust model, the mafic underplating is mainly limited to regions which are of inherited weakness. It means that most magma flows related to the mafic underplating are trapped into the weakest parts of the TNCO, for example, the Weihe and Taiyuan basins.
On the other hand, tiny velocity jump can be found in the northeastern Ordos block where no mafic lower crust is imaged as shown in Figure 15b. As mentioned above, the low-velocity zone (L1) beneath the Journal of Geophysical Research: Solid Earth northeastern Ordos is connected to an ongoing asthenosphere upwelling in the mantle part. So there may be a thick crust-mantle transition zone beneath the northern TNCO. We propose that mantle material melts at the base of the crust may give rise to further crystallization, thus resulting in mafic lower crust beneath the northeastern Ordos block and the Datong volcanic field in the future.

Relics of Cratonic Root and Passive Mechanical Barriers in the Development of the FWR
The FWR exhibits almost consistently low velocities in the middle crust from north to south as presented in Figure 14c. In contrast, notable heterogeneities are imaged in the upper mantle. Specifically, three highvelocity anomalies (named H1, H2, and H3) are resolved with uncertainties less than 0.09 km/s beneath the central to southern part of the FWR (Figures 13, 16, and 17). H1 is located adjacent to the east of the Weihe basin, and H2 is located beneath the Lingshi block, which is south to the Taiyuan basin and north to the Linfen basin. In the deep upper mantle, H3 is located beneath the central FWR. The mafic underplating processes beneath the Weihe and the Taiyuan basins are bounded by the H1 and H2 anomalies, respectively. Meanwhile, the low velocities in the upper mantle beneath the Datong volcanic field are bounded by H3 to the south.
These high-velocity anomalies are interpreted as the ancient cratonic remnants that survived the Phanerozoic evolution of TNCO lithosphere. Since the Precambrian amalgamation of the eastern NCC and western NCC, the TNCO has experienced dramatic and complex lithospheric mantle modification. The related fundamental tectonic events include the Jurassic-Cretaceous amalgamation of the NCC to the Siberian plate in the north, and the Triassic collision between the NCC and South China Block in the south (Kusky et al., 2014;Wu & Zheng, 2013). The northern and southern segments of the TNCO are more likely to be affected during the lithospheric collisions, while the south-central FWR could retain its Archean cratonic root to a certain degree. Major and trace element analysis for xenolith-bearing Neogene basalts in Hebi County, which is located at the central TNCO, also indicates that the relics of the Archean lithosphere Journal of Geophysical Research: Solid Earth may preserve locally at a relatively shallow depth (Zheng et al., 2001). Therefore, the high-velocity anomalies in the south-central FWR are probably the relics of Archean cratonic root.
Geologic investigations revealed that the FWR evolved as a cluster of trans-tensional basins which were sheared dextrally and extended in the ESE direction (Zhang et al., 1998). On the other hand, GPS measurements observed that the northeastern terminus of the NNE-SSW striking FWR is experiencing horizontal extension along N105°E direction (Shen et al., 2000). The discrepancy between the geologic and geodetic surveys might be due to the differential opening rates along the FWR (Shen et al., 2000). It is reasonable to speculate that the opening rate along FWR increases from south to north, as the high-velocity bodies of H1 and H2 may act as passive mechanical barriers that halt the northeastward advancing of the rift system. In addition, a synthetic view based on geologic data and remote sensing images suggests that the FWR has experienced multistage rifting processes (Zhang et al., 1998). The earliest extension occurred at the Weihe basin in the late Eocene to early Oligocene. When the extension process propagated to the central, and afterward to the northern parts of the FWR during the Oligocene-Pliocene, the rifting rate in the Weihe basin has been accelerated twice (Bao et al., 2013;Zhang et al., 1998). Therefore, we speculate that the H1 and H2 acted as rigid blocks and impeded the extension process on a lithosphere scale throughout the rifting evolution of the FWR.

Different Rifting Mechanisms Between the Southern and Northern Parts of the FWR
It is well known that the Cenozoic FWR developed from south to north (e.g., Bao et al., 2013;Jiang et al., 2013). The appearance of the large-scale velocity anomaly beneath the northern FWR may imply different rifting mechanisms between the southern and northern parts of the FWR. In addition, the difference of seismogenic depths between the southern and the northern FWR also implies the different rifting mechanisms in the two parts, and the decreasing trend of the seismogenic depths proposes that the extension of FWR is still an ongoing process (Song et al., 2012).  The common basaltic underplating in NCC is thought to take place before the Cenozoic (Fan et al., 1998;Liu et al., 2001). Specifically, the underplating beneath the Weihe and Taiyuan basins occurred much earlier than the inception of FWR rifting process as discussed before. Actually, Cenozoic magmatism activities are rare in the southern FWR in which the rifting process commenced (Xu & Ma, 1992). Therefore, it is not likely that the early extension of the FWR was induced by the upwelling of mantle materials. Considering the concordant ages of the rifting in the Weihei basin and the uplifting of the Tibetan Plateau, we speculate that the early uplift of the Tibetan Plateau enhanced regional extension stress field and caused the rifting process in the adjacent Weihe basin, since it is weaker compared to adjacent Archean blocks.
Because of the counterclockwise rotation of the Ordos block induced by the continuing extrusion of the Tibetan Plateau (Zhang et al., 1998), the rifting process in the southern TNCO propagated to the northern TNCO. However, the relics of cratonic root beneath the southern to central TNCO (H1 and H2) acted as mechanical barriers and blocked this process, thus caused the multistage passive extension in FWR, which has also been evidenced by geological data and SPOT imagery analysis (Zhang et al., 1998).
In contrast to the southern TNCO, the northern segment has experienced repeated volcanic eruptions during the Neogene and Quaternary. As discussed before, geochemical studies suggested that the Cenozoic basalts in the northern TNCO are derived from an asthenosphere source, and the northern TNCO has gone through progressive lithospheric thinning during the Cenozoic (e.g., Xu et al., 2004Xu et al., , 2005, which is evidenced by the strong low velocities beneath the northern TNCO in the uppermost mantle in this study (Figures 14,16,and 17). Combined with previous body wave tomography (Zhao et al., 2009), we propose that this LVZ is a robust feature which extends to at least the mantle transition zone. Although the origin of LVZ is still enigmatic, it could be either a mantle plume or an upwelling asthenosphere initiated by the subduction of the Pacific plate (Huang & Zhao, 2006). Naturally, the mantle upwelling plays a key role in the rifting in the north FWR and the thinned lithosphere facilitates the opening in the north TNCO. Overall, the rifting of the north FWR is due to the combined effects of the rotation of the Ordos block and the mantle upwelling, in contrast with the passive extension in the south, which can be illustrated in a 3-D visualization structure ( Figure 17).

Conclusions
A high-resolution 3-D crustal and uppermost mantle S wave velocity model of Fenhe-Weihe Rift system is built by jointly inverting the harmonic stripping receiver functions and Rayleigh wave phase and group velocity dispersions derived from ambient noise cross correlations. Along the FWR, significant heterogeneities are observed in both the crust and upper mantles. The lower crust beneath the Weihe and the Taiyuan basins are imaged as high-velocity layers. We interpret them as basaltic underplating that occurred before the inception of the FWR. The rare magmatic activities in the southern FWR indicate that the development of the southern FWR may be not triggered by the mantle upwelling but by the uplift of the Tibetan Plateau. In the southern-to-central FWR, discrete high-velocity anomalies are seen in the upper mantle, which may be the relics of ancient cratonic root and play as passive mechanical barriers in the development of FWR, consisting with the multistage extension processes in the southern FWR. Contrarily, the large-scale, low-velocity zone beneath the north TNCO is interpreted as ongoing asthenosphere upwelling, which may facilitate the opening of the northern FWR. In this scenario, the rifting mechanisms in the southern and northern parts of the FWR are different. The southern part has undergone a passive extension process due to the extensional forces caused by the uplift of the Tibetan Plateau, while the rifting in the northern part is due to the joint effects of the asthenosphere upwelling and the rotation of the Ordos block.