Crustal Deformation on the Northeastern Margin of the Tibetan Plateau from Continuous GPS Observations

Abstract

We installed 10 continuous Global Positioning System (GPS) stations on the northeast margin of the Tibetan Plateau at the end of 2012, in order to qualitatively investigate strain accumulation across the Liupanshan Fault (LPSF). We integrated our newly built stations with 48 other existing GPS stations to provide new insights into three-dimensional tectonic deformation. We employed white plus flicker noise model as a statistical model to obtain realistic velocities and corresponding uncertainties in the ITRF2014 and Ordos-fixed reference frame. The total velocity decrease from northwest to southeast in the Longxi Block (LXB) was 5.3 mm/yr within the range of 200 km west of the LPSF on the horizontal component. The first-order characteristic of the vertical crustal deformation was uplift for the northeastern margin of the Tibetan Plateau. The uplift rates in the LXB and the Ordos Block (ORB) were 1.0 and 2.0 mm/yr, respectively. We adopted an improved spherical wavelet algorithm to invert for multiscale strain rates and rotation rates. Multiscale strain rates showed a complex crustal deformation pattern. A significant clockwise rotation of about 30 nradians/yr (10⁻⁹ radians/year) was identified around the Dingxi. Localized strain accumulation was determined around the intersectional region between the Haiyuan Fault (HYF) and the LPSF. The deformation pattern across the LFPS was similar to that of the Longmengshan Fault (LMSF) before the 2008 Wenchuan M_S 8.0 earthquake. Furthermore, according to the distributed second invariant of strain rates at different spatial scale, strain partitioning has already spatially localized along the Xiaokou–Liupanshan–Longxian–Baoji fault belt (XLLBF). The tectonic deformation and localized strain buildup together with seismicity imply a high probability for a potential earthquake in this zone.

1. Introduction

Since the Middle Miocene, as a response to the India–Eurasia collision [1,2,3], the northeastern Tibetan Plateau has been undergoing lateral extrusion, which is partially accommodated by sinistral strike–slip movement along the HYF zone and crustal shortening toward the Tibet–Ordos transition zone of the northeastern margin along the Liupanshan area [4,5,6]. The northeastern part of the Tibetan Plateau is widely filled with Cenozoic sedimentary rocks [4,7,8]. To the north of the northeastern Tibetan Plateau, the Alxa Block (ALB) has undergone significant tectonic deformation since the Mesozoic and is a key region related not only to the Mesozoic deformation of Central Asia, but also to the Cenozoic development of the Tibetan Plateau [8,9]. The tectonic deformation west of the Liupanshan (LPS) is characterized by numerous active strike–slip and thrust faults [10]. The slip on these faults contributes to continuing uplift and seismic activity on the northeastern margin of the Tibetan Plateau [11]. East of the Liupanshan is the stable ORB, which tectonically belongs to the North China Block (NCB), a region that seems to share a common late Precambrian to Cenozoic geological history, but which has undergone regional extension during the Cenozoic [12,13,14].

The HYF and LPS are main boundary faults on the northeastern margin of the Tibetan Plateau, and are part of the most significant seismic hazard zone in the region [15]. Between these two main faults is the Xiaokou Fault (XKF), a left-lateral strike–slip fault which intersects the HYF and transfers displacement to the LPS [10,16]. East of the XKF is the Madongshan fold zone (MDFZ), located in the northern part of the Madongshan [10]. To the northeast of the HYF is the Xiangshan–Tianjingshan Fault (XTJSF), a left-lateral strike–slip fault with a large thrust component [3]. These two faults have similar slip features but different slip rates, which are 6.5 and 5.0 mm/yr, respectively, according to the geographical deflection of the Yellow River across the two faults [3]. On the south of the LPS, the Longxian–Baoji fault system (LXBJF) consists of four sub-faults, with the easternmost Qishan–Mazhao Fault (QSMZF) being the most active one, accommodating regional strain accumulation through strike–slip movement [17]. In this area, other major faults are the north fringe fault of the West Qinling Mountains (NWQF) and the East Kunlun Fault (EKF) [6]. Based on fault distribution, tectonic movement, and crustal deformation from geological and geodetic results, the study zone can be divided into five relatively rigid blocks, as shown in Figure 1 [18]. Besides the main boundary faults, there are some secondary faults within the blocks, which play an important role in the accommodation of strain distribution inside blocks. Three great earthquakes have occurred in this area with magnitudes larger than or equal to M_S 8.0, including the 1654 M_S 8.0 Tianshui earthquake, the 1879 M_S 8.0 Wenxian earthquake, and the 1920 M_S 8.6 Haiyuan earthquake.

