*Article* **The Current Crustal Vertical Deformation Features of the Sichuan–Yunnan Region Constrained by Fusing the Leveling Data with the GNSS Data**

**Yong Zhang 1,2, Caijun Xu 1,\*, Zhijiang Zheng 3, Hongbao Liang <sup>3</sup> and Shuang Zhu <sup>3</sup>**


**Abstract:** This study uses the least squares collocation method to fuse the leveling vertical deformation velocity in the Sichuan–Yunnan region with the GNSS observations of this region from 320 stations in the China Crustal Movement Observation Network (CMONOC) and the China Continental Tectonic Environment Monitoring Network (CMTEMN) from 1999 to 2017. Such fusion is to improve the accuracy of the vertical deformation rates in large spatial scales. The fused vertical deformation results show that: (1) the fused deformation field has a uniform spatial distribution, and shows detailed change characteristics of key regions; (2) the current vertical crustal motion in this region is featured by the contemporaneous occurrence of crustal compression, shortening and uplift and basin extensional subsidence; (3) most areas in this region experience uplifts, as the lateral push of the Qinghai–Tibet Plateau was blocked by the Sichuan Basin. The areas on the northwest side of the Longmenshan fault and the Lijiang-Xiaojinhe fault are dominated by uplifts, with the velocity of 1.5 mm/a–5.5 mm/a, and the region on the southeast side has slight uplifts, with the velocity of 1.0 mm/a–1.5 mm/a; (4) many areas have high gradient vertical deformation, especially the region close to the Wenshan fault and on the two sides of the Yarlung Zangbo fault that has the value of 3.0–4.0 <sup>×</sup> <sup>10</sup><sup>−</sup>8/a, deserving further attention to be paid to the long-term earthquake hazards.

**Keywords:** GNSS; leveling; Sichuan–Yunnan region; vertical deformation

#### **1. Introduction**

The Sichuan–Yunnan region (SYR) is located on the southeastern margin of the Qinghai–Tibet Plateau (QTP) in southwestern China, and is the edge of the area bearing the interaction between the Eurasian plate and the Indian plate. The northward subduction of the Indian plate underneath the Eurasian continent resulted in the northeastward extrusion of the QTP and pushed the plateau materials flow to the east, which, however, was blocked by the hard South China block. As the main channel of the lateral extrusion in the middle QTP, the SYR experienced particularly strong crustal motions, so this region has great topographic fluctuation and frequent seismic activity. There are many large-scale faults with great slip rates [1–7]. The material extrusion from the middle QTP together with the soft lower crust cause the strong compression of the Longmenshan area in the southeast of the Bayankala block, forming an active compression fault system consisting of the Longmenshan fault and the Minjiang fault, etc. Inside the SYR, the faults are characterized by the extrusion caused by strike-slip. From these faults, a high-speed left-lateral strike-slip fault system has developed, consisting of the Ganzi–Yushu fault, the Xianshuihe fault, the Anning River fault, and the Zemuhe fault, Daliangshan fault, Xiaojiang fault, Litang fault, Lijiang-Xiaojinhe fault, and East Kunlun fault (eastern segment). A corresponding rightlateral strike-slip fault system has also developed, consisting of the Zhongdian fault, the

**Citation:** Zhang, Y.; Xu, C.; Zheng, Z.; Liang, H.; Zhu, S. The Current Crustal Vertical Deformation Features of the Sichuan–Yunnan Region Constrained by Fusing the Leveling Data with the GNSS Data. *Remote Sens.* **2022**, *14*, 1139. https://doi.org/ 10.3390/rs14051139

Academic Editor: Alex Hay-Man Ng

Received: 23 January 2022 Accepted: 23 February 2022 Published: 25 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Red River fault, the Qujiang fault zone, the Nujiang fault zone, and the Jinsha River fault. The left-lateral strike-slip of the SYR together with the right-lateral strike-slip of the other blocks contribute to the conjugate strike-slip of the southern and western Yunnan areas. The faults experiencing left-lateral strike-slip are the Wanding fault, the Nanting River fault, the Daluo fault, the Jinghong fault, the Dien Bien Phu fault, and that experiencing the rightlateral strike-slip is the Longling-Lancang fault. Since the last century, 122 earthquakes of magnitude 6.0–6.9 and 23 earthquakes of magnitude 7.0–7.9 have occurred in the region, and the largest recorded earthquakes in the SYR are the 1833 Songming M 8.0 earthquake and the 2008 Wenchuan Mw8.0 earthquake.

