Geodetic constraints on three-component motion of the Ordos block (China) and their implications for lithospheric dynamics Open Access

Distributed deformation across broad networks of active fault systems accommodates crustal deformation driven by the India-Eurasia convergence (Holt et al., 2000; Thatcher, 2009; Zhang et al., 2004). While some regions appear to not experience active deformation, other regions experience localized strain and related higher seismic risks (Molnar and Dayem, 2010). The Ordos block, lying northeast of the northeastern Tibetan Plateau, is an example of such rigid block tectonics that contrasts with its periphery, which has experienced multiple stages of deformation controlled by outward growth of the northeastern Tibetan Plateau and subduction-driven rollback of the western Pacific slab (Fig. 1A; Molnar and Tapponnier, 1977; Northrup et al., 1995; Shi et al., 2020; Su et al., 2021). Recently, it has been suggested that the Ordos lithosphere is encountering regional modification by mantle convection caused by geodynamic processes, which contribute to observed surface deformation patterns by coupling with the brittle crust, resulting in dynamic topography and/or isostatic adjustment (Flament et al., 2013; Guo and Chen, 2017; Yu and Chen, 2016; Wu et al., 2022b; Yu et al., 2021; Guo et al., 2022; Kreemer et al., 2020). The Ordos block thus provides a rare natural laboratory to image surface deformation and examine the correlations between geodynamic processes, lithospheric thinning, and mantle flow.

Sitting within the India-Eurasia zone of N-S convergence, uplift of the Ordos block initiated in the Late Cretaceous and experienced two-stages of tectonic evolution during the Cenozoic (Ren et al., 2002; Shi et al., 2020). The first stage was accompanied by pure shear in its northwestern and northeastern margins, with NW-SE extension and basin inversion driven by rollback of the subducting Izanagi and Pacific plates and the spreading ridge in between (Fan et al., 2019; Liu et al., 2019; Northrup et al., 1995; Shi et al., 2015; Schellart and Lister., 2005; Su et al., 2021; Wu and Wu., 2019). During the second stage, the Shanxi and Hetao grabens (Fig. 1) developed, driving intensive mountain-building over three alternating episodes of shortening and extension resulting from northeastward growth of the Tibetan Plateau and Pacific subduction (Li et al., 2012; Shi et al., 2020; Ye et al., 1987; Zhang et al., 2006).

Contemporary horizontal deformation surrounding the Ordos block is mainly characterized by NE-SW–directed contraction and NW-SE–directed extension with attendant minor block rotations and fault slip (Chen et al., 2005; Li et al., 2021, 2023a), while strain rates inside the Ordos block are negligible (Cui et al., 2016; Li et al., 2020b; Middleton et al., 2017, 2018; Wang and Shen, 2020; Zheng et al., 2017). Geodetic and geologic observations have been used to propose a “book-shelf” model to explain the Cenozoic deformation within the North China Craton (NCC) via a large-scale WNW-ESE left-lateral shear zone induced from differential motion between the Amurian and the South China blocks (Fig. 1A; Xu and Ma., 1992; Zhang et al., 2018). Various kinematic models disagree on the net horizontal movement of the Ordos block compared to its surrounding regions. Results include anticlockwise rotation (Xu et al., 1993; Xu et al., 1994; Peltzer and Saucier, 1996; Chen et al., 2005; Cheng et al., 2021, 2015; Middleton et al., 2016; Zhang et al., 1995), simple shear (Shi et al., 2020), and “frame-wobbling” (Hao et al., 2021).

Contemporary vertical land motion (VLM) of the Ordos block has been unclear, hindering understanding of its kinematics and the potential importance of geodynamic processes, such as lithospheric thinning and mantle flow. Hao et al. (2016) presented a vertical velocity field of the Ordos block and its vicinities calculated from high-precision leveling observations, though the leveling routes are sparse in the Ordos block and the Qinling Orogenic Belt. In addition, only six campaign-mode Global Navigation Satellite System (GNSS) vertical velocities were used to constrain the velocity reference frame of the leveling observations, which are too sparse to remove systematic errors that accumulate with distance along the leveling routes over a large region. Moreover, they did not consider the effects of non-tectonic deformation from surface mass changes, such as may occur from changes in the distribution of hydrological loads. GNSS vertical velocities show that the Ordos block uplifts at a rate of ~1–2 mm/yr, while its surrounding regions, including the eastern Qinling Orogenic Belt and grabens, experience near zero vertical motion or subsidence (Fig. 1; Pan et al., 2021; Zhao et al., 2023). However, the geographic density of GNSS stations is too sparse to reveal the details of VLM within the Ordos block and its surrounding regions.