Previous studies have shown that the velocity field derived from geodetic observations could be explained by the elastic block model, indicating gradually declining slip rates from west to east on the HYF, with an average slip rate being about 5.0 mm/yr [11,18,19,20,21]. However, due to the sparse distribution of GPS stations, slip rates on the LPS and LXBJF are poorly constrained. The slip rate gradually declines from west to east on the north fringe fault of the West Qinling Mountains, and tends to be almost zero on the eastern margin [22,23]. Previous studies were mostly based on campaign GPS observations, with a very limited amount of continuous GPS (cGPS) data.

Since 2012, with financial support from the China Earthquake Administration, more cGPS stations have been deployed in some key areas, such as the LPSF zone and the XKF zone. We collected the data from all 58 cGPS sites on the northeastern Tibetan Plateau margin. In this study, we calculated a new velocity field and then calculated the strain field based on an improved multiscale spherical wavelet algorithm. In Section 2, we elaborated the details of the cGPS data processing, the fitting of GPS time series. In Section 3, we described an improved algorithm for strain rate estimation based on wavelet analysis. In Section 4 and Section 5, we presented the results of the estimated velocity and strain rate fields and analyzed their tectonic implications, respectively. In Section 6, we analyzed strain accumulation on different spatial scales and discussed the relationship between horizontal and vertical movements. In Section 7, we presented a concise conclusion.

2. GPS Data and Data Processing

2.1. GPS Observations

To investigate the slip rate of the LPSF and analyze the seismic hazard of the fault, we installed eight continuously recording GPS stations across this fault. We installed the cGPS profile according to the principle that the inter-site distance gradually increases with the distance from the fault trace, with an average inter-site distance of 20 km, in order to capture fault activity and strain accumulation as precisely as possible. Five stations, LP01 to LP05, are located in the west of the fault, and three stations, LP06 to LP08, are located in the east. We also established two cGPS stations across the XKF, with an inter-station distance of 10 km, in order to determine its slip rate and study the transformation of the slip pattern between the HYF and LPSF. These 10 stations formed a new regional network, named the Institute of Earthquake Forecasting Network (IEFNET), which was established on October 2012, following the Crustal Movement Observation Network of China (CMONOC) reference station standards, and has been operating since April 2013 [24]. Additionally, we collected 48 other cGPS observations in this region, from two networks. Twenty-two of them belong to CMONOC, a high-rate GPS network with sampling rates of up to 50 Hz which consists of 260 continuous stations and 2000 campaign stations installed by the China Earthquake Administration. This network was established in 2009 to benefit the monitoring and investigation of tectonic movement and environmental change in China. The remaining 26 stations belong to a regional network (REGENT) which was installed by local agencies with earthquake, meteorology, and surveying and mapping, in order to increase the spatial density of GPS sites in key tectonic areas and provide positioning services [24]. This regional network has been operational since January 2010, following the CMONOC reference station standards. These stations played an important role in supplying local positioning service. However, they were not employed for the study of tectonic movement. The distribution of the 58 cGPS stations used in this study covers the main tectonic area with a mean inter-site distance of about 110 km, and most GPS stations are located close to major seismogenic zones (Figure 2).

In the following, we obtained 58 cGPS observations over the period of 2010–2018. The velocities of the 10 newly installed cGPS stations are firstly reported. Then, the application of 26 cGPS stations of REGENT and the 22 cGPS stations, which have longer observational time span compared with previous study, of CMONOC to the study of crustal and tectonic deformation is described. Newer and denser cGPS observations are adopted to acquire a three-dimensional velocity field in order to elaborate regional tectonic deformation. The data from the analyzed continuously recording stations spans more than 2.5 years, with majority of them covering more than 5.0 years, which is long enough to reliably estimate crustal velocities.

2.2. GPS Data Processing