Accommodating the strongest crustal deformation and seismic activity in the QTP, the SYR has been a focus for studying the present-day crustal deformation model and tectonic evolution of the QTP. GPS observations show a large-scale southeastward lateral movement on the southeastern margin of the QTP [8,9]. The "rigid block extrusion mode" [10,11] and the "lower crustal flow mode" [12–15] were developed to constrain the internal deformation and lateral extrusion of the QTP, respectively. The "lower crustal flow model" believes that the middle and lower crusts of the QTP are weak, where the crust thickening and shortening on the eastern margin of the plateau are most likely to occur [14,16,17]. This model has been widely accepted as it can well explain the crustal deformation characteristics observed on the surface. Furthermore, seismologic observations [18–24] and magnetotelluric studies [25–29] have reported low-velocity and high-conductivity/lowresistivity layers in the middle and lower crusts, which confirm the existence of the middle and lower crust flow.

Although the "lower crustal flow model" can explain the deformation characteristics of the eastern margin of the QTP well, whether the lower crustal flow in the eastern QTP has entered the SYR and how it distributes are still debated in the geosciences field. Some scholars have suggested that the lower crust of the eastern margin of the QTP is unlikely to have large-scale ductile flows due to the constraints of faults and tectonic boundaries. There may be flows along two arc-shaped channels, which may have passed through the Jinsha River–Red River fault to southern Yunnan [29]. Some scholars believed that the lower crustal flow from the eastern QTP was blocked by the Lijiang-Xiaojinhe fault and stopped at the northern Yunnan [30]. Studying the crustal movements and the deep driving mechanism in the SYR helps to understand the dynamic process of the southeastern boundary of the QTP, and also helps to estimate the future seismic risks of the key structure in China.

The vertical crustal deformation is important to the study of regional subsidence, crustal evolution and fault motion. The research on seismic activity trends based on vertical deformation provides essential information for estimating the regional strong earthquake trends and detecting dangerous areas. The SYR is a seismically active zone, where several strong earthquakes occurred in the past two decades, including the 2008 Wenchuan Mw 8.0 earthquake, 2013 Lushan Mw 7.0 earthquake and 2017 Jiuzhaigou Mw 7.0 earthquake. These earthquakes brought great damages. Therefore, it is of great significance to analyze the vertical movement characteristics of the major structural faults in the region, and to further explore their correlation with strong earthquakes. The National Leveling Network, Earthquake Leveling Monitoring Network, China Crustal Movement Observation Network (CMONOC) and China Continental Tectonic Environment Monitoring Network (CMTEMN) form a relatively complete deformation observation network for the SYR. The rich archive leveling and GNSS observations can be used for studying the crustal vertical deformation characteristics.

The traditional precision leveling and GNSS were separately used to study the vertical crustal deformation [31–33], but their advantages cannot be fully explored. As shown in Figure 1, the leveling route in SYR was generally arranged along the main structural belt, resulting in some observation gaps in the middle part of the secondary blocks. The GNSS stations were uniformly distributed, covering the observation gaps in the leveling route. The two observations are spatially complementary, so combining them can solve the problem of lacking velocity reference when using only leveling observations, and slow down the error transferring and accumulation. It also can constrain the deformed surface, so as to obtain a regional vertical motion field with higher spatial resolution and reasonability. Here, we use the least squares collocation method to obtain the vertical subsidence rate of the SYR, combining high-precision leveling data [31,34] with the observation of 320 GNSS stations in the CMONOC and the CMTEMN from 1999 to 2017. On this basis, we analyze the vertical motion characteristics and dynamic mechanism of this region.

**Figure 1.** The tectonic setting of the Sichuan–Yunnan region (SYR). Triangles and squares denote GNSS stations and leveling points, respectively. Red dots denote the M ≥ 7.0 earthquakes occurred in this region since 1900. Profiles AA , BB and CC will be discussed in Section 3.1. Inset shows the location of the study area in China.