We here focus on the 3-D crustal movement of the modern Ordos block by taking advantage of geographically dense and extensive coverage of combined continuously operating GNSS and historical leveling networks. To focus on long-term motion, we estimate and subtract the contribution to VLM from short-term surface hydrological mass variations obtained from Gravity Recovery and Climate Experiment (GRACE) observations and from the effects of glacial isostatic adjustment (GIA) from the ICE-6G_D (VM5a) model of Peltier et al. (2018). The VLM predictions from these models are used to correct the geodetic data for transient processes that are not VLM attributable to long-term tectonic or other geodynamic processes. Following these corrections, we obtain an interpolated VLM map from corrected geodetic data by applying a robust network imaging (RNI) algorithm of Kreemer et al. (2020). Finally, we discuss the constraints that VLM and horizontal strain rates place on dynamic mechanisms at work in the Ordos block.

We collected data from 90 continuous GNSS stations operating between 2010 and 2022 from the Crustal Movement Observation Network of China (CMONOC I, II) and Shaanxi Earthquake Agency, China Earthquake Administration (Fig. 1A). These data are processed using the GAMIT/GLOBK10.6 software with a rigorous processing strategy (Herring et al., 2015; Hao et al., 2021) to obtain the daily free station position time series relative to the International Terrestrial Reference Frame (ITRF2014; Altamimi et al., 2016) (Fig. S1¹). We then estimate the vertical velocity at each GNSS station (GNSS VLM) by applying the robust MIDAS trend estimator (Fig. 1A; Fig. S3A; Blewitt et al., 2016), which is designed to be robust against outliers, seasonality, skewness, undocumented steps, and heteroscedasticity. The MIDAS algorithm estimates the trend of the position time series by calculating the median values of geodetic velocities from all data pairs with a time difference of one year and can produce accurate velocity with realistic uncertainties in blind tests. In order to provide a comparative analysis, we have employed the Hector software to estimate GNSS vertical velocities (Fig. S3A). This software allows for the fitting of GNSS vertical position time series with a model that incorporates offsets resulting from instrumental changes and earthquakes, as well as a trend, and a constant semi-annual and annual seasonal signal (Bos et al., 2013). The median difference of 0.07 mm/yr between these two GNSS vertical velocity fields suggests that both methods can give internally consistent estimations of the vertical velocities. The fitted amplitude and corresponding phase of annual periodic vertical displacements of each GNSS station are presented in Figure S3B. The MIDAS GNSS vertical velocities have uncertainties of 0.4–1.0 mm/yr (Fig. 1A).

Adding to the data published by Hao et al. (2016), we collected leveling measurements spanning 1000 km across the eastern Qinling Orogenic Belt. These measurements were taken twice, once within the interval of 1982–1985 and again in the interval 2005–2006, using the second-order leveling standard. In 2017, we conducted a third measurement using the first-order leveling technique. We also collected leveling data from the western Ordos block, covering a distance of 500 km. These measurements were taken twice, once from 1985 to 1986 and again in 2013, using the second- and first-order leveling standards separately. The uncertainties for the first- and second-order leveling are less than 0.45 mm/km and 1.0 mm/km, respectively. The misclosures from the forward and backward leveling of all routes are within ±0.5 mm/km. To reduce the systematic errors accumulating along the leveling routes, we operated forward and backward surveys along each route with the misclosure to be within ±0.4 mm/km (Lai et al., 2004; Hao et al., 2016). Moreover, we detect and remove unstable benchmarks by using the height differences time series between adjacent benchmark pairs. We show the leveling routes and surveying epochs in Figures 1B and 2A. We also show the root mean square (RMS) residuals of observation time in each epoch with the integer year removed to show the repeatability in the month of our observations, which can help reveal influences from seasonal non-tectonic deformation (Fig. 2B). The plot shows that more than 70% and 90% of stations have observational month RMS residuals smaller than 2 and 3 months, respectively. The months that leveling observations were made are mainly from April to August, which are consistent with the peak time of seasonal non-tectonic deformation induced from the GNSS position time series (Fig. S3B). This suggests that the seasonal non-tectonic signals can be removed effectively based on the multiple observations.

In the least squares adjustment, the inferred posterior uncertainty estimation of 1.26 mm/km is a little larger than the total mean square error of 1.0 mm/km in magnitude, which is estimated from misclosures and related to the length of the leveling circuit. This indicates that the assumption of the vertical crustal movement following a linear trend is reasonable. The slight discrepancy in between may result from the significant subsidence due to groundwater extraction within the grabens surrounding the Ordos block (Liu et al., 2023; Zhao et al., 2019, 2023). The inferred leveling vertical velocities are in the range of −6 mm/yr to 6 mm/yr with respect to the ITRF2014 origin (Fig. 1B), and are much greater than their uncertainties (0.4–1.6 mm/yr).

