Introduction

Southwest China (Fig. 1), which is located at the southeastern margin of the Tibet Plateau, is one of the regions with the most intense tectonic and seismic activity in Chinese Mainland due to the collision of the Eurasian and Indian plates (Royden et al. 2008; Shen et al. 2005; Yin and Harrison 2000). Furthermore, the climate in this region has a subtropical plateau monsoon that is characterized by abundant precipitation with a significant annual cycle (Li et al. 2015; Xiao et al. 2014), resulting in the impressive seasonal oscillation of crustal loading deformation. Following the development of satellite technology and data assimilation methods, such as global positioning system (GPS), gravity recovery and climate experiment (GRACE) and surface loading model (SLM), the elastic response of solid earth to mass variations can be determined from an independent perspective. Several researchers have conducted comparative analyses on the seasonal deformation obtained by the three methods at both global (Tregoning et al. 2009; Tesmer et al. 2011) and regional scales (Davis et al. 2004; Fu and Freymueller 2012; Nahmani et al. 2012; Ferreira et al. 2019; Liu et al. 2018), wherein the three technologies generally showed a good consistency. Sheng et al. (2014), Hao et al. (2016), Pan et al. (2018), and Zhan et al. (2020) found that the seasonal vertical deformation time series acquired by GPS was highly correlated with that acquired by GRACE in this region, indicating that the annual variation of land water is the main reason for the seasonal deformation. Zhan et al. (2017) compared the seasonal deformation captured by GPS and precipitation in this region, and reported strong temporal and spatial consistency. However, Gu et al. (2017) and Yan et al. (2019) noted that the amplitude of the mass loading signal extracted by GPS in this area was larger than that extracted by GRACE and SLM. Currently, quantitative analysis on the reasons for the differences between the three products is lacking in this region. In addition, owing to the strong tectonic activity, long-term vertical crustal movement in Southwest China is significant. Distinguishing the crustal tectonic movement and loading deformation in the vertical GPS time series is essential to understand the present tectonic activities in this area.

Fig. 1
figure 1

Distribution of continuous global positioning system (cGPS) stations (blue dots). The red line indicates the boundary of the block (Luo et al. 2012; Zhang et al. 2003); I: South China block; II: Sichuan–Yunnan rhombic block; II1: West Sichuan block; II2: Central Yunnan block; III: Southwest Yunnan block

Full size image

Following the implementation of the tectonic and environmental observation network in China (CMONOC II) and GRACE follow-on (GRACE-FO) missions, there have been 10 years of continuous GPS (cGPS) observations and time-variable gravity data. These long-term observations can allow for the objective comparison of the three methods in capturing seasonal loading deformation, especially with respect to understanding the systematic difference. In addition, these long-term data can provide relatively reliable long-term tectonic movement information. In this study, we systematically collected and processed three groups of vertical coordinate time series produced by GPS, GRACE, and SLM at 41 cGPS stations. The remainder of this paper is organized as follows. In “Data and methods”, the data and methods used are discussed. In “Results”, the consistency of the three products is quantitatively analyzed using the correlation coefficient (CC) and root mean square (RMS). In “Discussion”, the systematic differences in the seasonal loading deformation among the three products are discussed, and the long-term crustal motion is presented and analyzed. In “Conclusions”, the main conclusions drawn are summarized.

Data and methods

GPS data