#### **2. Data and Processing**

#### *2.1. Leveling Data*

The leveling data are from Hao et al. [31], containing the leveling vertical velocities of 3386 measuring points from 1977 to 2011. In addition, the velocities of 98 measuring points from 1990 to 2015 provided by Su et al. [34] are used as supplement. Therefore, the leveling data of 3484 measuring points in total during 1970s–2010s are used for this study. All leveling routes have the data accumulation of more than 2 periods to ensure the reliability of the results. A leveling network of several thousand kilometers usually takes many years to complete the first observation. To solve the time inconsistency, we use the linear dynamic adjustment method to obtain the average vertical deformation rates in a long period. For the rank deficiency of the normal equation, we adopt the rate of a small number of GNSS observation stations as the a priori constraint. The GNSS vertical velocity values (vertical Ý

Ý

Ý

Ý

Ý

Ý

motion velocity under ITRF2005) are from Wang [26]. Hao et al. [31] provided the data of eight GNSS sites, including CHDU, XIAG, KUNM, XNIN, JB08, JB09, JB27, JB35 and JB40, and Su et al. [34] provide the data of 4 GNSS sites, which are XIAG, KUNM, JB40 and SCJL. Theoretically, after the constrained adjustment, the benchmarks for the two datasets are the same. On the basis of Hao's results, this study uses the flux isostasy theory [35–37] to obtain the leveling results under the unified benchmark of the two datasets. The velocity uncertainty is estimated by the dynamic leveling adjustment. On this basis, we remove the gross errors and obtain the vertical motion vectors from the remaining 3278 measuring points (Figure 2a).

Ý

Ý

Ý

Ý

Ý

Ý

**Figure 2.** The vertical motion vectors of the SYR obtained (**a**) from the leveling data, (**b**) from the GNSS data, and (**c**) by fusing the leveling and GNSS results.

#### *2.2. GNSS Data*

This study collects the observations of 320 GNSS stations in the CMONOC and CMTEMN for the SYR from 1999 to 2017 and adopts unified models to process the global International GNSS Service (IGS) station data and the GNSS data by using the GAMIT /GLOBK10.61 software [38,39]. The steps are as follows: <sup>1</sup> simultaneously solve the satellite orbits, earth orientation parameters, station coordinates, tropospheric delay and horizontal gradient parameters. Use loosely constrained algorithms for the station positions (5 cm for IGS core stations and 10 cm for other stations); <sup>2</sup> fix the ambiguity; <sup>3</sup> set the satellite masking angle as 10◦; <sup>4</sup> set the sampling interval as 30 s, and reweight the observations according to the post-test phase residuals [40]; <sup>5</sup> calculate the earth tide correction, ocean tide correction and pole tide correction. The ocean tide load model adopts FES2004 [41]; <sup>6</sup> adopt GMF as the tropospheric mapping function [42]. The atmospheric tide load correction and non-tidal atmospheric load correction are not considered; <sup>7</sup> adopt the absolute antenna phase center correction model; <sup>8</sup> use GLOBK for benchmark transformation to obtain the station coordinate time series with respect to ITRF2014 [43]. During the transformation, 7 parameters, including the benchmark transformation estimates, translation, rotation and scale factor, are used to reduce the datum distortion caused by the unmodeled non-tidal atmospheric pressure load; <sup>9</sup> select 77 IGS stations distributed evenly around the world as the frame sites to achieve the self-consistency of the results from the multiphase observations. Finally, the epochs of the stations with the error bigger than 5 mm north, 5 mm east and 20 mm up are removed to obtain the precise coordinate sequence of the GNSS station.

The coseismic deformation and post-seismic deformation of large earthquakes are useful for understanding the long-term tectonic deformation. From 94 continuous GNSS stations in the SYR, we obtain the tectonic velocity by the following principles: for the observations from the stations in the near field of some large earthquakes, such as the 2008 Wenchuan M8.0 earthquake, only the pre-earthquake data are used to avoid the influence of post-earthquake deformation; a small number of stations with no pre-earthquake observations or only with the first-phase observation data are discarded; the data from the stations in the far field are corrected by the estimated coseismic displacement. We analyze the GNSS deformation time series to obtain reliable vertical velocity, during which only the coseismic deformation and unexplained step effects are considered. The vertical velocity is obtained by Formula (1).