The GAMIT/GLOBK package (version 10.6) was employed to estimate the daily solutions of station coordinates in a two-step approach [25,26]. In the first step, we analyzed doubly differenced GPS phase observations in daily session to estimate loosely constrained station coordinates together with orbital and Earth orientation parameters and their associated covariance matrices [27,28]. The sampling frequency was set to 30 s, and the cutoff of the elevation angle was 10° for the GPS raw data. The estimated orbital parameters included 15 elements for each satellite, while the Earth orientation parameters included pole position and rate along with the UT1 and its rate. The ionosphere-free combination was introduced in forming observation equations to eliminate the first-order ionosphere delay effects; however, the high-order ionosphere delay effects remained as errors in the observation equations, which can result in positioning errors of up to 15 mm in some areas with high electron density and taken into account in the precise positioning [29]. The second-order ionosphere effects were correlated with the total electron content and geomagnetic field [28]. Vertical Total Electron Content (VTEC) maps, from the Center for Orbit Determination in Europe (CODE), and the newly updated geomagnetic field model IGRF12 [30,31], from the international geomagnetic reference, were applied in order to mitigate their influences [24]. An automatic cleaning algorithm was applied to check raw data for better handling of poor data by repairing cycle slips and removing outliers [27]. IGS final satellite ephemeris was applied to calculate the position of GPS satellites [32]. The International Earth Rotation and Reference System Service (IERS) Bulletin B value was adopted for determining the initial values of the Earth orientation parameters. The IGS14 absolute phase center model was applied to correct the variation of phase center for satellites and receivers. The Vienna Mapping Function 1 (VMF1) mapping function was employed to model the troposphere delay effects by means of the dry and wet parts [33], with atmosphere pressure and temperature parameters extracted from station meteorological files or GPT2 grid files, to improve the precision of estimated station position on the vertical component [34]. The wet part of the troposphere delay effects was set as unknown hourly parameters in the final estimate of position coordinates, together with two daily gradient parameters in the east and north components, respectively. The ocean tidal loading correction was taken into account according to the FES2004 model [35], and the effects of polar tides and solid Earth tides were corrected according to the IERS2003 model [36]. The atmospheric loading correction was applied in the CM reference frame (center of mass of the solid Earth plus fluid loads) according to the IRES 2003 model. The solar radiation pressure was estimated using the Berne model [37]. In this step, we obtained the variance-covariance matrix, which includes bias-free and bias-fixed solutions. We also collected data from other CMONOC cGPS stations and about 80 ITRF2014 reference stations distributed all over the world to realize the reasonable reference frame [38]. The position time series plots of 58 cGPS stations relative to the ITRF2014 reference frame (Figure S1) are listed in the supplementary materials.

In the second step, we used the daily bias-fixed loosely constrained solutions of the estimated parameters and their covariance matrices as quasi-observation in GLOBK, which is a smoothing Kalman filter, to resolve position time series in the ITRF2014 by means of network adjustments. We constrained all the priori coordinates to 0.1 m to minimize the effects of unmodeled site position biases. The ITRF2014 reference frame was realized by minimal constraining of the 80 ITRF2014 core stations by a seven-parameter Helmert transformation, which used three orientation parameters, three translation parameters, and one scale parameter. We adopted four iterations to eliminate bad sites and to compute station weights for the reference frame stabilization.

2.3. Time Series Analysis

In order to fit coordinate time series, we employed a functional model that encompasses a constant velocity, an annual and a semiannual sinusoidal variation term, and possible offset parameters for each component of each station. The functional model can be written, after [39], as:

$$D_k\left(t_i\right)=D_k\left(t_0\right)+v_k{t}_i+{{\displaystyle \sum }}_{m=1}^2\left[a_{m}sin\left(2m\pi t_i\right)+b_{m}cos\left(2m\pi t_i\right)\right]+{{\displaystyle \sum }}_{j=0}^{j_0}{o}_{j}H\left(t_i-t_j\right)+e_i$$

(1)

where $k=e,n,u$ for the east, north, and up components, respectively, $D_k\left(t_0\right)$ is the station position at time $t_0$, $v_k$ is the linear rate of component $k$, $a_1$ and $b_1$ are the sine and cosine amplitudes at annual period, $a_2$ and $b_2$ are the sine and cosine amplitudes at semi-annual period, $H\left(t_i-t_j\right)$ is the Heaviside step function which is equal to 1 for $t_i\ge t_j$ and 0 otherwise, $o_j$ is the offset value at time $t_j$, and $e_i$ is the corresponding noise in the position time series.