We use the VLM rates at 35 continuously operating GNSS stations, which are each within a few kilometers from the nearest leveling benchmark, as prior constraints to define the reference frame. Twenty-five vertical velocity differences between these nearly collocated GNSS sites and leveling benchmarks are smaller than 1 mm/yr. The inferred median absolute deviation (MAD), expressed as MAD = median, where the former is GNSS observed and the latter is the imaged leveling vertical velocities, which is equal to 0.85 mm/yr (Fig. 3). This misfit is similar to the 1σ uncertainty of GNSS and leveling VLM, which suggests that there is close alignment between the GNSS and leveling vertical velocities.

To validate the alignment between GNSS and leveling VLM across the entire study area, we interpolate the leveling vertical velocities onto all 90 continuously operating GNSS stations by applying the RNI approach (Fig. S3A, the rightmost arrows). The inferred median difference and MAD between observed and imaged leveling vertical velocities on GNSS sites are 0.02 mm/yr and 1.3 mm/yr, respectively.

It has been noticed that geodetic data, e.g., GNSS and leveling observations, contain non-tectonic solid Earth deformation signals resulting from mass variations due to the redistribution of hydrological and non-tidal atmospheric mass loads (Blewitt et al., 2001; Pan et al., 2021; Zhang et al., 2021). To identify these non-tectonic deformation signals, we estimate the trend of VLM from the GRACE data as if it were completely explained by terrestrial water storage, atmospheric, and other non-tidal mass variations (Fig. 4A). We use the spherical harmonic coefficients of GRACE products and the AOD1B de-aliasing model (GAC) released by GFZ (Dobslaw et al., 2018). In our data processing, we remove the correlated errors in the original GRACE data applying a decorrelated filter algorithm proposed by Swenson and Wahr (2006). We also use a 350 km radius anisotropic Gaussian filter to compress gravity field high-frequency noise (Wahr et al., 1998). Additionally, to overcome the insensitivity of the GRACE satellites to the degree 1 and 2 components of the gravity field, we replace the degree 1 coefficients following the strategies of Swenson et al. (2008) and the C20 coefficients from the satellite laser ranging observations (Cheng et al., 2013). The VLM resulting from GIA is demonstrated by the ICE-6G_D (VM5a) model (Fig. 4B; Peltier et al., 2018). Noting the different time spans of GRACE and leveling data, we argue that the total magnitudes of these rates are much smaller than rates of the Ordos uplift and have relatively little impact on our analysis. These non-tectonic deformation signals are then subtracted from the geodetic VLM (Fig. 4C).

Vertical velocities resulting from surface mass variations including hydrology and atmosphere are in the range of −0.4 mm/yr to +0.4 mm/yr in our study area (Fig. 4A). The crust in southwest South China to northwest in the Alxa block (west of the Ordos block) subsided with gradually increasing subsidence rates southward. This suggests that surface mass gain resulting from precipitation was increasing there, but less so to the north. However, in the Ordos block and North China Plain, the crust uplifts with rates smaller than 0.4 mm/yr, indicating negative surface mass changes from decreasing precipitation and serious groundwater depletion (Zhang et al., 2021; Pan et al., 2021; Zhao et al., 2019). The amplitude of non-tectonic crustal uplift from water variations is likely underestimated in the North China Plain and grabens surrounding the Ordos block due to GRACE’s low spatial resolution of 300–400 km. So, we here mainly focus on the large-scale uplift of the Ordos block and surrounding regions.

Predicted VLM from GIA is 0.2–0.3 mm/yr in our study area with the rate decreasing from northwest to southeast (Fig. 4B). The variation of predicted VLM from GIA is less than 0.06 mm/yr in our study area and is less than the uncertainties in the VLM data. While the signals from both hydrological loading and GIA are small (Figs. 4A and 4B), we subtract them from our geodetic vertical velocities (Fig. 4C) to remove these non-tectonic signals from the dataset.

Considering that vertical crustal deformation in the Ordos block and its vicinities contains both potential long-wavelength (e.g., tectonic uplift) and short-wavelength (e.g., groundwater extraction, volcanic activity) signals, the corrected geodetic vertical rates (Fig. 4C) are imaged onto a 0.2° × 0.2° regular grid (Fig. 4D) using the RNI algorithm (Kreemer et al., 2020) based on the Delaunay triangulations. RNI is a variant of the GPS imaging algorithm of Hammond et al. (2016). Both algorithms have been proven to be efficient in capturing different wavelength signals based on stations as local as possible and, therefore, can be used to image all kinds of geospatial observations. The RNI methodology contains two steps. First, the vertical velocity at each geodetic site (v_i) is substituted with a median/despeckled velocity deduced from velocities at M local geodetic sites, which contain the site itself plus sites connected to it via Delaunay triangulation of the station locations, plus those sites that are closer to it than the median distance between it and the sites connecting to it. Second, we use the despeckled velocities and their uncertainties to calculate the weighted vertical velocity at a gridded cluster of N points, j, which is within the convex hull of all site locations. For each evaluation point on the grid, we construct a new Delaunay triangulation based on all of the GNSS site locations and the current evaluation point, and the weighted median vertical velocity at each evaluation point (v_j) is then calculated from the despeckled rates defined as the first step.