$$y(t\_i) = a + bt\_i + c\sin(2\pi t\_i) + d\cos(2\pi t\_i) + e\sin(4\pi t\_i) + f\cos(4\pi t\_i) + \sum\_{j=1}^{n\_\mathcal{S}} g\_j h(t\_i - t\_{\bar{\mathcal{S}}j}) \tag{1}$$

where *a* is the initial position, *b* is the rate of linear motion, *c* and *d* are the annual term coefficients, *e* and *f* are the semi-annual term coefficients, *g* is the magnitude of *j*th step caused by the replacement of equipment or the coseismic displacement and *h* is the Heaviside step function.

In order to ensure the reliability of the vertical velocity from the GNSS station, we use time series analysis to evaluate the observations in terms of observation time span, number of observation periods and station stability, etc. The selected GNSS stations should satisfy the following requirements: the weighted root mean square (WRMS) of the coordinate is less than 3.3 mm horizontally and 10 mm vertically; the observation time is longer than 7 years; the station has at least 4 observation periods. Finally, 226 GNSS stations with low observation noise and little interference from annual motion are selected. The linear motion rate is estimated without considering the annual motion. Figures 3 and 4 shows the time series and fitting results of the coordinates of some campaign and continuous stations, respectively. The residuals after fitting are small, with the horizontal WRMS smaller than 3 mm and the vertical WRMS smaller than 8.0 mm. These results demonstrate that the above processing scheme can obtain reliable horizontal and vertical movement trends of the stations.

In addition, in order to increase the spatial resolution of GNSS stations, we adopt the vertical velocity of 441 GNSS stations from 1998 to 2013 from Liang et al. [44]. Following the regional flux isostasy theory, we superimpose the calculated data on Liang's velocity field to obtain the GNSS vertical velocity in the SYR from 1990s to 2010s. Then, the gross error is removed to obtain the GNSS vertical motion vector (Figure 2b).

**Figure 3.** Time series and errors of the north (**left**), east (**middle**) and vertical (**right**) components of three GNSS campaign sites. The first two rows are for station G352, the third and fourth rows are for station GB01, and the fifth and sixth rows are for station H337. The red lines indicate the estimated velocities. The WRMS of every component are shown at the left up corner of error panels.

**Figure 4.** Time series and errors of three GNSS continuous sites: SCML (first two rows), YNDC (third and fourth rows), SCTQ (fifth and sixth rows). Others are the same as Figure 3.

#### *2.3. Fusion of the Leveling Results and the GNSS Results*

The 1970s–2010s leveling vertical velocity obtained by the linear dynamic adjustment method and the 1990s–2010s GNSS vertical velocity obtained using the time series analysis model are the annual average velocity of the long-time vertical motion and can reflect the background characteristics of the current vertical crustal movement in the SYR. In order to increase the spatial resolution of the deformation velocity by exploring the advantage of the above two results, we develop the least squares collocation method from the method of Vestol [45], which is based on multiquadric function [46–52].

The error equations for the two types of data are as follows

$$\begin{aligned} V\_i &= \sum\_{\substack{i=1\\m}}^m \mathbb{C}\_i Q\_i(\mathbf{x}, y) + \mathbb{S}\_i - l\_i \\ V\_j &= \sum\_{j=1}^m \mathbb{C}\_j Q\_j(\mathbf{x}, y) + \mathbb{S}\_j - l\_j \end{aligned} \Leftrightarrow \quad \begin{aligned} V\_1 &= A\_1 X + B\_1 S - L\_1 \\ V\_2 &= A\_2 X + B\_2 S - L\_2 \end{aligned} \int \Rightarrow V = AX + BS - L \quad \text{(2)}$$