Outliers were removed by using the interquartile range algorithm described by [40]. A noise model of time series can affect the estimation of the velocity field, especially the uncertainty of the estimated trend [41]. The white noise assumption that measurement error is random and uncorrelated from one epoch to its neighbor is improper for time series analysis. Time-correlated noise was confirmed to be present in the GPS position time series [42]. We employed the Hector software to estimate the power index of the noise model, which includes both white and power index models [43]. The results of the noise analysis show that the corresponding noise model is the white plus flicker noise model [41]. Therefore, to precisely estimate the velocity field, we applied a white plus flicker noise model to form the covariance matrix for the least squares adjustment of secular trends and acquired realistic uncertainties for the velocity estimates.

3. Strain Rate Calculation

Strain rate is independent of GPS reference frame, and provides a better way to characterize crustal deformation and fault strain accumulation [27]. We employed an improved spherical wavelet method, which can decompose total strain into various spatial scales compared to the method of Tape et al. [44], to acquire multiscale strain rates for an irregular distribution of GPS stations. The spherical wavelet theorem is evolved from infinite plane space to finite spherical space. The spherical wavelet basis function is derived from their mother function by means of the translation and scaling. We employed the difference of Gaussian (DOG) function to construct the wavelet basic function. The spatial distribution of the wavelet basic function was heterogeneous due to the spatial irregularity of the GPS stations. Figure 3 shows the shape of the normalized DOG wavelet basis function on different decomposition scales against the distance between the observational point and the location of the wavelet base function.

We determined the minimum decomposition scale 3 according to the radius of the study area. The maximum decomposition scale was determined by the density of the cGPS stations. We obtained the maximum decomposition scale by searching for the area with the smallest radius, which included three cGPS stations, at the location of the wavelet basic function. The decomposition scales 6, 7, and 8 correspond to spatial resolutions of 174.4, 87.2, and 43.6 km, respectively. Figure 4 displays the distribution of spatial resolution; most of the area attained a spatial resolution of 87.2 km. Therefore, we determined scales 3–7 as the decomposition scales to estimate the strain rates.

When we constructed the spherical wavelet basic function at various scales, we were able to build the relationship between the GPS-derived velocities and the coefficients of the wavelet basic function, which can be written as

$$v\left(\theta ,\phi \right)={{\displaystyle \sum }}_{k=1}^M\left[f_{k_{ew}}·g_k\left(\theta ,\phi \right)·\hat{\theta }+f_{k_{\mathrm{ns}}}·g_k\left(\theta ,\phi \right)·\hat{\phi }\right]$$

(2)

where $\theta$ and $\phi$ are the colatitude and longitude of the specified GPS site, respectively, $\widehat{\theta }$ is the unit vector of the north component and $\widehat{\phi }$ is the unit vector of the east component, $v\left(\theta ,\phi \right)$ is the observational velocity vector of the GPS site at $\left(\theta ,\phi \right)$, M is the total number of wavelet basis functions, $g_k\left(\theta ,\phi \right)$ is the k-th DOG wavelet basis function, and $f_{k_{ew}}$ and $f_{k_{ns}}$ are the k-th estimated parameters on the east and north components, respectively. As the DOG wavelet basis functions are derived from the translation and scaling of their mother functions, they are non-orthogonal, which results in rank deficiency for the coefficient matrix of the normal equation obtained from the observed equation. The Tikhonov criterion was employed to balance residual and smoothness to address this issue. Figure 5 shows regularization parameter curves in the fitting of the GPS observational velocities. The optimal regularization can be obtained with the minimum value of the ordinary cross-validation (OCV) function for the north and east components, respectively. Figure 6 displays the observed and modeled velocities with arrows of different colors. Our spherical wavelet model reconciles the observed GPS velocities well: 86% (50 of 58) of the residual velocity components are less than 1.0 mm/yr, indicating that our model can effectively reconstruct observed velocities.

The first derivatives are calculated for the wavelet basis functions, and the velocity gradient tensor at a specified point can be expressed as