To investigate the horizontal crustal motion of the Ordos block and the potential geodynamic mechanisms suggested by our VLM map, we model a robust horizontal strain rate field onto a 0.2° × 0.2° regular grid using the median estimation of local deformation (MELD) algorithm of Kreemer et al. (2020) based on the dense GNSS velocities published by Hao et al. (2021) (Figs. 5A). MELD uses a geographically variable spatial resolution that depends on the local density of contributing GNSS stations. The spatial resolution of our preferred strain rate model is, therefore, variable, but in most areas has a length scale of <60 km (Fig. S5B). To emphasize the point that the strain rates inside the Ordos block are not significant, in contrast to significant rates around the block perimeter (Fig. 5A), we plot the signal-to-noise ratios (SNR) of the second invariant strain rates in Figure 5C. As presented by Hao et al. (2021), horizontal deformation inside the Ordos block is negligible and much smaller than those in its surrounding regions especially the northeastern Tibetan Plateau (background in Fig. 5A). The deformation style is shown with background color in Figure 5B. We define it as , where and are the largest and smallest eigenvalue of horizontal strain rate tensor. The definition means that the values of 1.0 and −1.0 represent the pure extension and contraction, and 0.0 indicates strike-slip sense. The positive (negative) values imply extension (contraction). Additionally, we use the normalized principal strain rates, defined as to reveal extension (dark green arrows) and contraction (blue arrows) directions (Fig. 5A). We also calculate horizontal velocities with respect to the Ordos block reference frame at grid points based on the modeled robust rotation (Fig. S5A) and strain rates (Fig. 5A; green arrows in Fig. 5B).

Both GNSS and leveling VLM show significant crustal uplift in and surrounding the Ordos block, while regional subsidence is widespread (Figs. 4C, 4D, and 6; Liang et al., 2013; Pan et al., 2018, 2021; Zhao et al., 2023). These observations reveal the complexity of the VLM. Our geodetic VLM shows 0.5–4.0 mm/yr crustal uplift of the Ordos block with the maximum occurring within the northeastern corner (site 4 in Fig. 4C), and significant subsidence in the surrounding grabens, such as the Weihe and Hetao grabens (Fig. 4C). The grabens are areas that have serious groundwater extraction for industry, irrigation, and municipal uses (Figs. 4C, 4D, and 6A–6D). These VLM patterns are similar to those from GNSS vertical velocities alone (Pan et al., 2018; Pan et al., 2021; Zhao et al., 2023), except that the uplift of the northeastern Ordos block and Datong areas is revealed by our inclusion of leveling data (sites 4 and 8 in Figs. 4C and 4D). In the area between the Ordos block and the northeastern Tibetan Plateau (site 7 in Figs. 4C and 4D), the uplift rate is 1–3 mm/yr, while small-scale subsidence occurs near the eastern Haiyuan fault (site 9 in Figs. 4C and 4D; Su et al., 2018; Pan et al., 2021). This subsidence may correlate with the extensional tectonics in the eastern end of the Haiyuan fault, as well as the post-seismic deformation of the 1920 Haiyuan earthquake (Ou et al., 2020). In the Qinling Orogenic Belt, the crust uplifts in the west but subsides in the east (sites 5 and 6 in Figs. 4C and 4D), which is also confirmed by GNSS and interferometric synthetic aperture radar (InSAR) observations from previous studies (Liu et al., 2023; Zhao et al., 2023), although the rates are insignificant compared to the uncertainties of 0.4–1.6 mm/yr we find in our analysis. While the density of our geodetic observations in the Alxa block and Yinshan-Yanshan Mountains is sparse, the available data show crust uplift at a rate of 1–3 mm/yr in these regions.