where *Qi*(*x*, *y*) = ((*x* − *xi*) <sup>2</sup> <sup>+</sup> (*<sup>y</sup>* <sup>−</sup> *yi*) 2 ) 1.5 <sup>2</sup> + 1. *Vi* is the error equation for the leveling results. It is obtained by Helmert iterative adjustment, and the corresponding weight matrix is *P*1; *Vj* is the error equation for the GNSS results, and the weight matrix is *P*2. *l* is the observation, *C* is the kernel function coefficient, *Q*(*x*, *y*) is the kernel function, *S* is the random signal, *A* is the design matrix, *X* is the unknown parameter, *B* is the coefficient matrix and *L* is the observation Vector, *V* = *AX* + *BS* − *L* is the overall error equation, with the weight matrix *p*.

From Formula (2), we obtain

$$\begin{cases} \begin{array}{l} \mathcal{X} = \left(A^T \mathcal{C}\_{XX}^{-1} A\right)^{-1} A^T \mathcal{C}\_{XX}^{-1} L \\ \mathcal{S} = \mathcal{C}\_{SS} B^T \mathcal{C}\_{XX}^{-1} \left(I - A\hat{X}\right) \end{array} \tag{3} \end{cases}$$

where *CXX* = *C*ΔΔ + *BCSSB<sup>T</sup>*, with *C*ΔΔ and *CSS* denoting the covariance matrix of the observed value and of the random signal, respectively.

The mean square error of the unit weight is

$$
\begin{bmatrix}
\boldsymbol{\delta}\_{\Lambda}^{2(i)} \\
\boldsymbol{\delta}\_{\boldsymbol{S}}^{2(i)}
\end{bmatrix} = \begin{bmatrix}
n - 2tr\left(\boldsymbol{N}^{-1}\boldsymbol{N}\_{\Lambda}\right) + tr\left(\boldsymbol{N}^{-1}\boldsymbol{N}\_{\Lambda}\right)^{2} & tr\left(\boldsymbol{N}^{-1}\boldsymbol{N}\_{\Lambda}\boldsymbol{N}^{-1}\boldsymbol{N}\_{\boldsymbol{S}}\right) \\
tr\left(\boldsymbol{N}^{-1}\boldsymbol{N}\_{\Lambda}\boldsymbol{N}^{-1}\boldsymbol{N}\_{\boldsymbol{S}}\right) & n - 2tr\left(\boldsymbol{N}^{-1}\boldsymbol{N}\_{\S}\right) + tr\left(\boldsymbol{N}^{-1}\boldsymbol{N}\_{\S}\right)^{2}
\end{bmatrix} \begin{bmatrix}
\boldsymbol{V}^{T}\boldsymbol{P}\_{\Lambda}\boldsymbol{V} \\
\boldsymbol{\tilde{S}}^{T}\boldsymbol{P}\_{\tilde{\Sigma}}\boldsymbol{\tilde{S}}
\end{bmatrix} \tag{4}
$$

where

$$\boldsymbol{N} = \begin{bmatrix} \boldsymbol{A}^{T}\boldsymbol{P}\_{\boldsymbol{\Delta}}^{(i)}\boldsymbol{A} & \boldsymbol{A}^{T}\boldsymbol{P}\_{\boldsymbol{\Delta}}^{(i)}\boldsymbol{B} \\ \boldsymbol{B}^{T}\boldsymbol{P}\_{\boldsymbol{\Delta}}^{(i)}\boldsymbol{A} & \boldsymbol{B}^{T}\boldsymbol{P}\_{\boldsymbol{\Delta}}^{(i)}\boldsymbol{B} + \boldsymbol{P}\_{\boldsymbol{\mathcal{S}}}^{(i)} \end{bmatrix} \\ \boldsymbol{N}\_{\boldsymbol{\Delta}} = \begin{bmatrix} \boldsymbol{A}^{T}\boldsymbol{P}\_{\boldsymbol{\Delta}}^{(i)}\boldsymbol{A} & \boldsymbol{A}^{T}\boldsymbol{P}\_{\boldsymbol{\Delta}}^{(i)}\boldsymbol{B} \\ \boldsymbol{B}^{T}\boldsymbol{P}\_{\boldsymbol{\Delta}}^{(i)}\boldsymbol{A} & \boldsymbol{B}^{T}\boldsymbol{P}\_{\boldsymbol{\Delta}}^{(i)}\boldsymbol{B} \end{bmatrix} \\ \boldsymbol{N}\_{\boldsymbol{\mathcal{S}}} = \begin{bmatrix} \boldsymbol{\mathbf{0}} & \mathbf{0} \\ \boldsymbol{\mathbf{0}} & \boldsymbol{P}\_{\boldsymbol{\mathcal{S}}}^{(i)} \end{bmatrix}$$