$$\left[\begin{array}{c}{\epsilon }_{\theta }\\ {\epsilon }_{\theta \phi }\\ {\epsilon }_{\phi }\\ {\epsilon }_{\phi \theta }\end{array}\right]=\left[\begin{array}{cc}\frac{\partial G\left(\theta ,\phi \right)}{\partial \theta }& 0\\ \begin{array}{c}0\\ \frac{\partial G\left(\theta ,\phi \right)}{\partial \phi }\\ 0\end{array}& \begin{array}{c}\frac{\partial G\left(\theta ,\phi \right)}{\partial \theta }\\ 0\\ \frac{\partial G\left(\theta ,\phi \right)}{\partial \phi }\end{array}\end{array}\right]·\left[\begin{array}{c}{f}_{ns}\\ f_{ew}\end{array}\right]$$

(3)

where ${\epsilon }_{\theta }$, ${\epsilon }_{\theta \phi }$, ${\epsilon }_{\phi }$, and ${\epsilon }_{\phi \theta }$ are the components of the velocity gradient tensor, $\frac{\partial G\left(\theta ,\phi \right)}{\partial \theta }$ and $\frac{\partial G\left(\theta ,\phi \right)}{\partial \phi }$ are first-order partial differentials of the wavelet basis functions in the north and east components, respectively, $f_{ns}$ and $f_{ew}$ are estimated parameters for the north and east components, respectively. We constructed a uniform grid and used the wavelet basic function to calculate the two-dimensional velocity gradient tensor. The strain rate tensor and rotation rate tensor at a specified site can be obtained with the velocity gradient tensor. The strain rate tensor can be further used to calculate various strain components, such as second invariant, principal strain, dilatation, and maximum shear strain [45].

4. GPS-Derived Three-Dimensional Velocities

The ITRF2014 velocity solution, accounting for annual periodical signal, semi-annual periodical signal, offsets, and time-correlated noise contents, was derived through the combination of individual daily solutions. In order to highlight crustal deformation, we transformed the ITRF2014-fixed velocity field into a local reference frame by the Euler angular vector. We defined the Ordos-fixed reference frame through the Euler angular vector from five stations (YANC, QYSC, QYTQ, QYNX, and SNXY) over the interior of the ORB. The estimated angular rotation rate of the ORB relative to ITRF2014 was 0.3980 ± 0.0319°/Myr (Million year), with a pole at 77.4059 ± 0.7839°N and 179.6720 ± 1.5001°W. To validate the applicability of the estimated Euler angular parameters, we performed a comparison between the predicted velocities and the observed velocities. The mean discrepancy between the two was 0.9 mm/yr, which is on the same level as the uncertainties in the GPS velocities. The station velocities and their uncertainties, both relative to the ITRF2014 and to the Ordos-fixed reference frame (Table S1), are listed in the supplementary materials.

The lateral expansion in the northeastern Tibetan Plateau was obstructed by the ORB, resulting in large-scale crustal shortening, as can be seen from the velocity field shown in Figure 6. The pattern and magnitude of our velocities are in general agreement with previous results derived from both campaign-mode and cGPS stations [18]. However, our cGPS station profile across the LPSF provided a higher resolution kinematic pattern. The determined horizontal movement on the northeastern margin of the Tibetan Plateau was about 8.8 ± 0.7 mm/yr in the northeastern direction, which was an average from three stations in the southwestern corner (stations QHMQ, GSMA, and QHTR). The magnitude of the horizontal movement gradually decreased to 3.2 ± 0.6 mm/yr with an azimuth of 38° at station GSDX which was located at the center of the LXB and far away from the boundary faults. Furthermore, the closer to the LPS the GPS stations were, the smaller their velocities were. This indicates that the deformation was localized near the LPS because of the resistance of the relatively stable ORB [10]. From the LXB to the ALB, the magnitude of velocities decreased significantly, and the direction of velocities changed from northeast to north, indicating that the two boundary faults (the HYF and the XTJSF) are both presently active.