The second invariant strain rates (Fig. 5A) show that the Ordos block is rigid and tectonically stable compared to its surrounding regions, which is consistent with observations from seismicity and active tectonics (Xu et al., 1994; Chen et al., 2005; Zhao et al., 2017; Middleton et al., 2016, 2017, 2018). The contraction and extension directions surrounding the Ordos block are mainly NE-SW- and NW-SE-directed, respectively (Figs. 5A and 5D). They appear a little disordered within the Ordos block because of low SNR except in its northeastern corner (Figs. 5A and 5C; Cui et al., 2016; Hao et al., 2021; Li et al., 2020b; Middleton et al., 2017, 2018; Su et al., 2018; Wang and Shen, 2020). The contraction direction west of the Ordos block is mainly northeastward, kinematically consistent with shortening associated with the northeastward growth of the Tibetan Plateau related to extrusion from collision with the Indian Plate. In the Shanxi rift, lying east of the Ordos block, the orientations of contraction are much more complex, varying from NE-SW to W-E. These directions may be related to deformation of the lithosphere associated with thermally triggered subduction and subsequent rollback of the western Pacific slab as well as northeastward growth of the Tibetan Plateau (Northrup et al., 1995; Schellart and Lister., 2005; Ye et al., 1987; Zhang et al., 2006; Li and Kreemer., 2021; Shi et al., 2020; Su et al., 2021). The contraction directions along the Qinling Orogenic Belt vary from NE-SW in the west to NW-SE in the central and nearly W-E in the east (Fig. 5A; Cui et al., 2016; Middleton et al., 2018), which suggest complex geodynamic mechanisms.

Low rates of extension internal to the Ordos block mainly occur to the northeast and southwest of the block and along the northern Shanxi rift and Taihang Mountains (sites 1–3 in Fig. 5B; Northrup et al., 1995; Schellart and Lister, 2005; Hao et al., 2021; Li et al., 2020b; Middleton et al., 2017; Zhao et al., 2017). In the northeastern Tibetan Plateau, Qinling Orogenic Belt, and Yinshan-Yanshan Mountains, crustal horizontal deformation is mainly characterized by contraction associated with convergence (sites 4 and 5 in Fig. 5B; Su et al., 2018; Wang and Shen, 2020). Strike-slip mainly occurs along the western and eastern peripheries of the Ordos block (Figs. 5B and 5D; Zhao et al., 2017; Wang and Shen., 2020; Zheng et al., 2017; Hao et al., 2021), and a weak sinistral shear exists along the Hetao graben. There is no detectable strike-slip movement along the south of the Ordos block (Hao et al., 2021). Hao et al. (2021) reported a “frame wobbling” effect, which suggested that the rigid Ordos block settles steadily, whereas the surrounding blocks move in different directions and at different rates, therefore exerting different forces across the bounding deformation zones, such as the grabens and fold-thrust faults zones (Fig. 5D). For example, with respect to the stable Ordos block, the North China Plain to the east moves southward at a rate of 1.2 mm/yr causing dextral shear across the Shanxi rift (Hao et al., 2021), while the central Shanxi rift does not demonstrate extensional motion presently. The extension on the northern Shanxi rift is consistent with deformation of the tail end of the strike-slip fault along the northern boundary of the North China Plain.

Our VLM map shows that the Ordos crust experiences widespread uplifts, though the vertical velocities have large lateral variability in its regional patterns especially along a transect from northeast to southwest (Figs. 4C, 4D, and 6F). Within the northern Ordos block (site 1 in Figs. 4C and 4D), the crust uplifts at a rate of 0.5–1.0 mm/yr, less than the 1.5–2.0 mm/yr within the central and southern Ordos block (sites 2 and 3 in Figs. 4C and 4D). An exception to this pattern is the northeastern Ordos block, which experiences a more rapid uplift of 2.0–4.0 mm/yr (site 4 in Figs. 4C and 4D, Figs. 6E and 6F).

To investigate the observed NE-SW uplift differences, we introduce seismic tomography data from Wu et al. (2022a), azimuthal anisotropy data from Wu et al. (2022b), and shear-wave splitting observations (SKS) from Liu et al. (2021) (Fig. 7). Distinct NE-SW differences are also observed by seismic S-wave velocity and azimuthal anisotropy contrasts at different depths, SKS splitting orientations and lithospheric thickness variations (Figs. 7 and 8). Significantly lower seismic S-wave velocity and thinner lithosphere are imaged underneath the northeastern Ordos block, which is also characterized by more rapid crustal uplift (site 4 in Fig. 4C, Figs. 6F and 8), compared to the central and southwestern Ordos block (sites 2 and 3 in Figs. 4C and 4D).