The adaptive factor is

$$a^{(i)} = \mathfrak{d}\_{\mathbb{S}}^{2(i)} / \mathfrak{d}\_{\Delta}^{2(i)} \tag{5}$$

The two types of observations are divided into two groups, and we apply iterative calculation to each group using the adaptive least squares collocation method [53,54]. The specific process is as follows: <sup>1</sup> solve the least squares collocation configuration by formulas (2) and (3); <sup>2</sup> solve the posterior unit weight error *σ*ˆ (*i*) <sup>Δ</sup> , *σ*ˆ (*i*) *<sup>S</sup>* and adaptive factor *<sup>α</sup>*(*i*) by formulas (4) and (5); <sup>3</sup> determine whether |α(i) <sup>−</sup> 1| < <sup>ε</sup> works, where <sup>ε</sup> is a small constant, which is 0.001 in this paper. If the conditions are satisfied, the iteration stops; if the conditions are not satisfied, let *P*(*i*+1) <sup>Δ</sup> <sup>=</sup> *<sup>α</sup>*(*i*)*P*(*i*) <sup>Δ</sup> , modify the weight of the observed value and repeat steps <sup>1</sup> – <sup>3</sup> .

#### **3. Results and Discussion**

#### *3.1. The Current Vertical Crustal Deformation in the SYR*

The leveling vertical velocity ranges from −2.3 to 6.5 mm/a, and the error ranges from ±0.5 to ±2.1 mm/a, with an average error of ±1.2 mm/a (Figure 2a). The GNSS vertical velocity ranges from −3.6 mm/a to 4.9 mm/a, and the error ranges from ±0.1 mm/a to ±2.4 mm/a, with an average error of ±0.9 mm/a (Figure 2b). As Figure 2 shows, the two results are complementary in space. After fusion, the leveling observation increases the spatial density of the GNSS observation across the faults, and the GNSS observation makes the spatial distribution of the leveling observation more even. The vertical deformation velocity derived from these two datasets in the eastern margin of the QTP and its surrounding areas are consistent. Generally, the deformation in the SYR is dominated by uplifts. The uplift rates are large in eastern Tibet and the northern part of the rhombus block. Some regions have subsidence, including the Sichuan Basin and its periphery, the Jinghong–Pu'er– Jiangcheng area in southern Yunnan, and the Yunlong–Baoshan–Changning area.

In this paper, the least squares collocation method based on the multiquadric function is adopted to fuse the results of the two datasets with the adaptive Helmert iterative weighting. During the fusion, the two results constrain each other to ensure the comparability. The

velocity after fusion ranges from −3.1 mm/a to 6.3 mm/a (Figure 2c). We generate the contour and gradient maps from the fused vertical motion rates by the least squares collocation method (Figure 5). On the whole, most areas in the SYR show vertical deformation patterns dominated by uplifts. Specifically, the eastern Tibet and the northern part of the rhombus block in the SYR have rapid uplift, but the Sichuan Basin and its periphery, the Jinghong– Pu'er–Jiangcheng area in southern Yunnan and the Yunlong–Baoshan–Changning area show subsidence.

**Figure 5.** (**a**) Contour and (**b**) gradient map of the vertical motion velocity in the eastern margin of the Qinghai–Tibet Plateau and its surrounding areas. Solid black lines indicate the contours of vertical motion velocities in (**a**) and gradients in (**b**). The blue and red arrows are the vertical velocity after fusing the GNSS and the leveling data, respectively.

We combine the reported horizontal crustal movement features in the SYR [55] with our vertical deformation velocities to analyze the crustal movement features in the SYR: the vertical uplift is very obvious in the areas with strong horizontal crustal compression; the subsidence appears in the areas with extensive horizontal tension; the horizontal and vertical crustal movements are spatially coincided. These features reflect that the SYR was pushed and rotated by the QTP and blocked and squeezed by the South China block. In such an environment, the vertical crustal movement presents the contemporaneous occurrence of crustal compression, shortening and uplift and basin extensional subsidence. The detailed characteristics of the vertical deformation are as follows:


