A Crustal Deformation Pattern on the Northeastern Margin of the Tibetan Plateau Derived from GPS Observations

: The northeastern margin is a natural experimental ﬁeld for studying crustal extrusion and expansion mechanisms. The accurate crustal deformation pattern is a key point in the analysis of regional deformation mechanisms and seismic hazard research and judgment. In this paper, the present-day GPS velocity ﬁeld on the northeastern margin of the Tibetan Plateau was obtained from encrypted GPS observations around the Haiyuan–Liupanshan fault zone, combined with GPS observations on the northeastern margin of the Tibetan Plateau from 2010 to 2020. Firstly, we divided the study area into three relatively independent blocks: the ORDOS block, Alxa block, and Lanzhou block; secondly, the accurate fault distribution of the Haiyuan–Liupanshan fault zone was taken into account to obtain the optimal inversion model; ﬁnally, using the block and fault back-slip dislocation model, the inversion obtained the slip rate distribution, locking depth, and slip deﬁcit rate of each fault. The results indicate that the Laohushan Fault and Haiyuan Fault are dominated by the left-lateral strike-slip, while the Liupanshan Fault is dominated by the thrust dip-slip, and the Guguan–Baoji Fault has both left-lateral strike-slip and thrust dip-slip components. The maximum locking depths of the Laohushan Fault, Haiyuan Fault, Liupanshan Fault, and Guguan–Baoji Fault are 5 km, 13 km, 15 km, and 10 km, respectively, and the locking of the Haiyuan Fault is strong in the middle section and weak in the eastern and western section. The Haiyuan Fault is still in the post-earthquake stress adjustment stage. The slip deﬁcit rate decays from 3.6 mm/yr to 1.8 mm/yr from west to east along the fault zone. Combined with geological and historical seismic data, the results suggest that the mid-long-term seismic risk in the Liupanshan Fault is high.