The GPS observation network (Fig. 1) in Southwest China reached a certain density in the past decade following the implementation of CMONOCII. These continuous GPS stations had a 10-year observation time span from 2010 to 2020. We employed GAMIT/GlOBK to process the GPS data, in which the loading effects of ocean tides and solid earth tides were corrected using the corresponding models of FES2004 and IERS2003. We used antenna calibrations issued by National Geodetic Survey (https://geodesy.noaa.gov/ANTCAL/) for both GPS receiver and satellite antennas in the data processing. Then, the daily solutions were transformed to the ITRF2008 reference frame using some frame sites (Altamimi et al. 2011). Zhan et al. (2017) provide more detailed related information regarding this. Offsets caused by instrument updates, coseismic deformation, and unknown reasons were eliminated from the vertical coordinate time series. The effects caused by the mass loading of atmosphere, terrestrial water storage, and non-tidal ocean were not eliminated. Particularly, in the definition of ITRF2008, the origin of ITRF2008 was determined by the mean of the satellite laser ranging time series, which is sensitive to the center of mass of the whole earth system (CM). However, Dong et al. (2003) reports that in the long-term scale, ITRF is represented as a CM frame, while for short time scales or seasonal scales, it is characterized as a Earth’s center of figure (CF) frame.

Surface loading model data

The surface deformations caused by hydrological loading (HYDL), non-tidal atmospheric loading (NTAL), and non-tidal ocean loading (NTOL) were considered in the surface loading model for comparison with GPS datasets. Hydrological changes (including soil moisture, water equivalent of snow cover, and water content in the canopy) exert varying water pressure, which leads to crustal deformation (van Dam and Wahr 1998; Dong et al. 2002). Several studies (Nahmani et al. 2012; Ferreira et al. 2019) have suggested that variation in land water loading is the dominant factor affecting mass redistribution at the seasonal scale. Moreover, atmospheric pressure variation is also a non-negligible loading factor that results in crustal elastic deformation (van Dam et al. 1994). In this study, we used the International Mass Loading Service (IMLS) (Petrov 2015) to compute the surface deformation driven by mass redistribution. The numerical weather model, GEOS-FPIT, with a resolution of 0.5° × 0.625° × 1 h (time resolution was reduced to 3 h for loading computation) (Rienecker et al. 2008) was used to estimate the deformation resulting from HYDL and NTAL. Furthermore, the response of the elastic earth to changes in ocean bottom pressure is referred to as NTOL, and calculated via IMLS based on MPIOM06 (Jungclaus et al. 2013). In the following analysis, all the loading datasets were down-sampled to daily time series, and linked to the Center of figure (CF) reference frame for consistency with the GPS data. In addition, we combined HYDL, NTAL and NTOL to deduce the total loading deformation.

GRACE data

The time-variable gravity signals obtained by the GRACE satellite have been extensively used to monitor mass redistribution, such as variations in terrestrial water reserves, glacier mass changes, and the movement of the ocean and atmosphere. We employed the latest Release 06 (Level-2) GRACE/GRACE-FO solutions from 2010 to 2020, which was provided by the Center for Space Research, the University of Texas at Austin (Save et al. 2018), to estimate the vertical loading deformation in Southwest China, where these products were in the form of spherical harmonic coefficients up to a degree and an order of 60. The degree-1 terms estimated by Sun et al. (2016) were added. Moreover, because of the uncertainty of the C20 coefficients, we replaced them with the value estimated from analyzing the satellite laser ranging data (Cheng and Ries 2017). Given that the effect of non-tidal atmosphere and ocean would be determined from GPS datasets, to maintain data consistency, Level-2 GAC products (Bettadpur 2018), representing the sum loading of the non-tidal atmosphere and non-tidal ocean, were also employed. We used the equation derived by Kusche and Schrama (2005) to convert the harmonic coefficient to the surface vertical deformation:

$${\text{def}}. = R\sum\limits_{{n = 1}}^{{60}} {\frac{{h^{\prime}_{n} }}{{1 + k^{\prime}_{n} }}\sum\limits_{{m = 0}}^{n} {\left[ {\Delta C_{{nm}} \left( t \right)\cos m\lambda + \Delta S_{{nm}} \left( t \right)\sin m\lambda } \right]} P_{{nm}} \left( {\cos \theta } \right)} ,$$
(1)

where def. is the elastic deformation, R is the mean radius of the Earth, and Cnm and Snm are the residual values after removing the average spherical harmonic coefficients of 2010–2020, thus highlighting the periodic signal in this study. Furthermore, a simplified order-convolution smoothing approach (DDK5) (Kusche et al. 2009) was applied to weaken the uncertainty at high degrees (Wahr et al. 2004). Pnm represents degree n and order m values of the fully normalized Legendre functions. \(\theta\) and \(\lambda\) are the co-latitude and longitude, respectively. \(h^{\prime}_{n}\) and \(k^{\prime}_{n}\) are the load love numbers relating to the horizontal and radial directions, respectively. The love number adopted in this study was derived by Farrell (1972), where the gaps were repaired by linear interpolation; moreover, the degree of one love number was linked to the CF frame according to Blewitt (2003).

Time series processing method

The vertical displacement time series derived from GPS, GRACE, or SLMs usually contains different time–frequency signals, including noise, seasonal loading deformation, and long-term signals. In this study, we used ensemble empirical mode decomposition (EEMD) to decompose and reconstruct the time series into noise, seasonal, and long-term trend signals. EEMD is a time–frequency analysis method that self-adaptively and efficiently decomposes non-stationary and non-linear time series into a sequence of intrinsic mode functions (IMFs) and residual terms without parameter settings (Huang et al. 1998; Wu and Huang 2011). The general EEMD algorithm is expressed as follows:

$$S\left( t \right) = \sum\limits_{{i = 1}}^{m} {\overline{{x_{i} }} \left( t \right) + r_{m} \left( t \right),t = 1,2, \ldots ,}$$
(2)

where \(\overline{{x_{i} }} ~\) is the ith IMF, \(r_{m} \left( t \right)\) is the residual signal after m decomposition and \(S\left( t \right)\) is the original signal.

In this study, we employed the Hurst parameter (H) to classify the different time–frequency IMFs. The Hurst parameter is typically used for signal recognition. For example, Montillet et al. (2013) reconstructed the GPS-derived IMFs with H ≤ 0.5 into white noise signals. Peng et al. (2016) decomposed the GPS vertcial time series into several IMFs using EEMD, and then, summed the IMFs with 0 ≤ H < 1.1 as a colored noise signal. According to the statistics of GPS stations in mainland China, Peng and Dai (2016) found that H = 1.8 is the optimal boundary value between the trend and seasonal signals. In this study, we processed the time series as follows.

  1. 1.

    The time series was decomposed into several IMFs and residual terms using EEMD;

  2. 2.

    The Hurst parameter of each IMF was estimated;

  3. 3.

    IMFs with 0 ≤ H < 0.6 were reconstructed into the noise signal, and the sum of the IMFs with H > 1.8 and the residual terms were reconstructed into the long-term trend signal (Peng and Dai 2016). The rest of IMFs were reconstructed into the seasonal loading deformation signal.

  4. 4.

    The linear rates were estimated from the long-term trend signal using the least square method.

Figure 2 shows an example of time series processing using EEMD method. Compared with the function model, the result shows that EEMD has advantages in characterizing the signal of which amplitude changes with time.

Fig. 2
figure 2

Vertical coordinate time series at the YNRL station. a Raw time series; b intrinsic mode functions (IMFs) decomposed by ensemble empirical mode decomposition (EEMD); c reconstructed deformation signal. Gray point, red curve, and green curve denote GPS daily solution, reconstructed seasonal and trend signal, respectively. d Modeled result. Gray point is GPS daily solution, blue curve shows filter results using the linear + annual-cycle + semi-annual-cycle movement model, and black line is the linear trend

Full size image

Results

Components analysis of seasonal deformation

Variations in mass loadings, such as HYDL, NTAL, and NTOL, lead to seasonal crustal oscillations. We derived the temporal variations of different environmental mass loading components, which were estimated by SLM at 41 GPS stations. Figure 3 illustrates the stack-average variation curves (after removing the noise and long-term trend signals) of HYDL, NTAL, and NTOL that are represented by blue, red, and green lines, respectively. Furthermore, we calculated the annual amplitudes and phases of the different loading components using the least squares method. The results show that vertical deformation caused by NTOL (Fig. 3; green line) was the weakest, with an amplitude of less than 0.3 mm, as the research area is far from the coastline. The red curve in Fig. 3 shows that the deformation caused by NTAL, with an amplitude of 2.3 mm, exhibited significant seasonal fluctuation, and an upward peak appeared in late June or July. The upward peak of the HYDL signal (Fig. 3; blue line) was observed in March or April, approximately 3 months before NTAL. The amplitude of HYDL was 5.2 mm, approximately 2.3 and 17.3 times larger than those of NTAL and NTOL, respectively, indicating that the hydrological loading changes were the main factor causing the seasonal crustal deformation in this area.

Fig. 3
figure 3

Loading components in Southwest China. Stack-average deformation of 41 stations yielded by hydrological loading (HYDL), non-tidal atmospheric loading (NTAL), and non-tidal ocean loading (NTOL) represented by the blue, red, and green lines, respectively [inferred from the surface loading model (SLM)]

Full size image

In addition, we collected monthly surface precipitation data from the China Meteorological Science Data Center for further verification. The precipitation showed significant seasonal variations, more than 85% of the total rainfall occurred in the rainy season spanning from May to October, and the rest occurred in the dry season spanning from November to April of the succeeding year. As shown in Fig. 3, there was an evident negative correlation between rainfall and HYDL (blue line). That is, the ground sunk as the hydrological loading increased during the rainy season, and the ground surface rose following water loss during the dry season. In addition, a slight phase lag was clearly exhibited by the color stripe in Fig. 3, suggesting that the peak representing the deformation caused HYDL approximately 1.4–1.8 months later than precipitation. We attribute this lag to the fact that the loading deformation reflects the variations in cumulative land water, rather than instantaneous precipitation (Zhan et al. 2017; Birhanu and Bendick 2015).

In addition, as shown in Fig. 4, we calculated the deformations caused by the sum NTAL and NTOL, and HYDL using the GAC and GSM solutions of the GRACE products (Bettadpur 2018), respectively. The results indicate that the NTAL + NTOL obtained by SLM and GRACE showed high consistency in both amplitude and phase (Fig. 4a). The average amplitudes and phases of the 41 stations are 2.7 mm and the 171st day for SLM, and 2.8 mm and the 172nd day for GRACE, respectively. In terms of HYDL (Fig. 4b), the amplitude estimated by GRACE (6.5 mm) was larger than that estimated by SLM (5.2 mm), possibly owing to unmodeled components such as groundwater and deep soil moisture in SLM. The phase of GRACE is the 84th day, which is consistent with that of SLM (83rd day).

Fig. 4
figure 4

Phases and amplitudes of different loading components. a Non-tidal atmospheric loading (NTAL) + non-tidal ocean loading (NTOL); b hydrological loading (HYDL)

Full size image

Consistency analysis of the seasonal deformation estimated by GPS, GRACE, and SLM

To quantitatively evaluate the consistency of different time series products, we calculated the CC and reduction in the RMS of GPS/GRACE and GPS/SLM, where CC was sensitive to the signal phase, and RMS reduction could allow for overall assessment (Gu et al. 2017). The RMS reduction is expressed as follows:

$${\text{RMS}}_{{{\text{red}}}} = \frac{{{\text{RMS}}\left( {{\text{GPS}}} \right) - {\text{RMS}}\left( {{\text{GPS}} - {\text{GRACE}}/{\text{SLM}}} \right)}}{{{\text{RMS}}\left( {{\text{GPS}}} \right)}} \times 100.$$

Higher CC and RMS reduction values indicate more optimized consistency among the different techniques. Furthermore, we employed a function model, containing linear trends in addition to semiannual and annual harmonic signals, to characterize the deformation time series of GPS, GRACE, and SLM, where the model parameters were estimated using the least squares method (Table 1 and Fig. 2d).

Table 1 Annual amplitude and phase estimated by GPS, GRACE, and surface loading model (SLM)
Full size table

Figure 5 shows the time series (after removing the trend signal), which were estimated from GPS, GRACE and SLM at three sample GPS locations. The vertical time series among GPS, GRACE, and SLM exhibited consistent seasonal variations, suggesting that the seasonal oscillations in the GPS time series are mainly caused by mass redistribution. YNLC (top panel of Fig. 5), located in southwest of Yunnan province, yielded the highest CC: 0.90 and 0.89 for GPS/GRACE and GPS/SLM, respectively. The RMS reductions upon subtracting GRACE and SLM from the GPS signals were 55.8% and 47.5%, respectively. The annual amplitudes estimated from GPS, GRACE, and SLM by the least squares method were 11.6, 9.3, and 7.3 mm, respectively (Table 1). Figure 5 also shows the time series at SCYX station, which is located in southern Sichuan. This station, to some extent, represented the performance of most stations, where the CC was close to the average value (0.82 and 0.81) of 41 stations; the values were 0.82 and 0.81 for GPS/GRACE and GPS/SLM, respectively. The RMS reduction was 42.9% for GPS–GRACE and 39.7% for GPS–SLM. The annual amplitudes of GPS, GRACE, and SLM were 7.3, 5.7, and 4.8 mm, respectively (Table 1). Figure 5 (bottom panels) shows the stations of YNMZ with the poorest performance in the consistency analysis, which may be caused by an unknown local loading process (e.g., river or reservoir), rather than the relatively smooth signals captured by GRACE or SLM. The CC values of YNMZ were 0.61 and 0.66, respectively, and the RMS reductions were 17.0% and 24.8%, respectively, for GPS/GRACE and GPS/SLM owing to a visible phase delay in GRACE and SLM.

Fig. 5
figure 5

Seasonal vertical deformation time series derived by global positioning system (GPS) stations, gravity recovery and climate experiment (GRACE), and surface loading model (SLM) at YNLC, SCYX, and YNMZ stations

Full size image

Figure 6 illustrates the spatial distribution of CC at 41 stations for different seasonal deformation products. All the stations showed consistency among the three techniques. The CC of 27 and 26 stations exceeded 0.80 for GPS/GRACE (Fig. 6a) and GPS/SLM (Fig. 6b), respectively. The mean values were 0.82 and 0.81 for GPS/GRACE and GPS/SLM, respectively. Figure 6 also shows the spatial distribution of RMS reductions after subtracting the GRACE or SLM signals from the GPS time series. All stations showed a positive response after suppressing the seasonal signals. The mean RMS reductions for GPS–GRACE in Fig. 6a were 41.3% (> 40% for the 26 stations, 30–40% for the 9 stations, and 10–30% for the remaining 6 stations). Furthermore, the mean RMS reduction for GPS–SLM in Fig. 6b was 38.0% (> 40% for the 18 stations, 30–40% for the 16 stations, and 10–30% for the remaining 7 stations).

Fig. 6
figure 6

Correlation coefficient (CC) and root mean square (RMS) reduction at 41 stations. a Global positioning system (GPS)/gravity recovery and climate experiment (GRACE); b GPS/surface loading model (SLM). The color of the circle represents the RMS reduction, and the number in the circle is CC × 100

Full size image

Discussion

Differences among the three products

Although GPS, GRACE, and SLM show consistency to a certain extent in monitoring seasonal deformation, there are still evident systematic differences between the three products. For example, the amplitude observed by GRACE and SLM is smaller than that observed by GPS, and the phase lags behind the GPS (Table 1). Figure 7 shows the amplitudes and phases of the seasonal loading deformation calculated using GPS, GRACE, and SLMs at 41 stations. The average annual amplitude of GPS was 9.7 mm (ranging from 5.7 to 14.2 mm), which was 1.3 and 1.6 times higher than that of GRACE (7.4 mm; ranging from 5.3 to 9.8 mm) and SLM (6.1 mm; ranging from 4.0 to 8.6 mm), respectively. Furthermore, a systematic difference was observed in the phase. The mean phases of the three vertical displacement time series, representing the time when the upward peak appeared, were the 91st, 107th, and 113th days for GPS, GRACE, and SLM, respectively, indicating that there was a 16, 22, and 6 day phase lag between GPS and GRACE, GPS and SLM, GRACE and SLM, respectively.

Fig. 7
figure 7

Amplitudes and phase of total loading deformation estimated by global positioning system (GPS) stations, gravity recovery and climate experiment (GRACE), and surface loading model (SLM) at 41 stations

Full size image

As shown in Fig. 7, the amplitudes estimated by GPS were larger than those estimated by GRACE and SLM possibly due to the uncertainty of the three products. GPS obtains surface deformation with the point measurement mode, which is susceptible to local environmental loading (Nahmani et al. 2012). For example, Gu et al. (2017) noted that the GPS station of YNTH is affected by variations in the nearby river level. The seasonal deformation inferred by GPS is usually contaminated by the thermal expansion of the bedrock and concrete pillar (Fang et al. 2014; Yan et al. 2009). We employed the thermal expansion model derived by Yan et al. (2009) and the global temperature data provided by the National Centers for Environmental Prediction to estimate the thermoelastic deformation at 41 GPS stations. The results show that the average annual amplitude of the thermoelastic deformation was 0.4 mm. The seasonal signals may be affected by errors in GPS data processing. Zhan et al. (2017) indicated that the difference in annual amplitude under different reference frames was approximately 1 mm in this region. Moreover, the annual amplitude of GPS might be overestimated because of the influence of draconitic errors of satellite orbit, multipath effect, etc. (Gu et al. 2017; Rodriguez-Solano et al. 2014; Ray et al. 2007).

In terms of GRACE, the spatial smoothing filter is required to suppress the “north–south” strip error in the processing of GRACE geopotential coefficients. Abundant rainfall means a strong HYDL signal in Southwest China. Coincidentally, regions characterized by such high hydrological signals are susceptible to spatial filtering, causing signal leakage (Chen et al. 2006; Swenson and Wahr, 2011). Using the GSM and GAC monthly solutions of GRACE, we acquired the average annual amplitudes of HYDL and NTAL + NTOL under different filter radii, as shown in Fig. 8a. The results show that, by enhancing the filtering strength, the amplitude of the HYDL is gradually attenuated, indicating that the spatial filtering in this region may cause signal leakage to the adjacent area. The effect of spatial filtering on the annual amplitude of NTAL + NTOL was negligible (< 0.15 mm). In addition, the land water storage estimated by SLM only contains variations in shallow soil moisture, snow depth water equivalent, and plant canopy surface water, but lacks deep soil moisture and groundwater storage changes. Therefore, the land water storage will be underestimated, especially in regions with high HYDL (Scanlon et al. 2019; Mo et al. 2016). Then, the weakened signal attenuates the annual deformation amplitude. In general, other geophysical signals and data processing errors in GPS, the leakage-out errors in GRACE, and the underestimation of land water storage in SLM are jointly responsible for the difference in amplitude among the three products.

Fig. 8
figure 8

Reason for the difference among the three products. a Relation between the amplitude of gravity recovery and climate experiment (GRACE)-dependent loading deformation and filter radius; b simulative signals of hydrological loading (HYDL) and non-tidal atmospheric loading (NTAL) + non-tidal ocean loading (NTOL); c simulative signals of combined loading

Full size image

In addition to the amplitudes, a systematic phase difference among the three products is shown in Fig. 7. The seasonal deformation in this area is the result of the aforementioned signals of HYDL, NTAL, and NTOL, which we referred to as Shydl, Sntal, and Sntol, respectively, in this section. Because the NTOL is faint, we combined the signals of NTAL and NTOL and marked it as Sa&o in the following content. The seasonal loading deformation in this area was the superposition of Shydl and Sa&o, as shown in Fig. 8b. If Shydl is underestimated (e.g., leakage-out error in GRACE and unmolded components in SLM) or overestimated (e.g., GPS processing errors) when the two annual signals are superimposed, the amplitude of the synthesized signal will be weakened or amplified, respectively, and the phase will change at the same time. We simulated the two signals using a sine function, and set the amplitude (Aa&o) and peak (\(\varphi _{{{\text{a}}\& {\text{o}}}}\)) of Sa&o as 2.8 mm and the 172nd day, and set the peak of Shydl on the 84th day (these parameters were estimated by GRACE and SLM, as shown in “Results”). Then, we set the amplitude of Shydl to 5.2 mm (to simulate SLM), 6.5 mm (to simulate GRACE), and 9.4 mm (to simulate GPS), and combined them with Sa&o to obtain the superimposed signal (Stotal), as shown in Fig. 8c. The results show that the phase delays between GPSsimu/GRACEsimu, GPSsimu/SLMsimu, and GRACEsimu/SLMsimu were 8, 14, and 5 days, respectively. This indicates that approximately 50% (8/16 days), 64% (14/22 days), and 83% (5/6 days) phase difference between GPS/GRACE, GPS/SLM, and GRACE/SLM are caused by the underestimation or overestimation of Shydl. The remaining unexplained difference may be caused by the effects of local loading, the thermal expansion of the crust, and the error in the data processing.

Long-term crustal movement in Southwest China

As an elastic body, variations in the surface loading cause crustal deformation, primarily in the radial direction (Yan et al. 2019). Owing to the environmental mass redistribution (including hydrology and atmosphere), the vertical deformation in Southwest China tends to exhibit a downward tendency in addition to seasonal oscillation. Figure 9a shows the long-term surface mass changes and the resulting deformation estimated by GRACE from 2010 to 2020. In particular, we did not consider the impact of glacial isostatic adjustment in data processing, as Pan et al. (2016) and He et al. (2017) (including their citations) show that the correction model in this region is still controversial, owing to insufficient space geodetic data to restrain it. We assumed that other nontectonic geophysical factors (e.g., polar motion (King and Watson 2014), mantle anelasticity (Ding and Chao 2017)) did not influence loading deformation. Figure 9a shows that, in the study area, except for the northwest, the material in other regions increased at a rate of approximately 5–20 mm/year (in equivalent water height), with a peak increment in the eastern part of Sichuan. Furthermore, long-term vertical loading deformation rates at 41 GPS stations were estimated with the deformation time series acquired by GRACE. The rate ranges from − 0.79 to + 0.1 mm/year, with an average of − 0.58 mm/year, indicating that most areas were under subsidence due to the increased mass loading (Fig. 9a, white vector).

Fig. 9
figure 9

Three-dimensional deformation in Southwest China. a Long-term surface mass variations (in equivalent water height, the background color), the linear vertical motion rate estimated by GRACE (white vector) and GPS (blue and red vector); red line represents the block boundary. b Vertical crustal motion rates by subtracting the rates of GRACE from the linear rates of GPS. F1: Anninghe fault; F2: Xiaojiang fault; F3: Red river fault; F4: Jinshajiang fault; F5: Xiaojinhe fault. I: South China block; II: Sichuan–Yunnan rhombic block; II1: West Sichuan block; II2: Central Yunnan block; III: Southwest Yunnan block. c Horizontal velocity field provided by Wang and Shen (2020)

Full size image

In addition to mass loading, vertical motion is controlled by tectonic activities. The GPS records the vertical displacement caused by tectonic activity and mass variation simultaneously. We obtained the crustal motion rates by subtracting the mass loading deformation rates estimated by GRACE from the linear trend rates of the GPS (Fig. 9b). The corrected rates at 41 stations range from − 3.2 to 2.9 mm/year. As shown in Fig. 9b (three anomalous stations are not presented as they might be affected by local environmental loads), stations in the Sichuan–Yunnan rhombic block (block II in Fig. 9b) that are bounded by the Red river, Anninghe, Xiaojiang, and Jinshajiang faults are uplifted. In particular, the uplift at the Central Yunnan block (block II2 in Fig. 9b) and its eastern boundary are significant, with a upward rate of 1–3 mm/year, which is consistent with the value estimated by precise leveling (Su et al. 2018; Hao et al. 2014). In contrast with the uplifting points in the northeastern wall of the Red river fault, most of the points in the southwestern wall (block III in Fig. 9b) are under subsidence, which may indicate that the Red river fault is the constraining factor for vertical crustal motion in this area.

In addition, to analyze and explain the characteristics of vertical deformation in Southwest China, we collected the GPS horizontal velocity (Fig. 9c) provided by Wang and Shen (2020). The results show that the horizontal velocities rotate at the eastern Himalayan syntaxis; a part of which flows along the southeast direction to the east boundary of Sichhuan–Yunnan rhombic block. Affected by the Xiaojiang fault, the movement component along the fault causes a sinistral strike-slip motion, and the component perpendicular to the fault is blocked by the stable South China block, possibly causing crustal compression and uplift at the eastern boundary of the Central Yunnan block. Similarly, the horizontal crustal movement is affected by the northwest (NW) striking Red river fault. The velocity component along the fault direction causes a dextral strike-slip, and the component perpendicular to the fault direction uplifts the NE wall of the fault. Furthermore, the horizontal velocity in Southwestern Yunnan block (Fig. 9b, block III) featured with extensional movement in EW direction, which coincides with the signs of subsidence in the vertical direction.

Conclusions

The surface environmental mass redistribution and tectonic activity cause vertical deformation. Using GPS, GRACE, and SLM, we gained insight into the seasonal and long-term vertical deformation in Southwest China, and discussed the systematic differences in the seasonal scale among the three products. The main conclusions drawn are as follows.

  1. 1.

    The seasonal deformation in Southwest China is mainly caused by the annual variations of HYDL, of which the annual amplitudes are approximately 2.3 and 17.3 times that of NTAL and NTOL, respectively.

  2. 2.

    The CC of seasonal deformation captured by the three techniques are 0.82 for GPS/GRACE and 0.81 for GPS/SLM. After subtracting the results of GRACE and SLM from the GPS time series, the average RMS reductions are 41.3% and 38.0%, respectively.

  3. 3.

    The annual amplitude and phase estimated by GPS, GRACE, and SLM show systematic differences. The annual amplitude estimated by GPS is 1.3 and 1.6 times those of GRACE and SLM, respectively. The uplift peak of the GPS was approximately 16 and 22 days earlier than that of GRACE and SLM, respectively. The underestimation or overestimation of HYDL signals can explain approximately 50%, 64%, and 83% phase lag between GPS and GRACE, GPS and SLM, GRACE and SLM, respectively.

  4. 4.

    The vertical crustal motion in Southwest China is block-dependent. The Central Yunnan block is uplifting with rate in 1–3 mm/year, and the Southwest Yunnan block is undergoing subsidence with rate in − 2 to 0 mm/year.