The vertical movement at most GPS sites was characterized by uplift, except for several sites which were situated on localized sedimentation layer (Figure 7). The average uplift rate (relative to the ITRF2014 reference frame) was 1.6 mm/yr in the northeastern Tibetan Plateau. Using the leveling data which have been observed repeatedly between 1970 and 2012, Li et al. [18] estimated the vertical motion in the northeastern Tibetan Plateau with a maximum uplift rate of 7.0 mm/yr. Using the same data as Li et al. [18], we obtained the average uplift rate of 3.7 mm/yr in this area. The discrepancy between our GPS and leveling results was about 2.0 mm/yr, which may be attributed partly to the different time spans of observation. The change in vertical velocities across the HYF and the NWQF was small, while significant change was observed across the LPS. The two faults were both dominated by lateral slip, which generated horizontal crustal movement. However the latter was dominated by extrusion slip, which generated deformation horizontally and vertically [46]. There were a few stations in the ALB, showing small but consistent uplift which is insignificant relative to their uncertainties. The vertical movements in the Qinling Block (QLB) were quite diffuse, indicating a heterogeneous tectonic deformation pattern. The vertical movements of the LXB was characterized by uplift, with an average rate of 1.0 ± 0.9 mm/yr, and a gradual increase from northwest to southeast. Although the average rate was insignificant relative to its uncertainty, the uniform uplift feature at the majority of GPS stations can be observed at the interior of LXB. Five stations in the ORB were uplifting, with an average rate of 2.0 ± 1.1 mm/yr and a maximum rate of 3.7 ± 1.5 mm/yr, which is a larger than that of LXB.

To qualitatively analyze the deformation process from the LXB to the ORB, we incorporated our velocities with the campaign-mode GPS velocities of Ye et al. [47], which adopted the same data processing strategy with this study and the observational period being 1999–2017. We took a velocity profile to identify the change of velocities (Figure 8). A noticeable velocity decrease from northwest to southeast interior of the LXB was identified on the horizontal component. We estimated a reduced velocity of 5.3 ± 0.7 mm/yr within the range of 200 km west of the LPSF. However, no significant differential motions were identified in the ORB. This deformation pattern across the LPS is similar with that observed across the LMSF before the 2008 Wenchuan M_S 8.0 earthquake [48,49].

5. GPS-Derived Strain Rates

By performing a multiscale spherical wavelet strain calculation, we obtained a two-dimensional strain rate tensor and rotation rate tensor. Dilatation rate is the first invariant of the strain tensor, which is equal to the sum of the maximum and minimum strain rate vectors, and is independent of coordinate system, indicating the areal change of the tiny deformation unit. The formula of the second invariant is $\sqrt{e_{ee}^2+e_{nn}^2+2e_{ne}^2}$, which reflects the magnitude of total strain rates. Maximum shear strain rate is equal to the difference between the maximum and minimum strain rate and is related to regional tectonic shear deformation. By these strain rate components, we can determine comprehensive horizontal surface deformation. The rotation rate is the antisymmetric part of the velocity gradient tensor, representing the rotation feature of the very small deformation unit.

Figure 9 shows the second invariant and corresponding uncertainties. The maximum magnitude of the uncertainties is about 8 nstrain/yr (10⁻⁹/year), much smaller than the second invariant at the significant strain partitioning area. The remarkable strain rates were located in the main fault zone and tectonic transformation area. In the intersection area between the HYF and Lenglongling Fault (LLLF), a wide region of strain rates were obtained with a maximum magnitude of about 45 nstrain/yr. Along the XKF, located at the intersection region between the HYF and LPSF, localized strain rates were obtained with a magnitude of about 50 nstrain/yr. Along the NWQF and the EKF, significant strain rates were also obtained with various magnitudes. Additionally, a large strain rate (about 70 nstrain/yr) was obtained in the Weihe Basin (WHB), possibly due to the shortening between stations SNHX and XIAA. However, these two stations are undergoing significant subsidence, which may have affected the estimation of the horizontal velocities.

The principal strain rates showed a wide range of crustal deformation behaviors in the northeastern margin of the Tibetan Plateau (Figure 10). Across the HYF, XTJSF, and LPSF, the magnitude of the strain rate vectors gradually decreased, and tended to zero in the ALB and the ORB. Li et al. [18] also obtained negligible strain at the interior of the ALB and the ORB using combined cGPS and campaign-mode observation on the period of 1999–2017. The HYF and LPSF, as the main boundary faults in this region, played an important role in the adjustment of the strain distribution [6]. Several great earthquakes have occurred along these two boundary faults, including the 1622 Guyuan M_S 7.0 earthquake and the 1920 Haiyuan M_S 8.6 earthquake.