fault to the Xianshuihe fault and the Kangding–Yajiang–Litang–Batang area. The maximum velocity, 4 mm/a, is observed at the eastern section of the Ganzi–Yushu fault, but it rapidly drops to about 2 mm/a at the junction with the Xianshuihe fault, and then gradually rises to 3 mm/a~4 mm/a along the Xianshuihe fault, showing a segmentation feature. The results of the Gongga Mountain area, a part of the southeastern section of the Xianshuihe fault, are consistent with those of Wang [26] (1.8 mm/a~3.6 mm/a) and Wang et al. [57] (about 4 mm/a), but are smaller than those of Su et al. [34] (5 mm/a~6 mm/a) and Wang et al. [58] (7.4 mm/a). The overall uplift also confirms that during the plastic flow, the material under the QTP crust was compressed near the boundary due to the blockage of the hard Sichuan Basin and exerted great upward fluid pressure on the elastic upper crust [26,59].


To further explore the vertical deformation features of some important sections in the SYR, we extract three profiles, AA , BB and CC , in the eastern Bayankala block, the eastern part of the Sichuan–Yunnan border and the southern Yunnan area, respectively (Figure 1). Vertical motions along these three profiles are shown in Figure 6. Profile AA starts from the Aba area of Sichuan and trends NW–SE to Chengdu, with a total length of about 300 km, crossing the Longriba and Longmenshan faults. Profile BB starts from Yanyuan and trends NE–SW to Yuexi, with a total length of about 370 km, crossing the Anning River and the Mabian-Ebian fault. Profile CC' starts from Mojiang and trends NE–SW to Qujing, with a total length of about 380 km, crossing the Red River and Xiaojiang faults. The vertical deformation characteristics reflected by the profiles are as follows:

**Figure 6.** Three profiles, AA , BB and CC , showing the vertical motions in the SYR. Profile locations are shown in Figure 1.


and Yunnan experienced compression and rising caused by the eastward movement of the northwestern Sichuan block. There may be risks of regional strong earthquakes.

(3) Profile CC shows that the vertical movement rate of the east side of the Xiaojiang fault is 0.5 mm/a~1.5 mm/a, and that of the west side is 0.5 mm/a~2.0 mm/a, which are not significantly different. The vertical movement rate on the south side of the Red River fault zone is 0.5 mm/a~1.5 mm/a, and the vertical difference between the two sides is not significant.

#### *3.2. Kinetic Mechanism*

The SYR is located on the southeastern margin of the QTP, where the lateral extrusion of the QTP caused by the northward push of the Indian plate was blocked by the stable South China block, resulting in the plateau eastern margin continuously uplifting [31,44]. Such a movement pattern has seen no great changes since the Pliocene [17,60]. The collision between the Indian plate and the QTP plate not only caused the crust vertical movement, but also controlled the fault motion mode and the horizontal movement of each block. Most of the material from the QTP moved southward. The rheological strength of the middle and lower crust in the SYR is weaker than that of other crusts. The crust thickness in the region is 20~30 km less than that of the Sichuan Basin and South China block, and the elevation difference is about 3000 m, causing lateral pressure difference. Under this pressure difference, the weak rheological layer in the middle and lower crust flowed, which enhanced the clockwise rotation of the lithosphere in the southeastern QTP around the eastern Himalayan tectonic junction. Such rotation caused the deformation of the upper brittle crust and the movement of the faults, which finally showed as a huge left-lateral strike-slip arc fault structure [61].

The material from the eastern margin of the QTP flowed to three directions. Some moved northeastward and interacted strongly with the Ordos Basin. Some moved to western Sichuan and strongly coupled with the Longmenshan fault, and the other moved to Yunnan after a sharp turn. In the Longmenshan area north to the Xianshuihe fault, the material accumulated in the middle and lower crust after being blocked by the hard crust of the Sichuan Basin, resulting in sharp crust thickening and strong surface uplifts, forming the towering Longmenshan Mountains [14,16,62–64]. In the SYR south to the Xianshuihe fault, the crust was relatively soft and the material flow channels in the middle and lower crusts were open. Affected by the drag of this pipeline flow, the SYR had a clockwise rotation, and the crust gradually thickened on a scale of thousands of kilometers. A large-scale of ground surface had uplifted, forming the broad southeastern margin of the QTP, featured by similar geomorphology and blurred boundary [14,16,17,63,64].