Seismic velocity contrast, electrical heterogeneity, and seismic anisotropy variations have revealed complex structural properties beneath the Ordos block (Fig. 7; Dong et al., 2014; Gao et al., 2018; Liu et al., 2021; Lv et al., 2021; Wu et al., 2022a, 2022b; Yu et al., 2021). It has been suggested that the heterogeneous thinning of the lithosphere underneath the northern and northeastern Ordos block and its vicinities results from thermomechanical erosion by vertical mantle upwelling (Li et al., 2018, 2020a; Lv et al., 2021; Sun and Liu., 2023; Wu et al., 2022b; Zhang et al., 2022). The resulting significant contrast of lithospheric thickness is consistent with small-scale mantle convection especially if the viscosity of the asthenosphere is reduced by melts and volatiles from dehydration, which can further modify the lithosphere (Fig. 8; Conrad et al., 2011; Jiang et al., 2013; Li et al., 2020a; Sun and Liu., 2023). Other previous studies have reported horizontal asthenospheric mantle flow around the Ordos block, which may result in the ongoing reworking of the Ordos lithosphere resulting from thermal-chemical erosion caused by horizontal flow (Fig. 8; Guo et al., 2022; Yu et al., 2021; Zhang et al., 2022; Lv et al., 2021). Given that the negative Bouguer gravity anomalies in the Ordos block suggest isostatic adjustment in the crust (Wang et al., 2014), we, therefore, suggest that the observed differential VLM between the northeastern and southwestern Ordos block may reflect vertical re-equilibration resulting from the heterogeneous thinning and reworking of the lithosphere caused by mantle convection.

The India-Eurasia collision and post-collisional convergence may also contribute to the overall NE-SW contraction and differential uplift in this region (Figs. 5D and 8; Guo et al., 2015; Tian et al., 2021; Su et al., 2021). However, it has been suggested that VLM induced from horizontal shortening, with the assumption of an incompressible crust and upper mantle, is smaller than 0.2 mm/yr in the Ordos block (Zhao et al., 2023). This is much less than the observed geodetic VLM, which exceeds 1.0 mm/yr within most of the Ordos block. Meanwhile, the low-temperature thermochronological data and sedimentary records from both northern and southern parts of the Ordos block reveal that Eocene rifting of the Shanxi rift system is attributed to the far-field effect of subduction of the Izanagi and western Pacific plates and the spreading ridge in between (Clinkscales et al., 2021, 2020; Fan et al., 2019; Su et al., 2021; Wu and Wu., 2019), whereas the Miocene rifting results from the eastward crustal extension of the Tibetan Plateau (Cheng et al., 2021; Ren et al., 2002; Su et al., 2021; Xu et al., 1993; Ye et al., 1987). The synchronous rifting in the northern and southern parts of the Shanxi rift suggests that the northern and southern parts of the Ordos block were likely under the same tectonic setting, sharing the same geodynamic crustal process. We thus propose that the deep processes might have played a dominant role in the observed large-scale Ordos block uplift contrast, whereas the short-wavelength variations are related to regional tectonism in the crust. We recognize that the geodetic VLM is that of a singular (present-day) time shot, which may not fully represent that of long-term crustal motions within the geological time scale. However, the ongoing lithospheric thinning and reworking around the perimeter of the Ordos block and the presence of negative Bouguer gravity anomalies suggest that mantle convection beneath the Ordos block may be responsible for its observed differential uplift. Additionally, we cannot exclude the influence from flexural loading from northeastern Tibet, which requires further study.

Crustal uplift rates are 2–4 mm/yr in the northeastern Ordos block and Datong volcanic area (DTV) (sites 4 and 8 in Figs. 4C and 4D), and are 2–3 mm/yr larger than those in the other parts of the Ordos block (Figs. 6D and 9A). We also find significant horizontal extension, with dilatational rates larger than 1.5 × 10⁻⁹/yr (sites 1 and 3 in Fig. 5B; Figs. 6A, 6E, 6F, and 9B), accompanied by a pattern of crustal contraction in its surrounding regions (Fig. 9B). The high SNR of second invariant strain rates, larger than three (dark green box in Fig. 5C), suggests that these deformation patterns are well resolved in the dense GNSS data (Hao et al., 2021). We thus consider this 3-D geodetic deformation anomaly as the Datong volcanic anomaly.

This area also exhibits widespread seismic low-velocity anomalies in the middle-to-lower crust and upper mantle, revealing a deep mantle upwelling associated with the western Pacific slab subduction (Figs. 7A–7C; Lei., 2012; Guo et al., 2016; Li et al., 2018; Cai et al., 2021; Wu et al., 2022b). We, therefore, attribute this geodetic anomaly to mantle upwelling and extension in the DTV and surrounding regions (Liu et al., 2004; Jiang et al., 2013; Lei, 2012; Li et al., 2018; Wang et al., 1989; Wu et al., 2022a, 2022b). Similar 3-D deformation anomalies have been identified in geodetic data, such as that in the Eifel volcanic area in Germany (Kreemer et al., 2020). A volcanic related 3-D crustal deformation anomaly also occurs in the Xiyang-Pinding volcanic area (XPV) (Figs. 5B and 9B), which has slower uplift compared to that in the DTV but is also above the mantle with low seismic velocity anomalies (Figs. 7A–7C; Cai et al., 2021; Guo et al., 2016; Wu et al., 2022b). The 3-D deformation field in the Hebi volcanic area shows differences compared to the DTV and XPV, but it is unreliable because of the lack of data resulting in large uncertainties in the VLM map (Figs. 4D and 5B).