The determined orientation of the principal strain rates along the HYF is consistent with the pattern of pure strike–slip fault [46], the axes of the strain rates being oriented close to 45° to the fault trace, and even swinging around to maintain that angle as the fault zone changes strike [46]. However, the magnitudes of the strain rates decrease from west to east. Along the XKF, the magnitude of extension reduces to about 40 nstrain/yr, and the azimuth of extension is different from that of other segments of the HYF. Along the LPSF, the magnitude of shortening was small, indicating insignificant strain accumulation near the fault and significant strain accumulation in a large area, quite similar to what was observed in the LMSF before the 2008 Wenchuan earthquake [48,49].

Significant positive dilatation (about 60 nstrain/yr) was observed around the LXBJF, with the direction of extension being perpendicular to the fault trace. This indicates that the effective slipping from the four faults of the fault is a mixture of normal faulting and strike–slip faulting. Using geological data and earthquake focal mechanism solutions, Li. et al. [17] concluded that the late Quaternary activity of the QSMZF is left-lateral with normal component, with a strike–slip rate of about 0.5–1.0 mm/yr. In the LXB, the negative dilatation was the first-order feature, while a positive dilatation belt with a north-northeast direction was mapped in the central part of LXB. The Zhuanglanghe–Maxianshan Fault (ZLMXF) is the main active fault in this positive dilatation belt. The Zhuanglanghe Fault (ZLHF) is dominated by thrust slipping at a rate of 3.9 mm/yr as estimated from previous GPS observations [50]. Strain rates are diffuse in the QLB, indicating a complicated pattern of regional tectonic deformation.

Large shear strain rates were obtained at the location of large second invariant (Figure 11). Along the HYF, the magnitude of the shear rates was varied, indicating spatial differences in the slip rates and the locking depth [19,20]. Based on a elastic model, the inversion result from geodetic data showed that on the western part of the HYF, from the Tuolaishan segment (TLS) to the Laohushan segment (LHS), the locking depth of the fault was 15 km, while on the eastern segment the locking depth was shallower, at about 5 km [19,20].

A significant clockwise rotation can be observed around the Dingxi (Figure 12). In this area, we observed a clockwise rotation of 5–7°/Myr, which was in agreement with the result of Li et al. [18]. Paleomagnetic data [51] indicated a significant clockwise rotation of 20–30°/Myr near Lanzhou. Our estimates of rotation rates based on cGPS data show that this region was still undergoing a clockwise rotation. However, the magnitude of this rotation was smaller by a factor of four compared to the paleomagnetically-determined rotation. Along the HYF, we observed an anticlockwise rotation belt with an average rate of 10 nradians/yr. If we used the absolute value of the rotation rate to map the magnitude which was adopted by Li et al. [18], we can obtain the similar results with them, indicating that our rotation rate along the Haiyuan fault agreed with the result of Li et al. [18].

6. Discussion

6.1. Strain Buildup at Various Spatial Scales

GPS strain rate estimated by multiscale spherical wavelet approach can be decomposed into different spatial scales which is dominated by scales of structures. We decomposed total strain rates into scales 3–6 and scale 7, as shown in Figure 13. At scales 3–6, several areas with significant strain accumulation can be identified. A distinct strain rate belt can be observed along the XLLBF at scale 7. This indicates that strain rate, along the HYF, NWQF, and EKF, was significant at the large spatial scale. However, the strain rate was spatially localized along the XLLBF. According to the seismicity and tectonic deformation near the XLLBF, Du et al. [52] concluded this belt has a high probability for a future earthquake with the magnitude of 7.0. The spatially localized strain suggested a high seismic probability along the XLLBF.

6.2. Strain Partitioning at Different Fault Activity Styles

On the northeastern margin of the Tibetan Plateau, except for several main faults, there were numerous strike and thrust faults, indicating complicated historical tectonic deformation. Geodetic data from GPS and InSAR has shown that the left-lateral slip gradually decreased from west to east across the Haiyuan Fault, transformed to shortening on the Liupanshan Fault, and led to multiple thrust fold belts and secondary strike faults [53]. The heterogeneous distribution on space is related to the active styles of the faults. The pattern of strain rates was consistent with present-day geological structures.

6.3. Decoupling between Uplift of the Ordos Block and Horizontal Extrusion on the Northeastern Tibetan Plateau