Introduction
The northeastern margin of the Tibetan Plateau is the frontal region of the uplift development that has expanded into the continental interior, the newest and forming part of the plateau, and a natural experimental field to study the continental interior deformation of China. The extrusion and deformation of the crust have resulted in the formation of complex tectonic zones and seismic activity zones; a series of historical earthquakes occurred on the active faults in the region and margins (Figure 1), such as the 1920 Haiyuan M8.5 earthquake and the 1927 Gulang M8.0 earthquake, which occurred on the HYF [1] ( Figure 1). According to the back-slip theory, a fault can accumulate strain accumulation near the fault surface due to the relative motion of the footwall and hangingwall blocks when the fault is locked and produce a slip deficit phenomenon [2]. Ref. [3] obtained the crustal motions from the Caucasus mountains to the Adriatic Sea and north-south from the southern edge of the Eurasian plate to the northern edge of the African plate at 189 GPS sites. Ref. [4] analyzed the complex kinematic and deformation fields in the North Aegean Sea and adjacent regions using GPS observations from 1993 to 2009. Ref. [5] studied the crustal deformation in the Northern Aegean Sea using GPS data. Ref. [6] used the GPS velocity field to explore the detailed complexity of Aegean tectonics. Therefore, Many scholars have conducted large numbers of studies on this region using three methods: geological methods, GPS, and InSAR. However, due to different spatial resolutions of data, station density, observation duration, and model selection, the slip rates and deformation characteristics show significant discrepancies (Table 1).  [7] 3.3-9.2 Geological results [8] Zhongwei section: 8-10, reduced eastward to 4-6 [9] [10] 5.8 [11] 4.5 ± 1.1 Many scholars have conducted large numbers of studies on this region using three methods: geological methods, GPS, and InSAR. However, due to different spatial resolutions of data, station density, observation duration, and model selection, the slip rates and deformation characteristics show significant discrepancies (Table 1).
In summary, the slip rates of the HYF obtained by geological methods are 3-10 mm/yr, with an average rate of 5-6 mm/yr; the slip rates of the HYF obtained derived from GPS observations and block deformation analysis are 2-9 mm/yr; and the slip rates of the LPSF are 1-3.5 mm/yr. Derived from InSAR technology, the average rate of the HYF is 3-5 mm/yr, and the LPSF is about 3 mm/yr.
In order to obtain a more accurate crustal deformation pattern of the faults on the northeastern corner of the Tibetan Plateau, this paper integrated the observations from encrypted GPS stations around the Haiyuan-Liupanshan fault zone, with the GPS observations on the northeastern margin of the Tibetan Plateau from 2010 to 2020, based on the block-fault back-slip dislocation model to obtain the locking depth, slip rate deficit distribution, and motion characteristics of the faults.  Figure 1 shows the distribution of active faults and earthquakes. The northeastern margin of the Tibetan Plateau is located on the plate that is extruded eastward by the collisional contraction between the Indian and Eurasian plates at the orientation of N20 • E [25]. The Haiyuan-Liupanshan fault zone is the boundary zone between the Lanzhou block and Alxa block on the northeastern corner of the Tibetan Plateau, which occupies a significant function for studying the mechanism of contraction and expansion of the Tibetan Plateau into the interior of the Chinese continent [26]. The HYF and LHSF are the junction faults between the Lanzhou and Alxa blocks, which are NWW, with the top of the arc convex to the NE (Figure 1), and its motion since the Late Cenozoic has been successively thrust left-lateral and strike-slip, and this strike-slip displacement has been absorbed by the LPSF to the east with east-west crustal shortening [27]; and earthquakes with a magnitude of more than seven have not occurred since 1920. The LPSF is located in the boundary zone between the Lanzhou block and the ORDOS block, which is a thrust fault and has left-lateral strike-slip characteristics [26]. Three earthquakes with a magnitude of above seven have occurred from 1200 to 1700 [28], and the earthquake risk should not be underestimated in the LPSF. The GGBJF is connected with the LPSF to the north and the northern edge of the West Qinling Fault in the south, which trends NNW ( Figure 1). It can be considered as the southern extension of the LPSF, but its movement pattern is completely different from that, so its movement characteristics also deserve attention. Figure S1 shows the earthquake distribution with magnitudes of more than 3.0, and the time span is from 1976 to 2020. The data is from the National Earthquake Data Center of China (https://data.earthquake.cn/ accessed on 24 September 2022). We can find that many earthquakes have occurred in the northeast margin of the Tibetan Plateau, and the focal depths of most earthquakes are less than 25 km. networks. One is the China Continental Environmental Tectonic Monitoring Network (CMONOC), which includes 459 GPS stations in our research area, and the observation time span is from 2010 to 2020. In addition, to better constrain the near-field deformation of the HYF and LPSF, we installed the Institute of Earthquake Forecasting Seismic Forecasting Network (IEFNET) in 2013. The IEFNET includes ten GPS stations, eight stations with an average spacing of 20 km across the middle of LPSF, and another two stations with a spacing of 10 km across the Xiaokou Fault. The distribution of GPS stations of different observation networks can be found in Figure S1. Richer GPS observations were used to obtain the fault slip rate, accurately capture the fault activity characteristics, and discuss the transition of fault activity patterns in the intersection region of the HYF and LPSF. The station distribution covers the main faults, with most stations located near major seismic zones, providing data assurance for obtaining a more accurate regional three-dimensional GNSS velocity field and data basis for obtaining a more accurate crustal deformation pattern in this region. Figure 2 shows the GPS stations and corresponding horizontal motions with respect to the ORDOS reference block. The GPS observations were obtained from two observation networks. One is the China Continental Environmental Tectonic Monitoring Network (CMONOC), which includes 459 GPS stations in our research area, and the observation time span is from 2010 to 2020. In addition, to better constrain the near-field deformation of the HYF and LPSF, we installed the Institute of Earthquake Forecasting Seismic Forecasting Network (IEFNET) in 2013. The IEFNET includes ten GPS stations, eight stations with an average spacing of 20 km across the middle of LPSF, and another two stations with a spacing of 10 km across the Xiaokou Fault. The distribution of GPS stations of different observation networks can be found in Figure S1. Richer GPS observations were used to obtain the fault slip rate, accurately capture the fault activity characteristics, and discuss the transition of fault activity patterns in the intersection region of the HYF and LPSF. The station distribution covers the main faults, with most stations located near major seismic zones, providing data assurance for obtaining a more accurate regional three-dimensional GNSS velocity field and data basis for obtaining a more accurate crustal deformation pattern in this region. GAMIT/GLOBK software, version 10.6, was employed to estimate the daily solutions of station coordinates [29,30]. The ORDOS reference frame was defined by five stable GPS stations, and the ITRF2014 velocity field was transformed into the ORDOS reference frame. The specific strategy of the GPS data process can be found in the literature [31]. Figure 2 shows the GPS velocity field with respect to the ORDOS block regional framework.

Block Division Scheme
According to the previous division model of blocks [32], the northeastern margin of the Tibetan Plateau is divided into the ORDOS block, Alxa block, and Lanzhou block with the boundaries of the LHSF, HYF, LPSF, and GGBJF, and block boundaries coincide with the active faults ( Figure 3). GAMIT/GLOBK software, version 10.6, was employed to estimate the daily solutions of station coordinates [29,30]. The ORDOS reference frame was defined by five stable GPS stations, and the ITRF2014 velocity field was transformed into the ORDOS reference frame. The specific strategy of the GPS data process can be found in the literature [31]. Figure 2 shows the GPS velocity field with respect to the ORDOS block regional framework.

Fault Geometry Model
The LHSF is the boundary zone between the Alxa block and the Lanzhou block, and it intersects with the west section of the HYF; the GGBJF is the boundary zone between the ORDOS block and the Lanzhou block, which is an extensional fault of the LPSF and intersects the south section of it. Both the LHSF and GGBJF are highly active faults and have frequent historical earthquakes ( Figure 1) and have a strong influence on the near-field observations of the Haiyuan-Liupanshan fault zone. In order to improve the accuracy of model inversion and better fit the GPS observations, the LHSF and GGBJF were added to the model, and the elastic deformation caused by their locking was considered. The LHSF and HYF are the boundary zone between the Alxa block and the Lanzhou block, and the LPSF and GGBJF are the boundary zone between the ORDOS block and the Lanzhou block.
Accurate determination of the geometric model of the fault could make the model better fit the observed GPS velocity field, thus reducing the residuals of the velocity field to obtain better results. When modeling the fault surface (Figure 4), the dip angle of the LHSF and HYF was set to 60° with an SSW tendency, and the dip angle of the LPSF and GGBJF was set to 45° with an SW tendency. The fault was set to be completely locked ( = 1) at 0.1 km and free-slipped ( = 0) below 25 km. The fault node parameterization type was set to a one-dimensional depth profile type in which the node value is a function of depth [33], which was expounded in detail in Section 6.2. In the inversion, two surface nodes were set along the strike, and six nodes were set along the depth direction for the LHSF. The HYF had eight surface nodes and six depth nodes, the LPSF had six surface nodes and six depth nodes, and the GGBJF had three surface nodes and six depth nodes. The fault surface node meshing was based on the actual surface traces of the active fault. Additionally, by setting surface and depth nodes along the strike and dip of the fault, the fault is constructed into a fault surface composed of many small grids. The locking fraction at each small grid is iteratively calculated using an inversion algorithm, thereby obtaining a continuous distribution of the locking fraction of the entire fault surface. The number of nodes on the fault surface and the depth are the optimal values obtained through extensive trial calculations. All these settings are aimed at obtaining more accurate crustal deformation characteristics.

Fault Geometry Model
The LHSF is the boundary zone between the Alxa block and the Lanzhou block, and it intersects with the west section of the HYF; the GGBJF is the boundary zone between the ORDOS block and the Lanzhou block, which is an extensional fault of the LPSF and intersects the south section of it. Both the LHSF and GGBJF are highly active faults and have frequent historical earthquakes ( Figure 1) and have a strong influence on the near-field observations of the Haiyuan-Liupanshan fault zone. In order to improve the accuracy of model inversion and better fit the GPS observations, the LHSF and GGBJF were added to the model, and the elastic deformation caused by their locking was considered. The LHSF and HYF are the boundary zone between the Alxa block and the Lanzhou block, and the LPSF and GGBJF are the boundary zone between the ORDOS block and the Lanzhou block.
Accurate determination of the geometric model of the fault could make the model better fit the observed GPS velocity field, thus reducing the residuals of the velocity field to obtain better results. When modeling the fault surface ( Figure 4), the dip angle of the LHSF and HYF was set to 60 • with an SSW tendency, and the dip angle of the LPSF and GGBJF was set to 45 • with an SW tendency. The fault was set to be completely locked (φ = 1) at 0.1 km and free-slipped (φ = 0) below 25 km. The fault node parameterization type was set to a one-dimensional depth profile type in which the node value is a function of depth [33], which was expounded in detail in Section 6.2. In the inversion, two surface nodes were set along the strike, and six nodes were set along the depth direction for the LHSF. The HYF had eight surface nodes and six depth nodes, the LPSF had six surface nodes and six depth nodes, and the GGBJF had three surface nodes and six depth nodes. The fault surface node meshing was based on the actual surface traces of the active fault. Additionally, by setting surface and depth nodes along the strike and dip of the fault, the fault is constructed into a fault surface composed of many small grids. The locking fraction at each small grid is iteratively calculated using an inversion algorithm, thereby obtaining a continuous distribution of the locking fraction of the entire fault surface. The number of nodes on the fault surface and the depth are the optimal values obtained through extensive trial calculations. All these settings are aimed at obtaining more accurate crustal deformation characteristics.

Figure 4.
A simplified map and boundary conditions of the faults model. The dip angle of the LHSF and HYF was 60°, and the dip angle of the LPSF and GGBJF was 45°. The fault was completely locked ( = 1) at 0.1 km and free-slipped ( = 0) below 25 km. The LHSF had two surface nodes and six depth nodes, the HYF had eight surface nodes and six depth nodes, the LPSF had six surface nodes and six depth nodes, and the GGBJF had three surface nodes and six depth nodes.

Inversion Principle
The model of interseismic deformation inversion of fault is the block-fault back-slip dislocation model, which is calculated by the TDEFNODE program. Constrained by the GPS velocity field, the model can obtain the block rigid motion parameters of the elastic lithosphere, the uniform strain rate inside the block, and the surface elastic deformation caused by the fault locked at the block boundary by inversion [2,[34][35][36]. The rigid motion of the block is represented by the Euler pole, the uniform strain rate inside the block is represented by the spherical formula [37], and the slip deficit rate on the fault surface is based on the elastic half-space model [38]. The detailed calculation formula is shown in McCaffrey (2002). One pure dynamical scalar is used to represent the degree of locking on the fault surface [34,36]. The locking fraction has no unit, = 0 indicates the fault is completely creeping, = 1 indicates the fault is completely locked, and between 0 and 1 indicates that the fault is partially locked; its mathematical expression is = 1 − / , , which represents the short-term creep rate of the fault, and represents the long-term slip rate of the fault. The fault geometry is represented by an irregular grid formed by three-dimensional nodes (longitude, latitude, depth) on the fault surface. The value of each node is calculated in the inversion. In order to obtain the linear change in value between adjacent nodes, the locking fraction of each divided sub-fault surface is obtained by bilinear interpolation. The slip deficit rate on the fault surface can be expressed by the product of the locking fraction and the long-term slip rate .
The goodness-of-fit of the model to the observed data can be expressed by the chi-square value (sum of chi-square/degree of freedom) [34]: In the formula, is the number of observations, is the number of parameters to be estimated, f is the degree of freedom; is the residual error of the observation, is the mean square error of the observation, and is the weight of the observation. When

Inversion Principle
The model of interseismic deformation inversion of fault is the block-fault backslip dislocation model, which is calculated by the TDEFNODE program. Constrained by the GPS velocity field, the model can obtain the block rigid motion parameters of the elastic lithosphere, the uniform strain rate inside the block, and the surface elastic deformation caused by the fault locked at the block boundary by inversion [2,[34][35][36]. The rigid motion of the block is represented by the Euler pole, the uniform strain rate inside the block is represented by the spherical formula [37], and the slip deficit rate on the fault surface is based on the elastic half-space model [38]. The detailed calculation formula is shown in McCaffrey (2002) [35]. One pure dynamical scalar φ is used to represent the degree of locking on the fault surface [34,36]. The locking fraction φ has no unit, φ = 0 indicates the fault is completely creeping, φ = 1 indicates the fault is completely locked, and φ between 0 and 1 indicates that the fault is partially locked; its mathematical expression is φ = 1 − V c /V, V c , which represents the short-term creep rate of the fault, and V represents the long-term slip rate of the fault. The fault geometry is represented by an irregular grid formed by three-dimensional nodes (longitude, latitude, depth) on the fault surface. The φ value of each node is calculated in the inversion. In order to obtain the linear change in φ value between adjacent nodes, the locking fraction of each divided sub-fault surface is obtained by bilinear interpolation. The slip deficit rate on the fault surface can be expressed by the product of the locking fraction φ and the long-term slip rate V.
The goodness-of-fit of the model to the observed data can be expressed by the chisquare value χ 2 n (sum of chi-square/degree of freedom) [34]: In the formula, n is the number of observations, m is the number of parameters to be estimated, f is the degree of freedom; r i is the residual error of the observation, σ i is the mean square error of the observation, and f is the weight of the observation. When n approaches 1, we considered that the model fits the observed data well. The weighted root-mean-square error (W rms ) is defined as (McCaffrey, 2002) [35]:

Fault Slip Rate
The slip rate of the faults quantitatively describes the relative motion of the hangingwall and footwall blocks of the fault, which is of great significance for extracting the characteristics of the deformation field. The TDEFNODE uses the inversed block motion parameters to calculate the slip rate of the boundary fault [39].
The slip rates of the LHSF, HYF, LPSF, and GGBJF are shown in Tables 1 and 2 and Figure 5. The results of the best-fit model indicate that the LHSF is controlled by a left-lateral strike-slip at a rate of~3.5 mm/yr, which is consistent with geology [11] and geodesy [12,15,17,18,24], approaching the result of 4.1 ± 0.4 mm/yr [21]. The HYF is controlled by a left-lateral strike-slip at a rate of 2.8-3.3 mm/yr, and the average rate is about 3 mm/yr, which is close to the results of 3.3-9.2 mm/yr obtained by geology [7] and 4.5 ± 1.1 mm/yr [11], is close to the GPS, InSAR, and other geodetic methods [12,15,17,18,23,24], and is also consistent with the rates of 3.9 ± 0.4 mm/yr, 3.7 ± 0.4 mm/yr, and 3.6 ± 0.4 mm/yr for the western, middle, and eastern sections of the HYF [21]. The average contraction rate is 0.6 mm/yr, and the maximum rate is 1.4 mm/yr in the middle of the fault. Overall, the slip rate of the HYF obtained by the geological method is 2-6.5 mm/yr [10,40]. Derived from the GPS, InSAR, and other geodetic methods, the slip rate of the HYF obtained by different scholars using the block model and dislocation model is 1.2-9 mm/yr; our results are within the range of the above results. The LPSF is controlled by a thrust contraction at a rate of 1.7-1.9 mm/yr and a slight left-lateral strike-slip at a rate of 0.1-0.9 mm/yr, which is close to the result rates of 1.3 ± 0.1 mm/yr and 1.7 ± 0.1 mm/yr for the north section and south section obtained by [19], at 1.4 mm/yr and 1.2 mm/yr [22]. The GGBJF has both left-lateral strike-slip and thrust extrusion components. The strike-slip rate is 1.3-1.4 mm/yr and the contraction rate is 1.2-1.3 mm/yr, which is consistent with the result of [22]. The slip rate of the HYF appears consistent with the results of [21], while the slip rate of the LPSF differs significantly. Similarly, the slip rate of the LPSF appears consistent with the results of [22], while the HYF differs significantly. We believe that the reasons for these differences are as follows. Firstly, this paper considered the impact of the activities of the near-field faults. Additionally, this paper included GPS observations from the near-field of the faults. In terms of model settings, they only set free nodes and locking fractions decreasing along the depth, while we set the locking fraction as a function of the depth. Finally, when constructing the fault, we considered the actual surface trace and three-dimensional geometry of the active fault. These differences are the reasons for different results.   To analyze the impact of fault node density along the fault strike direction on the inverted slip rates, we further increased the density of the nodes along the actual trace of the faults. The distance between adjacent nodes was set at 3-6 km. The inversion was then carried out again while ensuring that other parameters remained unchanged. The resulting fault slip rates are shown in Figure S3. Although there are differences in the inverted slip rates in some areas compared to the previous model, the Haiyuan Fault and the Liupanshan Fault still maintain their characteristic left-lateral strike-slip and thrust compression, respectively. This difference may be due to the different orientations between the newly added fault nodes and the original ones. Considering the limited spatial distribution density of the GPS stations, we believe that dividing dense fault nodes may not effectively improve the final results. The optimal inversion model should not only consider the surface nodes of the fault but also the maximum locking depth, fault dip angle, parameterization type of nodes, global restrictions of parameters, and smoothing methods. Using the residual GPS velocity field as the criterion, a large number of trials were carried out to obtain the optimal inversion result. Ref. [8] proposed that the total amount of the left-lateral strike-slip component of the HYF should be balanced with the crustal shortening in the LPSF. Based on the inversion results of the model, there is a difference of about~1 mm/yr. We speculate that this part of the differential motion may be converted into strain accumulation in this region. Refs. [41,42] suggested that the north section and the south section of the LPSF both have the characteristics of thrust movement, and the thrust component of the north section is higher than the south section, which is consistent with our results.
To analyze the impact of fault node density along the fault strike direction on the inverted slip rates, we further increased the density of the nodes along the actual trace of the faults. The distance between adjacent nodes was set at 3-6 km. The inversion was then carried out again while ensuring that other parameters remained unchanged. The resulting fault slip rates are shown in Figure S3. Although there are differences in the inverted slip rates in some areas compared to the previous model, the Haiyuan Fault and the Liupanshan Fault still maintain their characteristic left-lateral strike-slip and thrust compression, respectively. This difference may be due to the different orientations between the newly added fault nodes and the original ones. Considering the limited spatial distribution density of the GPS stations, we believe that dividing dense fault nodes may not effectively improve the final results. The optimal inversion model should not only consider the surface nodes of the fault but also the maximum locking depth, fault dip angle, parameterization type of nodes, global restrictions of parameters, and smoothing methods. Using the residual GPS velocity field as the criterion, a large number of trials were carried out to obtain the optimal inversion result.

Fault Locking Slip Distribution
According to fault geometry and optimal model results in Section 3.3, the distribution of interseismic locking degree and slip deficit rate of each fault on the northeastern margin  Table 3). In Figure 6, the locking degree of each fault is not uniformly distributed in the depth. The locking degree of the LHSF is weak, the distribution is uniform and is strongest within 1 km, and the maximum locking fraction is~0.9, decreases to 0.302 at a depth of 3-4 km, and completely creeps below 5 km. The HYF is generally in the creeping stage, but the locking degree is stronger than the LHSF. The locking fraction of the west section gradually increases along the depth direction, is strongest at a depth of 0-5 km, gradually decreases below 5 km, and almost completely creeps below 10 km. The east section exhibits the same behavior as the west section. The locking in the middle section is the strongest and is evenly distributed with a depth of 7-8 km, gradually decreases to a depth of 13 km, and completely creeps below 15 km. The HYF had an 8.5 magnitude earthquake in 1920, and the focal depth is generally believed to be 17-20 km [14]. Our results showed that the fault is still relatively weak at present. We considered that the fault is still in the post-earthquake adjustment stage under the influence of the earthquake, and its movement is still dominated by post-earthquake creep. In the LPSF, the overall locking is the highest, with the south section having stronger locking than the north section. The fault is completely locked within a depth of 5 km and is still strongly locked within 10 km. The locking fraction decreases gradually at a depth between 10 km to 15 km and completely creeps below 15 km. Therefore, it is speculated that the fault may have a large strain accumulation. The locking fraction of the GGBJF decreases gradually from the north section to the south section. The locking fraction is stronger when the depth is above 5 km, gradually decreases between 5-10 km, and completely creeps below 10 km. Compared with previous studies, the locking depth of the LPSF and GGBJF obtained in this paper is slightly lower, which may be related to different models and observations, but the overall trend is consistent. The locking characteristics of each section of the fault zone will be discussed in detail in Section 5.3.

Fault Locking Slip Distribution
According to fault geometry and optimal model results in Section 3.3, the distribu tion of interseismic locking degree and slip deficit rate of each fault on the northeaster margin of the Tibetan Plateau were obtained (Figures 6 and 7, Table 3). In Figure 6, th locking degree of each fault is not uniformly distributed in the depth. The locking de gree of the LHSF is weak, the distribution is uniform and is strongest within 1 km, an the maximum locking fraction is ~0.9, decreases to 0.302 at a depth of 3-4 km, and com pletely creeps below 5 km. The HYF is generally in the creeping stage, but the lockin degree is stronger than the LHSF. The locking fraction of the west section gradually in creases along the depth direction, is strongest at a depth of 0-5 km, gradually decrease below 5 km, and almost completely creeps below 10 km. The east section exhibits th same behavior as the west section. The locking in the middle section is the strongest an is evenly distributed with a depth of 7-8 km, gradually decreases to a depth of 13 km and completely creeps below 15 km. The HYF had an 8.5 magnitude earthquake in 1920 and the focal depth is generally believed to be 17-20 km [14]. Our results showed tha the fault is still relatively weak at present. We considered that the fault is still in th post-earthquake adjustment stage under the influence of the earthquake, and its move ment is still dominated by post-earthquake creep. In the LPSF, the overall locking is th highest, with the south section having stronger locking than the north section. The fau is completely locked within a depth of 5 km and is still strongly locked within 10 km The locking fraction decreases gradually at a depth between 10 km to 15 km and com pletely creeps below 15 km. Therefore, it is speculated that the fault may have a larg strain accumulation. The locking fraction of the GGBJF decreases gradually from th north section to the south section. The locking fraction is stronger when the depth i above 5 km, gradually decreases between 5-10 km, and completely creeps below 10 km Compared with previous studies, the locking depth of the LPSF and GGBJF obtained i this paper is slightly lower, which may be related to different models and observations but the overall trend is consistent. The locking characteristics of each section of the fau zone will be discussed in detail in Section 5.3.  A quantitative analysis of the strain accumulation on the fault surface is necessary to understand the seismic activity of the fault, which can be reflected by the locking degree of the fault; the stronger the locking degree, the easier the strain accumulation. The rate of strain accumulation can be quantitatively described by the slip deficit rate. The slip deficit rate could be calculated by the product of the fault long-term slip rate V and the locking fraction φ [35]. Figure 7 and Table 3 show the slip deficit rate distribution of the faults.  A quantitative analysis of the strain accumulation on the fault surface is necessary to understand the seismic activity of the fault, which can be reflected by the locking degree of the fault; the stronger the locking degree, the easier the strain accumulation. The rate of strain accumulation can be quantitatively described by the slip deficit rate. The slip deficit rate could be calculated by the product of the fault long-term slip rate V and the locking fraction [35]. Figure 7 and Table 3 show the slip deficit rate distribution of the faults.
For the same fault, the area where the slip deficit rate is significantly quick is the strongly locked region of the fault. The rate of slip deficit gradually decreases with depth along the fault surface until it reaches 0 mm/yr. This indicates that the fault surface no longer impedes the relative motion of the two blocks, and the fault is completely creeping. [34] pointed out that the 'slip deficit rate' means that the energy not released by the long-term slip rate of the fault through the creep slip is converted into strain energy to accumulate on the fault surface, so we used the slip deficit rate to quantitatively analyze the speed of strain accumulation on the fault surface. Similar to the distribution of the locking, the slip deficit rate of each fault is also unevenly distributed in depth, but the magnitude and trend of the slip deficit rate are similar to the locking distribution. The slip deficit rate decreased from 3.6 mm/yr at the LHSF to 1.8 mm/yr at the GGBJF. Figure 7 shows that the LHSF and HYF are mainly controlled by a left-lateral slip deficit rate, among which the LHSF has the largest slip deficit rate at a rate of 3.6 mm/yr, but its depth is the shallowest, distributed from the surface to 5 km. The strike-slip deficit rate of the HYF is large, with a maximum rate of 3.2 mm/yr, and the depth in the middle section is the deepest, with the west section showing an increasing slip deficit rate with depth, while the east section shows a decreasing slip deficit rate toward the surface, with  For the same fault, the area where the slip deficit rate is significantly quick is the strongly locked region of the fault. The rate of slip deficit gradually decreases with depth along the fault surface until it reaches 0 mm/yr. This indicates that the fault surface no longer impedes the relative motion of the two blocks, and the fault is completely creeping. [34] pointed out that the 'slip deficit rate' means that the energy not released by the long-term slip rate of the fault through the creep slip is converted into strain energy to accumulate on the fault surface, so we used the slip deficit rate to quantitatively analyze the speed of strain accumulation on the fault surface. Similar to the distribution of the locking, the slip deficit rate of each fault is also unevenly distributed in depth, but the magnitude and trend of the slip deficit rate are similar to the locking distribution. The slip deficit rate decreased from 3.6 mm/yr at the LHSF to 1.8 mm/yr at the GGBJF. Figure 7 shows that the LHSF and HYF are mainly controlled by a left-lateral slip deficit rate, among which the LHSF has the largest slip deficit rate at a rate of 3.6 mm/yr, but its depth is the shallowest, distributed from the surface to 5 km. The strike-slip deficit rate of the HYF is large, with a maximum rate of 3.2 mm/yr, and the depth in the middle section is the deepest, with the west section showing an increasing slip deficit rate with depth, while the east section shows a decreasing slip deficit rate toward the surface, with the distribution in a depth range from the surface to 15 km. The LPSF is mainly controlled by the thrust dip-slip deficit rate, with a maximum rate of 1.9 mm/yr, but its distribution depth of the slip deficit is relatively uniform, ranging from the surface to a depth of 15 km. The GGBJF has both left-lateral strike-slip deficit rate and thrust dip-slip deficit rate components. The rate gradually decreases from the north section to the south section at a maximum rate of 1.8 mm/yr and is distributed within a depth range of 0-10 km.

Block Motion
The optimal parameter model showed that the present model can better explain the surface GPS observations and obtain the rigid rotation parameters and internal strain parameters of the Alxa blocks and Lanzhou blocks. Figure 8 shows that the GPS velocity residuals of this model are randomly distributed and consistent with a Gaussian distribution. Most of the residuals are within 1 mm and have a random direction. Some individual points have residuals between 2~3 mm, but they are far from the fault zone in this model; this may be due to the existence of active faults with strong activity in the area of larger residuals, such as the West Qinling Fault, which has no significant effect on the actual inversion results of this model. Among 156 residual distributions, 12 residuals exceeded 1 mm, 2 residuals exceeded 2 mm, and 92% of the residuals are within 1 mm. This indicates that the fitting degree between this model and GPS observation data is good, and no obvious signals were missed. trolled by the thrust dip-slip deficit rate, with a maximum rate of 1.9 mm/yr, but its distribution depth of the slip deficit is relatively uniform, ranging from the surface to a depth of 15 km. The GGBJF has both left-lateral strike-slip deficit rate and thrust dip-slip deficit rate components. The rate gradually decreases from the north section to the south section at a maximum rate of 1.8 mm/yr and is distributed within a depth range of 0-10 km.

Block Motion
The optimal parameter model showed that the present model can better explain the surface GPS observations and obtain the rigid rotation parameters and internal strain parameters of the Alxa blocks and Lanzhou blocks. Figure 8 shows that the GPS velocity residuals of this model are randomly distributed and consistent with a Gaussian distribution. Most of the residuals are within 1 mm and have a random direction. Some individual points have residuals between 2~3 mm, but they are far from the fault zone in this model; this may be due to the existence of active faults with strong activity in the area of larger residuals, such as the West Qinling Fault, which has no significant effect on the actual inversion results of this model. Among 156 residual distributions, 12 residuals exceeded 1 mm, 2 residuals exceeded 2 mm, and 92% of the residuals are within 1 mm. This indicates that the fitting degree between this model and GPS observation data is good, and no obvious signals were missed. Since the ORDOS block is a rigid block, the model set it as the reference block without calculating its internal strain, so the motion of each block obtained is relative to the motion of the ORDOS block. The obtained rigid rotation parameters and internal strain parameters of each block are shown in Tables 4 and 5. Relative to the ORDOS block, both the Alxa block and the Lanzhou block make clockwise rotations, and the former has a larger rotation rate than the latter. The Euler pole of the Alxa block is far from the block itself and is located on the northeastern of the ORDOS block; the Euler pole of the Lanzhou block is far from the inner part of the block and is located on the Since the ORDOS block is a rigid block, the model set it as the reference block without calculating its internal strain, so the motion of each block obtained is relative to the motion of the ORDOS block. The obtained rigid rotation parameters and internal strain parameters of each block are shown in Tables 4 and 5. Relative to the ORDOS block, both the Alxa block and the Lanzhou block make clockwise rotations, and the former has a larger rotation rate than the latter. The Euler pole of the Alxa block is far from the block itself and is located on the northeastern of the ORDOS block; the Euler pole of the Lanzhou block is far from the inner part of the block and is located on the south of the ORDOS block. The motion trend of the block is consistent with the collision of the Indian plate in the N20 • E direction [25] with the Eurasian plate, resulting in an eastward extrusion of the northeastern margin of the Tibetan Plateau. However, when the relatively stable ORDOS block is encountered, the movement direction is deflected to the south, showing clockwise movement along the fault zone. Compared with the relatively stable Alxa block, the internal strain rate of the Lanzhou block is larger and dominated by compressive strain and supplemented by tensile strain, which is due to the lateral escape of the crust blocked by the relatively stable Alxa block and the hard ORDOS block [43]; the principal compressive strain rate . ε 1 = 14.34 nanostrain/yr and the tensile strain rate . ε 2 = 6.05 nanostrain/yr. Compared with the results of previous studies, the magnitude of the kinematic parameters of the blocks differs due to the different block divisions, but the kinematic trend of the block is consistent [41,44,45]. Table 4. Euler parameters and the internal strain rate of each block relative to the ORDOS block. The angular velocity of the rigid body rotation is positive with a counterclockwise rotation.

Eulerian Vector
Internal Strain Rate (Nanostrain/yr) In order to verify whether adding the LHSF and GGBJF to the model improved the data fitting degree, four sets of experiments were done for comparison. In the first group, the HYF and LPSF were modeled; in the second group, the LHSF, HYF, and LPSF were modeled; in the third group, the HYF, LPSF, and GGBJF were modeled; and in the fourth group, the LHSF, HYF, LPSF, and GGBJF were modeled and analyzed. When modeling the faults, only the motion of the surrounding circumscribed blocks was considered, and only the data within the corresponding blocks were used while keeping the other model parameters unchanged. The modeling results are shown in Table 6. The results show that the addition of the LHSF and GGBJF can improve the overall fitting effect of the GPS observations, reducing both the χ 2 n and W rms of the data within the Alxa and Lanzhou blocks. The effect of adding the GGBJF in the experiment is slightly greater than the LHSF, which may be due to the existence of shallow creeps near the LHSF and the greater activity of the GGBJF, as well as the more abundant near-field observations. When modeling the fault zone, both the LHSF and GGBJF were added in the inversion finally.

Optimal Model Parameter Settings
Parameter settings in the model mainly include selecting reference blocks, block properties (rigid or elastic), fault parameterization type, the limit of locking depth, constructing faults (node, fault slip type, locking depth, fault dip angle, node index), making pseudofaults (closure block), etc. The locking depth is the depth at which the fault slips completely freely, and the locking depth was set to 25 km in this model. Previous research showed that the ORDOS block is stable and rigid, so the ORDOS block was set as the reference block in the inversion, and its internal strain was not calculated. The Lanzhou block and Alxa block were set as elastic blocks, and their internal strains were calculated. The geometrical characteristics of the faults had been studied extensively, and some results on the dip angle of the major faults are given in Table 7. The dip angles of the LHSF and HYF are between 41 • and 70 • , and the optimal dip angle was 60 • by a grid search; the dip angle of the LPSF is between 20 • and 70 • , and the optimal dip angle was 45 • by a grid search; and the dip angle of the GGBJF is between 45 • and 88 • , and the optimal dip angle was 45 • by a grid search. The initial locking depths of the faults can also be calculated by a grid search. For the LHSF and HYF, the initial locking depth is 15 km by searching in the range of 10-20 km; for the LPSF and GGBJF, the initial locking depth is 15 km by searching in the range of 15-25 km. With other parameters unchanged, the model performed inversion at intervals of 5 km at locking depths between 15 km and 40 km. Comparing the χ 2 n of the results, it was concluded that 25 km was the optimal choice for this model. In addition, according to the catalog of earthquakes above M3-magnitude in mainland China and its surrounding areas [50] and the study on the precise location of earthquakes [51], the focal depth of most earthquakes on the northeastern margin of the Tibetan Plateau is less than 25 km. Therefore, the final locking depth was set as 25 km.
The fault node settings were given in Section 3.3. The locking fraction between 0.1 km and 25 km was inverted by the block-fault back-slip dislocation model (0 ≤ φ ≤ 1), and the fault locking fraction was set based on [33], who set the node φ values as a function of depth, which is expressed as follows.
In Formula (3), G is the shape parameter that constrains the decay of the locking fraction, Z1 is the upper boundary of the transition region, the depth where the fault is completely locked, and Z2 is the lower boundary of the transition region, the depth where the fault slips freely. When z ≤ Z1, the fault is completely locked; when z ≥ Z2, the fault is completely creep-slipped; and in Z1 < z < Z2, the locking fraction of the fault node decays as a function of Formula (3).

Distribution Characteristics of the Degree of Locking and Slip Deficit Rate in Different Segments of the Fault Zone
The results show that the LHSF and HYF are controlled by a left-lateral strike-slip, the LPSF is controlled by a thrust dip-slip, and the GGBJF has both left-lateral strikeslip and thrust dip-slip components. As shown in Figures 6 and 7, the locking in the LHSF is uniformly distributed, with a strong locking within 1 km, followed by a gradual decrease in the locking fraction to a complete creep slip below 5 km with a slip deficit rate of 3.4 mm/yr, which is consistent with the results of [22,40,42]. The 1920 Haiyuan 8.5 magnitude earthquake in the HYF formed a 200 km long left-lateral strike-slip rupture zone at the surface [1]. The overall HYF is in the creep-slip stage, with an overall locking depth of 15 km, but the locking fraction between the depth of 13 km and 15 km is quite small, so the locking depth is considered to be 13 km. The locking of the western section is weak, with the locking fraction gradually increasing along the strike, and is strong within 5 km, gradually decaying to a depth of 10 km, and the slip deficit rate is 2.8-3.2 mm/yr. The middle section has the strongest locking and relatively uniform distribution with depths reaching 7-8 km, gradually decaying to 13 km and completely creeping below 13 km, and its slip deficit rate is 2.4-2.8 mm/yr; the east section has the same expression as the western, except that its locking fraction gradually decays in the depth along the strike, and the slip deficit rate is~2.4 mm/yr. In summary, the locking of the HYF is similar to the arc-shaped distribution (Figure 6), with a characteristic of weak locking in the west and east sections and strong locking in the middle section, and the slip deficit rate decreases from 3.2 mm/yr to 2.4 mm/yr from west to east. In general, the strain accumulation capacity of the HYF is relatively low. Combined, the earthquake recurrence interval is 800-1600 years [1], and the focal depth is generally considered to be 17-20 km, so we considered that the HYF is still in the post-earthquake stress adjustment stage, which better explains the absence of earthquakes of a 6.0 magnitude or higher in this fault since 1920. The locking of the LPSF is, overall, the strongest and most evenly distributed, with slightly higher locking in the southern section compared to the northern. The high-value locking fraction areas are distributed almost across the entire fault surface within 10 km (Figure 6), the fault is completely locked within 5 km depth (φ = 1), a strong locking fraction (φ ≥ 0.5) still exists within 5-10 km and gradually decreases between 10 km and 15 km and completely creeps below 15 km, and the overall slip deficit rate is 1.7-1.9 mm/yr, reflecting the characteristic that the south section is higher than the north section. The overall locking of the GGBJF is weak, and the locking fraction gradually decreases along the strike at the depth, with the strongest locking near the south section of the LPSF where the locking degree is strong at 5 km, the locking fraction gradually decreases between 5 km and 10 km in depth and completely creeps below 10 km, and its slip deficit rate decreases from 1.8 mm/yr to 0.2 mm/yr from the north section to the south section. The north-south segmented movement characteristics of the LPSF are clearly evident [52]. The south section of the LPSF refers to the GGBJF in this paper; in this way, our results show that the north-south section of the LPSF has obvious motion characteristics, and the north section is more active than the south section. In terms of the locking degree, the LPSF is easier to strain accumulation, but in terms of the slip deficit rate, it accumulates strain slowly; however, no earthquake with a magnitude of 7.0 or higher had occurred in the LPSF in the past 1400 years [53], so it will still reach high strain accumulation after long-term accumulation, which will lead to a strong earthquake, such as the Wenchuan earthquake.
The results also found that there was a weak right-lateral strike-slip phenomenon with a rate of 0.1-0.2 mm/yr in the top area of the HYF when it was converted from the sinistral strike-slip region at the tail end of the east section to the north section of the LPSF. If we consider the explanation of this phenomenon, we believe that this area is located on the junction transition area of three blocks, while the Alxa block and Lanzhou block are rotating clockwise around the ORDOS block; the Lanzhou block is also slightly rotating clockwise around Alxa block. However, due to the obstruction of the hard ORDOS rigid block to the east, the motion of the block boundary in the transitional area changes, resulting in a right-lateral strike-slip. There is another explanation: the displacement along a strike-slip fault must be absorbed and adjusted at a certain part along the fault, often at both ends of the fault, so that the whole fault can reach a new balance in energy; the contraction area at the end of the HYF is this kind of tectonic deformation [26]. Therefore, we can also assume that the energy accumulated by the left-lateral strike-slip at the end of the HYF may be released by the right-lateral strike-slip at the top of the north section of the LPSF to achieve a relative energy balance. Secondly, the phenomenon of weak left-lateral extension in some sections of the HYF in our results can also be explained by this theory.

Conclusions
In this paper, we collected the observations of encrypted GPS stations around the Haiyuan-Liupanshan fault zone, integrated the GPS velocity field of the northeastern margin of the Tibetan Plateau and neighboring areas from 2010 to 2020, and obtained the present-day GPS velocity field of the northeastern margin of the Tibetan Plateau. We divided the region around the fault zone into independent blocks, taking into account the fine fault distribution of the Haiyuan-Liupanshan fault zone to obtain the optimal inversion model. A block-fault back-slip dislocation model was used to invert the slip rate, locking depth, and the slip deficit rate of the fault. Combining the results of historical seismic data, regional geological data, and the motion characteristics of each fault, we discussed and analyzed the locking degree and slip deficit rate distribution characteristics of different sections of the fault zone. The main conclusions are as follows: (1) The LHSF is mainly controlled by a left-lateral strike-slip at a rate of~3.5 mm/yr.
The fault is completely locked within a 1 km depth and completely creeps below 5 km, so the locking degree of this fault is low. Combined with the slip deficit rate of 3.2-3.6 mm/yr, we considered that the seismic risk of this fault is low. (2) The HYF is mainly controlled by a left-lateral strike-slip at a rate of 3.0-3.2 mm/yr.
The western and eastern sections of the fault are weakly locked, with strong locking at a depth of 5 km, and completely creeps at a depth of 10 km or less. The whole fault completely creeps below 13 km, and the overall locking degree showed the distribution characteristics of strong in the middle section and weak in the eastern and western sections. Combined with the slip deficit rate of 2.4-3.2 mm/yr, we believe that there is a certain amount of strain accumulation in the middle section, but combined with the regional geological and historical seismic data, we believe that this fault is still in the post-earthquake stress adjustment stage, the seismic risk is considered low. (3) The LPSF is mainly controlled by a thrust dip-slip at a rate of 1.7-1.9 mm/yr. The overall degree of locking in this fault is the strongest and uniformly distributed, it is completely locked at a depth of 5 km, there is still strong locking at a depth of 10 km, and the fault completely creeps at a depth of 15 km or less; the slip deficit rate of this fault is 1.7-1.9 mm/yr. Although its rate of strain accumulation is slow, when combining the analysis of geology, historical seismic data, and the activity characteristics, there is still a risk of a moderate to strong earthquake after a long period of strain accumulation, which should be continuously observed and analyzed. (4) The GGBJF has both left-lateral strike-slip and thrust dip-slip components, its strikeslip rate is 1.3-1.4 mm/yr and its dip-slip rate is 1.2-1.3 mm/yr, the locking fraction decreases gradually from north to south, locking is strong in the range of 5 km in depth, it completely creeps below 10 km, and is combined with a slip deficit rate of 1.8 mm/yr and decays to 0.2 mm/yr from west to east. The seismic risk of this fault is considered to be low.