Seasonal and long-term vertical land motion in Southwest China determined using GPS, GRACE, and surface loading model

Owing to the intense tectonic activity and significant seasonal surface mass change, Southwest China is characterized by noticeable vertical land movement. We determined the vertical movement of Southwest China using 10 years of data from 41 continuous global positioning system (GPS) stations, gravity recovery and climate experiment (GRACE), and surface loading model (SLM). The annual variation in hydrological loading is the main factor causing the seasonal oscillation of surface deformation in Southwest China. Seasonal deformations captured by GPS, GRACE, and SLM are consistent to a certain extent, and the correlation coefficients between GPS/GRACE, and GPS/SLM are 0.82 and 0.81, respectively. After deducting the results yielded by GRACE and SLM from the GPS time series, the average reductions in root mean square were 41.3% and 38.0%, respectively. However, some systematic differences were observed among the annual amplitudes and phases of the seasonal deformations among the three products. For example, the average amplitudes estimated by GRACE and SLM were 7.4 mm and 6.1 mm, respectively, which were smaller than the amplitude deduced from GPS (9.7 mm). Furthermore, mean phase delays of 16, 22, and 6 days were observed between GPS/GRACE, GPS/SLM, and GRACE/SLM. The data processing errors and local geophysical signals in GPS and the underestimation of land water storage in GRACE and SLM were jointly responsible for the systemic differences. The simulated data show that the misestimating of hydrological loading can explain approximately 50%, 64%, and 83% of the phase delays between GPS/GRACE, GPS/SLM, and GRACE/SLM, respectively. In addition, we obtained long-term vertical crustal motion rates by subtracting the loading deformation rates estimated by GRACE from the linear rates of the GPS. The vertical crustal motion in this region is block-dependent. The Central Yunnan block and its eastern boundary are uplifted; meanwhile, the Southwest Yunnan block, which features stretching in the horizontal direction, appears to be subsiding. The aforementioned results can provide data support for the study of water resource utilization and geodynamics in Southwest China.


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), Su and Zhan Earth, Planets and Space (2021) 73:131 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  (Luo et al. 2012;Zhang et al. 2003); I: South China block; II: Sichuan-Yunnan rhombic block; II 1 : West Sichuan block; II 2 : Central Yunnan block; III: Southwest Yunnan block 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.
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.

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:// geode sy. 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 downsampled 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 C 20 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: where def. is the elastic deformation, R is the mean radius of the Earth, and C nm and S nm 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). P nm represents degree n and order m values of the fully normalized Legendre functions. θ and are the co-latitude and longitude, respectively. h ′ n and k ′ 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: (1) where x i is the ith IMF, r m (t) is the residual signal after m decomposition and S(t) 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 GPSderived IMFs with H ≤ 0.5 into white noise signals.  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,  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.

The time series was decomposed into several IMFs
and residual terms using EEMD; 2. The Hurst parameter of each IMF was estimated; 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 . The rest of IMFs were reconstructed into the seasonal loading deformation signal. 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.

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 longterm 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.
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).

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: 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). 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  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. 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).
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.
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 S hydl , S ntal , and S ntol , respectively, in this section. Because the NTOL is faint, we . The color of the circle represents the RMS reduction, and the number in the circle is CC × 100 Su and Zhan Earth, Planets and Space (2021) 73:131 combined the signals of NTAL and NTOL and marked it as S a&o in the following content. The seasonal loading deformation in this area was the superposition of S hydl and S a&o , as shown in Fig. 8b. If S hydl 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 (A a&o ) and peak ( ϕ a&o ) of S a&o as 2.8 mm and the 172nd day, and set the peak of S hydl on the 84th day (these parameters were estimated by GRACE and SLM, as shown in "Results"). Then, we set the amplitude of S hydl to 5.2 mm (to simulate SLM), 6.5 mm (to simulate GRACE), and 9.4 mm (to simulate GPS), and combined them with S a&o to obtain the superimposed signal (S total ), as shown in Fig. 8c. The results show that the phase delays between GPS simu / GRACE simu , GPS simu /SLM simu , and GRACE simu /SLM simu 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 S hydl . 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). 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 Fig. 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 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 II 2 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.  Southwest Yunnan block. c Horizontal velocity field provided by Wang and Shen (2020)