#### *3.3. Distribution of Lower Crustal Flows*

Deep geophysical exploration revealed that the middle and lower crusts on the eastern margin of the QTP have seismic low-velocity [21,23,28], high Poisson's ratio [19,27] and high conductivity/low resistivity [25,27,28]. These conditions favor partial melting, especially if there are salt-bearing fluids. Partial melting reduces the strength of the middle and lower crusts and promotes laminar flow [61].

In this study, the finite difference method was adopted to calculate the effective elastic thickness of the lithosphere in the SYR from the EIGEN6C4 Bouguer gravity anomaly and SIO V15.1 topographic data (Figure 7). The effective elastic thickness of the lithosphere reflects the relation between the mechanical strength of the lithosphere [65,66] in the geological time scale (>1 Ma) and the composition of crust–mantle matter, thermal state, and the crust–mantle coupling extent [67,68]. Combined with the reported results in seismology and magnetotelluric fields, we infer that the extruded material from the QTP was blocked by the Sichuan Basin, and accumulated in the Litang–Daocheng–Lijiang area, which showed as prominent surface uplifts and large-scale internal crustal deposits with low velocity and high conductivity/low resistivity. A small part of the material may have moved northward through the Kangding–Daofu area in the Xianshuihe fault, but most of

the material moved southward. After being blocked by the southern Yunnan block, a part of the material flowed to Tengchong in the southwest and the rest flowed southeast to the Panzhihua–Dongchuan area.

**Figure 7.** The effective elastic thickness of the lithosphere in the SYR, calculated from the EIGEN6C4 Bouguer gravity anomaly and SIO V15.1 topographic data.

The vertical crustal movement features in the SYR (Figure 2) suggest that the southeastern margin of the QTP was blocked by the Sichuan Basin, so most of the material flowed southward along the arc-shaped channel into southern Yunnan, some of which was blocked by the Lijiang-Xiaojinhe fault. The rest continued to move southward, but was blocked by the southern Yunnan block. Some crossed the Red River fault and flowed to Pu'er in the southwest, and the rest crossed the Xiaojiang fault and flowed to Liupanshui in the southeast. In the SYR, due to the blockage of the Lijiang-Xiaojinhe fault, the material crossed the Jinsha River–Red River fault and the Xiaojiang fault in the southwest and southeast directions, respectively.

#### **4. Conclusions**

Fusing the precise leveling data of the 1970s–2010s and the GNSS data of the 1990s–2010s by the least squares collocation method, we obtain the background field of the vertical crustal movements in the SYR, which reflects the features of the current crustal vertical movements in this region. These features coincided with the neotectonic movement background that the plateau and mountains uplifted and basins subsided. We come to the following conclusions:


ment in the SYR presents the contemporaneous occurrence of crustal compression, shortening and uplift and basin extensional subsidence.


**Author Contributions:** Conceptualization, Y.Z. and C.X.; methodology, Y.Z. and Z.Z.; software, Z.Z., H.L. and S.Z.; validation, Y.Z. and Z.Z.; formal analysis, Y.Z.; writing—original draft preparation, Y.Z. and C.X.; writing—review and editing, Y.Z., C.X. and Z.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported by the National Natural Science Foundation of China under Grants No. 41721003 and No. 42074007.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data used in this manuscript can be archived at the Figshare website (https://doi.org/10.6084/m9.figshare.18871442.v1, last accessed on 1 February 2022) for free public access.

**Acknowledgments:** We appreciate the GPS observations from the Crustal Movement Observation Network of China and the China Continental Tectonic Environment Monitoring Network. We thank Ming Hao at the Second Monitoring Center, China Earthquake Administration, for providing the leveling data. All figures were constructed using Generic Mapping Tools (GMT) version 4.5.1 [69].

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