GPS-derived horizontal velocities revealed that there was a large area of crustal shortening deformation on the northeastern margin of the Tibetan Plateau, while the shortening was insignificant from the LPS to the ORB. Due to the limited distribution of GPS stations, we cannot precisely delineate the specific boundary of the deformation. However, we can confirm that the LPS was at the forefront of the deformation zone, and the ORB was quite stable with no relative deformation at its interior. The structures of the LPS area were characterized by numerous active thrust and strike–slip faults that suggest thin-skinned deformation there [8,16].

A deep seismic reflection line across the Tibet–Ordos transition zone revealed a complex crustal structure underlying the transition zone and adjacent areas [14]. The thrust deformation on the northeastern margin of the Tibetan Plateau appeared to stop just east of the Xiaoguanshan Fault (XGSF), marking the easternmost edge of the strata [16], indicating that the XGSF is at the frontier of the thrust deformation of the Tibetan Plateau (Figure 14). The overall movement of the ORB on the eastern side of the XGSF was decoupled with the deformation mode of the Tibetan Plateau.

According to previous results of tectonic model across the LPSF [2,54], we could assume that the northeastern part of the Tibetan Plateau should have a larger uplift rate compared with the ORB. However, our results indicated an inconsistent pattern. The first-order vertical crustal deformation is characterized by uplift in the northeastern margin of the Tibetan Plateau. The LXB is uplifting at an average rate of 1.0 mm/yr, while the ORB is uplifting at a rate of 2.0 mm/yr, indicating that the lateral expansion of the Tibetan Plateau had insignificant effect on the uplift pattern of the stable ORB. The vertical crustal deformation derived from GPS observations in NCB (except for the North China Plain area but along the Yanshan and Taihangshan orogenic belts) was characterized by uplift with a rate of 2.0–3.0 mm/yr, similar to the value obtained for the ORB. Therefore, we can infer that the uplift pattern of the ORB was coupled with the NCB but decoupled with the deformation of the northeastern Tibetan Plateau. The similarity of the stratigraphic sequence present in the ORB to that of the rest of the larger NCB indicated that the entire region has undergone a common depositional and tectonic evolution [16].

7. Conclusions

Using GPS data from 58 cGPS stations over the period of 2010–2018, we presented an updated three-dimensional velocity field on the northeastern margin of the Tibetan Plateau and its surrounding areas. We employed a color noise model to represent the statistical characteristic of coordinate time series, and solved for reasonable velocity uncertainties. We adopted an improved spherical wavelet algorithm to estimate multiscale strain rates and rotation rates and analyzed the characteristic of strain accumulation. Some conclusions can be drawn as follows:

The GPS-derived velocities showed a noticeable decrease in magnitude and a change in orientation from the southwestern QLB to the western edge of LPS. The total velocity decrease from the northwest to southeast in the LXB, on the horizontal component, was 5.3 mm/yr within the range of 200 km west of the LPSF. The direction and magnitude of the velocities are diffuse near the boundary of the LPSF. The tectonic deformation and localized strain buildup, together with seismicity jointly imply a high probability for a potential large earthquake at this belt.

The vertical crustal deformation was characterized by uplift in the northeastern margin of the Tibetan Plateau. The LXB and the ORB are both uplifting, with rates of 1.0 and 2.0 mm/yr, respectively. With the new GPS observations, we derived a more rapid uplift rate of the ORB compared to the LXB, refuting the detached structural model of the northeastern margin of the Tibetan Plateau.

Multiscale components of strain rate tensor show a complex deformation pattern on the northeastern margin of Tibet Plateau. Positive dilation was observed around the LXBJF. A significant clockwise rotation was obtained around the Dingxi. The HYF and LPSF played an important role in controlling regional strain partitioning. The distribution of the principal strain rates along the HYF is in agreement with the pattern of strike–slip faults. However, along the LPSF, the direction of principal strain rates is in disorder, indicating that the LPSF is undergoing complicated activities. Strain partitioning demonstrated various spatial scales in different regions, such as a wide range strain near the LLLF and a localized strain near the XKF. The strain rate was spatially localized along the XLLBF. The tectonic deformation and localized strain buildup together with seismicity imply a high probability for a potential large earthquake at this belt.