With the ongoing India-Eurasia post-collisional convergence, the mantle and crustal materials in the northeastern Tibetan Plateau are extruded outward (Tapponnier et al., 2001; Molnar and Tapponnier, 1977). The predominantly NW-SE–oriented fast directions of the seismic anisotropy in the lower crust and upper mantle are generally parallel to the direction of large-scale faults and tectonic boundaries, e.g., the Haiyuan fault (Figs. 7B–7E). These observations suggest that lithosphere of the northeastern Tibetan Plateau deforms coherently and is extruded southeastward. However, directions of fast polarization of seismic anisotropy systematically deviate from absolute plate motion (APM) directions, which is dominantly northeastward or eastward (Fig. 7A; Kreemer, 2009; Li and Kreemer, 2021). APM is the surface expression of a dynamic mantle and is the net-result of the various forces acting on the lithosphere (Li and Kreemer, 2021). This difference is attributed to the presence of sub-asthenospheric mantle flow, which leads to anisotropy caused by shearing due to differential motion between the lithosphere and such flow (Li and Kreemer, 2021; Silver and Holt, 2002). The inferred flow is eastward-directed at 1–2 cm/yr depending on an APM frame with an intermediate global net-rotation underneath East Asia (Kreemer, 2009; Li and Kreemer, 2021; Jolivet et al., 2018).

The eastward expansion of the northeastern Tibetan Plateau is blocked by the stable Ordos block, resulting in the significant crustal horizontal contraction and anticlockwise rotation, and uplift across its margins (Figs. 4D and 5A; Hao et al., 2021; Li et al., 2020b). Low rate deformation and the presence of seismic high-velocity within the western and southwestern Ordos crust indicate that in these regions it appears not to deform (Wang and Shen, 2020; Hao et al., 2021; Li et al., 2020b). In contrast, the Qinling Orogenic Belt is characterized by relatively low-seismic velocity anomalies at depths from 150 km to 300 km, and a strong anisotropic layer having large SKS-wave splitting delay times (Fig. 7E; Guo and Chen, 2017; Lv et al., 2021; Yu et al., 2021; Liu et al., 2021; Wu et al., 2022b), which together indicate the presence of weak and possibly flowing asthenosphere underneath this region.

The direction of flow of crustal materials along the Qinling Orogenic Belt, however, is still controversial. Some studies suggest that crustal materials in the northeastern Tibetan Plateau are extruded eastward along the Qinling Orogenic Belt and/or beneath the southern Ordos block (Molnar and Tapponnier, 1977; Clark and Royden, 2000; Royden et al., 1997), while other recent studies assert that the mechanically strong middle-to-lower crust underneath the central segment of the Qinling Orogenic Belt impedes the eastward movement of the channelized crustal flow from the northeastern Tibetan Plateau (Guo and Chen, 2017; Jiang et al., 2016; Song et al., 2018).

Our three-component geodetic observations of crustal deformation provide a means to distinguish between competing geodynamic hypotheses. The VLM maps show that the crust in the western Qinling Orogenic Belt and the northeastern Tibetan Plateau uplifts at a rate between 0.5 mm/yr and 2.0 mm/yr, particularly in the region west of ~106°E longitude (sites 5 and 7 in Figs. 4C and 4D; Figs. 6C and 6D). The surface uplift rate decreases and transitions to subsidence near the east boundary of the Hannan-Micang (HNMC) crustal-scale high seismic velocity body (known as “Hannan Dome”) (Figs. 4C and 4D; Hu et al., 2006; Jiang et al., 2016; Song et al., 2018). To the east of this boundary, the crust subsides at ~0.5 mm/yr (site 6 in Figs. 4C and 4D; Fig. 6D), although it falls within the range of data uncertainty (0.4–1.6 mm/yr). The weak horizontal deformation along the whole Qinling Orogenic Belt is contractional with varying directions (site 5 in Fig. 5B; Cui et al., 2016; Li et al., 2020b; Su et al., 2018) and sits amid the collision between the NCC, South China, and eastward extrusion of the northeastern Tibetan Plateau. The modeled GNSS horizontal velocity field relative to the Ordos block, inferred from the strain rate model, indicates no eastward components in the east of the HNMC with respect to the stable Ordos block, but rather slow eastward movements in its west (green arrows in Figs. 5B and 5D).

To further examine the lithospheric deformation and crustal flow beneath the Qinling Orogenic Belt corresponding to the expansion of the northeastern Tibetan Plateau, we compare horizontal geodetic surface motions relative to the APM reference frame of Kreemer (2009) with directions of fast polarization of seismic anisotropy. APM directions along the Qinling Orogenic Belt align with directions of fast polarization of seismic anisotropy and subparallel to the structural axis of the orogeny (Figs. 7A–7E), which reveals coherent deformation. This deformation pattern reflects the lithospheric fabric created between the N-S collision of the NCC and Yangtze Craton during the Triassic (Guo and Chen., 2017; Yu and Chen, 2016; Wu et al., 2022b). Moreover, the seismic anisotropy in the upper 120 km only contributes a small portion of the observed splitting time (~30%), suggesting that the lower crust and upper mantle beneath the Qinling Orogenic Belt are relatively strong (Wu et al., 2022b). The observation of widespread high seismic velocity also offers support for the presence of a strong lower crust and upper mantle beneath the Qinling Orogenic Belt, suggesting low-viscosity crustal materials beneath the northeastern Tibetan Plateau not flowing eastward into the Qinling Orogenic Belt, especially to its eastern segment (Fig. 7; Guo and Chen, 2017; Jiang et al., 2016; Song et al., 2018). Together with the slow 3-D crustal movement, we suggest instead that there was a negligible migration of crustal materials along the Qinling Orogenic Belt especially to the east of 106°E longitude. However, the presence of eastward lower crustal flow beneath the western Qinling Orogenic Belt is consistent with a seismically low-velocity structure in the middle crust (Fig. 7B; Zhao et al., 2021), and other geophysical observations such as the EW-trending anomaly of tectonic stress and distinct loading ratios calculated from gravity measurements (Wang et al., 2018).

We present VLM and strain rate maps for the Ordos block and its vicinities based on continuous GNSS and new leveling observations. Considering that geodetic velocities are contaminated by non-tectonic deformation resulting from short-term surface mass changes and GIA, we estimate the related VLM of −0.4 mm/yr to +0.4 mm/yr from GRACE observations and 0.2–0.3 mm/yr from a GIA model. We then remove these non-tectonic VLM from geodetic observations and obtain an imaged VLM map by applying the RNI algorithm. Our VLM map reveals that uplift rates for most of the Ordos block are 1.5–2.0 mm/yr, while its northeast uplifts at the rate of 2–4 mm/yr and its northwest uplifts at the minimum rate of ~1 mm/yr. In the intersections between the southwestern Ordos block and the northeastern Tibetan Plateau, the crustal uplift rates are 1–3 mm/yr. The crust uplifts at a rate of 1–3 mm/yr in the Alxa block and Yinshan-Yanshan Mountains. Large-scale subsidence occurs in the eastern Qinling Orogenic Belt at a rate of ~−1 mm/yr, while the crustal uplift rate is 0.5–1.5 mm/yr in the western Qinling Orogenic Belt. Significant uplift of 2–4 mm/yr occurs in the Datong volcanic area adjacent to the northeastern Ordos block.

We use the combined constraints from these geodetic networks, seismic velocity, and anisotropy and other geophysical data to argue that uplift and extension of the northeast Ordos block and Datong volcanic area are consistent with mantle upwelling and lithospheric thinning. Combining with variations in the W-E component of horizontal GNSS velocities and exhibition of widespread high seismic velocity in the crust and upper mantle, we argue the migration of crustal materials along the Qinling Orogenic Belt especially east of 106°E longitude is negligible.

We thank science editor Wenjiao Xiao, associate editor Zheng-Xiang Li, and two anonymous reviewers for their insightful feedback. We thank Wenjun Zheng for his helpful discussions. We thank all the staff who participated in the leveling surveys. The Global Navigation Satellite System (GNSS) and leveling vertical velocities used in this paper are attributed to Li et al. (2023b). GNSS horizontal velocities are from Hao et al. (2021). Gravity Recovery and Climate Experiment solutions and GAC models were used to correct for vertical velocities for hydrological loading (Dobslaw et al., 2018). The glacial isostatic adjustment model ICE-6G_D (VM5a) is from Peltier et al. (2018). Seismic tomography and anisotropy data are from Wu et al. (2022a, 2022b) and Liu et al. (2021). Figures were plotted using the Generic Mapping Tools (Wessel et al., 2013; https://www.generic-mapping-tools.org/). This work was funded by the Shanxi Taiyuan Continental Rift Dynamics National Observation and Research Station (NORSTY20-01, NORSTY2022-02), National Natural Science Foundation of China (U22B6002, 42002234, and 41888101), National Key Research and Development Program of China (2017YFC1500100), Lhasa National Observation and Research Station of Geophysics (NORSLS20-03), Shanxi Key Research and Development program (2021SF2-03), and Postdoctoral Fellowship Program of China Postdoctoral Science Foundation (GZC20230043).