Analyzing the Seasonal Vertical Displacement Fluctuations Using the Global Navigation Satellite System and Hydrological Load: A Case Study of the Western Yunnan Region

Abstract

The non-tectonic deformation caused by hydrological loads is an important influencing factor in GNSS vertical displacement. Limited by the temporal and spatial resolution of global models and model errors, the hydrological load results calculated by traditional methods are difficult to meet the high temporal and spatial resolution requirements of small to medium-scale regions. This paper introduces the idea of the remove–restore method, assimilates regional high-resolution hydrological data, and obtains higher temporal and spatial-resolution hydrological load results. Subsequently, utilizing data from 12 CORS observed in the western Yunnan region between January 2018 and December 2020, the quantitative relationship and variation characteristics between GNSS vertical displacement and hydrological load displacement were analyzed in detail. Furthermore, the annual signals of both were extracted using the SSA method for comparative analysis. After removing the effects of atmospheric load and non-tidal ocean load, the average correlation coefficient between GNSS vertical displacement and hydrological load displacement is 0.84, with an average reduction of WRMS (%) reaching 37.17%. The average correlation coefficient of the annual signals between GNSS vertical displacement and hydrological load deformation is 0.94, with an average reduction of WRMS (%) reaching 46.5%, indicating that the contribution of hydrological load to the GNSS non-tectonic vertical displacement annual signal is close to 50%. The research results provide scientific support and important references for studying surface tectonic deformation by removing non-tectonic deformations such as hydrological loads from GNSS vertical displacement. Additionally, it helps to explore the mechanisms of interaction between water storage migration and surface deformation.

1. Introduction

Surface mass changes can induce noticeable load responses, causing vertical surface deformation and changes in the Earth’s gravity field [1]. Studying surface load deformation can provide insights into changes in atmospheric, oceanic, and terrestrial water surface mass and their interactions, which are of great significance for global climate change and geodynamics research [2]. Analyzing methods for load-induced vertical deformation caused by atmospheric, hydrological, and non-tidal oceanic surface mass loads mainly include global surface load models and GRACE (Gravity Recovery and Climate Experiment) models [3]. Both global surface load models and GRACE models, limited by their spatial resolution, primarily represent information at spatial and long-wave scales and are unable to reveal the local characteristics of regional load-induced vertical deformation. Therefore, the current research hotspot in the field of geodesy is how to integrate high-resolution and high-precision regional models to enhance the applicability of load-induced vertical deformation monitoring methods at small to medium scales [4,5,6].

GNSS (Global Navigation Satellite System) technology provides an effective means for real-time continuous monitoring of surface deformation, allowing the reflection of vertical surface deformation through the time series of coordinates obtained from continuous observations at monitoring sites [7]. Previous studies have shown that after subtracting long-term changes caused by tectonic movements, the annual variations in the GNSS vertical displacement time series are mainly caused by surface mass migration [8]. Therefore, to ensure the reliability of load-induced vertical deformation monitoring results, most studies validate them using GNSS displacements, primarily comparing the annual variation characteristics of the two sets of results. Studies on atmospheric and non-tidal oceanic loads have shown that non-tidal oceanic loads provide better correction effects on GNSS vertical displacement in coastal areas [9,10] and studied the seasonal variations in vertical displacement caused by atmospheric loads, which can reach 18–20 mm [11,12]. In studies on hydrological loads, the deformation caused by hydrological load obtained from surface load models has been found to have varying effects on GNSS vertical seasonal changes, ranging from millimeters to centimeters [13,14,15]. Several scholars have conducted comparative studies between hydrological load displacement obtained from GRACE model data and seasonal changes in GNSS vertical displacement. The majority of the results indicate a good overall correlation and consistency between the two [16,17,18,19,20,21]. The results show that the maximum annual amplitudes of GNSS continuous stations induced by atmospheric, hydrological, and non-tidal oceanic loads can reach approximately 4 mm, 7–8 mm, and 2–3 mm, respectively. These loads can explain about 40% of GNSS vertical seasonal changes [22]. In summary, it is evident that different surface mass loads (atmospheric, hydrological, and non-tidal oceanic) have varying effects on GNSS vertical seasonal changes in different regions.

The western Yunnan region (23.5° N–27° N, 97.5° E–101.5° E), located on the southeastern side of the Qinghai–Tibet Plateau, is a significant geological research area composed of multiple blocks. It contains several large and complex fault zones, making it one of the most seismically active regions in mainland China [23]. Some scholars have studied the relationship between GNSS vertical displacement and hydrological load deformation in Yunnan Province using GRACE models and GNSS data, suggesting that hydrological loads are one of the main factors causing seasonal changes in GNSS vertical motion in the Yunnan region [24]. However, GRACE can effectively distinguish hydrological load changes within a range of approximately 300 km, but it cannot effectively discern the impact of hydrological loads on GNSS continuous station local small-scale ranges [25]. Therefore, this study uses a remove–restore method combined with regional high-resolution models to calculate the high resolution of hydrological load impacts in the western Yunnan region, aiming to reveal the characteristics of hydrological load migration on crustal vertical deformation at small to medium scales. This paper provides data processing methods for monitoring surface vertical deformation using hydrological loads and GNSS. It mainly quantitatively compares hydrological load displacement with GNSS-monitored vertical displacement and analyzes the seasonal and annual variation characteristics of both, discussing the impact of hydrological load on GNSS vertical displacement in western Yunnan.

2. Date

2.1. Data of CORS Network

The paper utilized data from 12 CORSs (Continuous Operating Reference Stations) in western Yunnan from January 2018 to December 2020. The spatial distribution of each station is indicated by the red circles in Figure 1 [26]. Additionally, to obtain high-precision CORS coordinate change time series, 15 IGS stations were selected, and their data were processed alongside CORS data. Data processing was performed using the GAMIT/GLOBK software (version 10.75). For the issues of ambiguity, ionospheric delay, and tropospheric delay, the LAMBDA method, LC_AUTCLN, and ‘Saastamoinen+GPT2w+estimation’ processing modes were respectively adopted [27]. Daily GNSS observation data errors for each station were corrected, and the GAMIT software (version 10.75) was used to obtain single-day regional loose solutions for the stations. Network adjustments were conducted using GLOBK software (version 10.75) to obtain the station coordinate changes under the ITRF (International Terrestrial Reference Frame) 2014 framework [28]. The GNSS calculations in this paper have already removed the effects of solid tides, ocean tides, and atmospheric tides. However, non-tidal loads caused by changes in surface mass, such as atmospheric, hydrological, and non-tidal oceanic, have not been removed.

2.2. Atmospheric Pressure Data

The atmospheric pressure data used in this study were obtained from the reanalyzed ERA-interim surface pressure product data from the ECMWF, with a spatial resolution of 0.25° × 0.25° [29]. The data were processed to obtain weekly values, covering the period from 3 January 2018 to 30 December 2020. Using the date of each Wednesday to represent the week, for example, 3 January 2018 represents the first week of 2018.

2.3. SLA Data

Daily global sea level anomalies data were sourced from AVISO, with a spatial resolution of 0.25° × 0.25°. This dataset integrates altimetry measurements from multiple satellites, such as TOPEX/Poseidon, Jason-1/2, and Envisat, and undergoes necessary geophysical corrections [30,31]. The data were also processed to obtain weekly values, covering the period from 3 January 2018 to 30 December 2020.

2.4. Hydrological Model

(1)

The GLDAS Global Model

To obtain soil moisture variations, this study utilized the GLDAS V2.1 NOAH model provided by NASA, with a spatial resolution of 0.25° × 0.25° [32]. This model considers soil moisture content from 0 to 200 cm depth, as well as surface water in vegetation canopies and snow water. The data have a temporal resolution of every three hours.

(2)

The CLDAS regional Model

The regional high-resolution atmospheric pressure data used in this study were obtained from the soil moisture analysis product provided by CLDAS (the China Meteorological Administration Land Data Assimilation System) V2.0. This dataset covers the Asian region (0°–65° N, 60°–160° E) and provides soil moisture analysis products with a spatial resolution of 0.0625° × 0.0625° latitude-longitude grid, divided into five layers vertically (0–5 cm, 0–10 cm, 10–40 cm, 40–100 cm, 100–200 cm). The real-time product has a lag of 1.5 h, and the near-real-time product has a lag of 2 days. This dataset was developed using various observational data sources, such as ground-based and satellite observations. It exhibits superior quality in China, with higher spatial and temporal resolutions.

3. Methods

(1)

The method of spherical harmonic approximation.

The changes in surface mass can be represented by the variation in EWH (equivalent water height). According to the theory of Earth’s load deformation, the normalized spherical harmonic expansion of the EWH variation is [33,34]:

$$\Delta h_w(\phi ,\lambda )=a{\displaystyle {\sum }_{l=1}^L{\displaystyle {\sum }_{m=0}^l\left[\Delta C_{lm}^q\cos m\lambda +\Delta S_{lm}^q\sin m\lambda \right]}}{\overline{P}}_{lm}(\sin \phi )$$

(1)

where $(\phi ,\lambda )$ represents the geocentric latitude and longitude, while $\Delta C_{lm}^q,\Delta S_{lm}^q$ is the load sphere harmonic coefficient with degree $l$ and order $m$. ${\overline{P}}_{lm}(\sin \phi )$ is the associated Legendre function with degree $l$ and $m$ order the changes in vertical displacement caused by surface mass loads can be calculated using the following equation [22,33]:

$$H(\phi ,\lambda ,t)=3\frac{{\rho }_w}{{\rho }_e}\frac{GM}{\gamma R}{\displaystyle {\sum }_{l=2}^L\frac{h_l^{\prime }}{2l+1}}{\displaystyle {\sum }_{m=0}^l\left[\Delta C_{lm}^q\cos m\lambda +\Delta S_{lm}^q\sin m\lambda \right]}{\overline{P}}_{lm}(\sin \phi )$$

(2)

where ${\rho }_e\approx 5.5\times {10}^3 \mathrm{kg}\cdot {\mathrm{m}}^{-3}$; $G$ is the gravitational constant; and $\gamma$ is the average ground gravity.

(2)

The method of loading Green’s functions.

The GNSS vertical displacement caused by surface mass loads should satisfy the dynamic equation of load deformation:

$$H(\phi ,\lambda ,t)={\displaystyle {\int }_0^{2\pi }d{\lambda }^{\prime }{\displaystyle {\int }_0^{\pi }{\rho }_w\Delta h_{w}G(\psi )R^2}}\sin {\lambda }^{\prime }d{\phi }^{\prime }$$

(3)

where $t$ represents time; $({\phi }^{\prime },{\lambda }^{\prime })$ represents the load point to be integrated on the ground; ${\rho }_w\approx {10}^3 \mathrm{kg}\cdot {\mathrm{m}}^{-3}$; $\psi$ represents the spherical angular distance between the calculated and load points; and $R$ represents the average earth radius.

The radial load Green’s function $G(\psi )$ is as follows [34]:

$$G(\psi )=\frac{R{h}_{\infty }^{\prime }}{2M\sin (\psi /2)}+\frac{a}{M}{\displaystyle \sum _{n=0}^N(h_l^{\prime }-h_{\infty }^{\prime })}{P}_n(\cos \psi )$$

(4)

where $h_l^{\prime }$ is the radial load LOVE number [35], $M$ is the mass of the Earth, and $P_n$ is the Legendre function.

(3)

Calculation of regional hydrological load based on the remove–restore method.

In gravity field research, the characteristic of the Earth’s gravity field is that the long-wave component predominates (greater than 90%), while the perturbation mass generated by terrain and crustal disturbances contributes relatively small amounts of medium and short-wave components, with the short-wave effect particularly minimal, varying from meters, decimeters, and centimeters in large, medium, and small mountainous areas, respectively. Global gravity field models mainly contain information on long-wave components. In refining regional gravity field models, the remove–restore method, combined with locally measured gravity data, can introduce information on medium- and short-wave components, thus enhancing the accuracy and resolution of regional gravity field models. Similar to this principle, this paper, based on the remove–restore method, combines global hydrological models with regional high-resolution data to calculate high-precision, high-resolution regional hydrological load results, thereby enhancing its applicability in medium- and small-scale regions.

As shown in Figure 2, the process involves the following steps: ① Download and process GLDAS global model data, Using Formula (1), the global hydrological model is expressed in terms of spherical harmonic coefficients. ② Utilize Formula (2) to estimate the reference EWH and hydrological load reference values, serving as contributions from the far field in the global background field. ③ Acquire high-resolution hydrological data for the study area from CLDAS. ④ Densify the grids of reference EWH and local high-resolution EWH into the same 1′ × 1′ grid. ⑤ Subtract the reference EWH grid from the local high-resolution EWH grid to obtain the residual EWH grid. ⑥ Utilize Formula (3) to calculate the hydrological load caused by residual EWH, introducing finer near-field load source contributions. ⑦ Sum the load impacts of residual EWH with the hydrological load reference values, obtaining a regional hydrological load with finer spectrum information at a resolution of 1′ × 1′.

From the perspective of the spectrum, the impact of the load can be divided into “near-zone impact” and “far-zone impact”. The remove–restore method not only preserves the long-wave information of the global model, that is, the far-zone impact, but also introduces the mid-short-wave information of the regional model, which is a refined expression of the near-zone impact. This method not only utilizes regional high-resolution data to improve the accuracy and spatiotemporal resolution of hydrological load results but also effectively suppresses truncation errors generated during spherical harmonic expansion.

4. Results

4.1. The Impact of Hydrological Load on Vertical Displacement in the Western Yunnan Region

In this study, the global hydrological model provided by GLDAS and the regional high-resolution hydrological data provided by CLDAS were averaged on a weekly basis to obtain results with a temporal resolution of one week. To standardize the baseline, the average of the 52 weekly values for the year 2018 was subtracted as the reference, resulting in a grid of weekly hydrological changes.

Figure 3 shows the model values of vertical deformation in the western Yunnan region calculated based on the GODAS global model, representing the reference vertical deformation. The grid results for each week are represented by the Wednesday of each week. From the figure, it can be observed that vertical deformation in the western Yunnan region exhibits significant seasonal characteristics, with the impact of hydrological load on vertical deformation ranging from −8 to 8 mm, reaching the centimeter scale. The results obtained based on the GODAS global model primarily reflect spatial and long-wave scale information. However, due to limitations in the resolution of the global model and truncation errors in spherical harmonic expansion, the short-wave information is not accurately captured. As a result, the spatial difference of the hydrological load deformation field is not pronounced, and it fails to reveal local characteristics in the region.

Figure 4 depicts the refined hydrological load impact obtained after assimilating regional high-resolution data from CLDAS, representing the impact of residual hydrological load on vertical displacement. The magnitude of vertical deformation induced by residual hydrological load ranges from −2 to 2 mm. This result does not include global background field information but provides a refined representation of near-field load sources, thus better revealing local spatial characteristics. By adding the reference hydrological load impact (Figure 3) to the residual hydrological load impact (Figure 4), we obtain the vertical deformation field of high-resolution regional hydrological load calculated based on the remove–restore method. This method assimilates regional high-resolution hydrological model data, making the spectral information of the calculation results more refined. It also effectively suppresses truncation errors in the traditional spherical harmonic approximation method, thereby enhancing the applicability of hydrological load results in medium- and small-scale regions.

4.2. The Comparison with the Vertical Displacement Seasonal Fluctuations from Hydrological Load and the GNSS

Figure 5 shows the results of the vertical displacement time series for the four selected CORS. The red dots represent the GNSS vertical displacement time series calculated by GAMIT/GLOBK software, indicating the changes in the geodetic height at CORS, after gross error detection and step processing. The blue lines represent the linear trend component. Furthermore, this study utilized the least squares method to perform low-frequency reconstruction of the nonlinear variation in GNSS vertical displacement, which better expresses its seasonal characteristics while suppressing the influence of high-frequency noise in the original time series, as depicted by the black curves. The vertical displacement at most stations ranges from −25 to 25 mm.

Figure 6 shows the vertical displacement at the XIAG station caused by hydrological load, atmospheric load, and non-tidal ocean load, where the hydrological load is based on the remove–restore method. All three types of surface mass loads exhibit significant seasonal characteristics. The impacts of hydrological load and atmospheric load are relatively large, both reaching the centimeter scale. The influence of non-tidal ocean loads is smaller, within the range of −2 to 2 mm.

To quantitatively compare and analyze hydrological load with GNSS vertical displacement, this study employed RMS (%), R two metrics, as defined in Equations (5) and (6):

$$WRMS (\%)=\frac{WRM{S}_{GNSS}-WRM{S}_{GNSS-load}}{WRM{S}_{GNSS}}$$

(5)

where $WRM{S}_{GNSS}$ is the WRMS of the GNSS vertical displacement. WRMS (%) can reflect the hydrological load influences on the GNSS’s vertical displacement.

The Pearson correlation coefficient (R) is calculated as follows:

$$R=\frac{Cov(X,Y)}{\sqrt{Var(X),Var(Y)}}$$

(6)

where $X=(x_1,x_2,\cdots ,x_N)$ represents the GNSS vertical displacement;$Y=(y_1,y_2,\cdots ,y_N)$ represents the vertical deformation caused by hydrological loads. The R values range from −1 to 1, indicating strong negative and positive correlations between the periodic phases of the two series.

To better compare and analyze the seasonal variation relationship between GNSS vertical displacement and hydrological load displacement, we performed a load correction, i.e., the atmospheric and non-tidal ocean load deformation was removed from the LSF-reconstructed GNSS vertical displacement. Figure 7 shows the comparison between the GNSS vertical displacement and the vertical displacement caused by hydrological load at XIAG, YNLJ, YNSD, and YNLC stations. The red curve represents the GNSS vertical displacement without load correction, while the blue curve represents the GNSS vertical displacement after load correction. It can be observed from the figure that the phases and amplitudes of all CORS are relatively consistent. These CORS are located in areas of significant soil moisture variation, indicating that hydrological load deformation can effectively explain the seasonal variation of GNSS vertical displacement in the western Yunnan region.

It is worth noting that in the western Yunnan region, the influence of non-tidal ocean load on atmospheric load has a compensatory effect (as shown in Figure 6), so the amplitude difference of GNSS vertical displacement before and after load correction is not significant. The specific effects can be seen in the statistical confidence of

Table 1. Correlation coefficient and contribution rate WRMS (%) between vertical time series from GNSS and hydrological load vertical displacement in the western Yunnan region from 2018 to 2020.

CORS	Without Load Correction		Atmospheric Load and Non-Tide Ocean Load Correction
	Correlation Coefficient	WRMS (%)	Correlation Coefficient	WRMS (%)
XIAG	0.68	25.33	0.78	33.60
YNCX	0.84	31.18	0.85	33.59
YNGM	0.79	23.67	0.84	34.86
YNJD	0.71	25.78	0.80	32.48
YNLC	0.72	26.63	0.82	35.96
YNLJ	0.79	34.80	0.86	41.99
YNRL	0.77	33.74	0.86	43.28
YNSD	0.75	29.87	0.85	38.38
YNTC	0.81	33.56	0.88	39.13
YNYA	0.78	30.98	0.86	37.13
YNYL	0.77	27.73	0.85	33.53
YNYS	0.78	31.82	0.86	38.75

. The correlation coefficient calculation results of the GNSS vertical displacements and hydrological load displacements for 12 CORS are shown in Table 1. When the atmospheric and non-tidal ocean load effects are not deducted from GNSS vertical displacements, the average correlation is 0.77.the smallest correlation being at the YNGM, with a value of 0.68. And the largest correlation is at the YNCX, with a value of 0.84. After correcting for atmospheric and non-tidal ocean load effects, the average correlation between the two becomes 0.84, with the smallest correlation still at the YNGM, with a value of 0.78, and the largest correlation being at the YNTC, with a value of 0.88. This indicates that most CORS have a strong correlation between GNSS vertical displacement and hydrological load displacement, and after correcting for atmospheric and non-tidal ocean load effects on GNSS vertical displacement, the correlation between the two can be significantly improved, verifying the reliability of the methods and theories presented in this paper.

To further illustrate the consistency of seasonal variations between GNSS vertical displacements and hydrological load deformations, this paper quantitatively evaluates whether hydrological load deformations can effectively correct non-tectonic deformations in GNSS vertical displacement by using the WRMS (%) obtained from subtracting hydrological load displacement from GNSS vertical displacement. If the WRMS (%) is positive, it means that hydrological load displacement can effectively correct non-tectonic deformations in GNSS vertical displacement. From Table 1, it can be seen that the values for all CORS are positive, indicating that hydrological load displacement can effectively remove non-tectonic deformations for Yunnan CORS. When the atmospheric and non-tidal ocean load effects are not deducted from GNSS vertical displacement, the average WRMS (%) for all CORS is 29.59%, which increases to 37.17% after load correction. This indicates that the contribution of hydrological loads to the seasonal vertical displacements of CORS stations in western Yunnan is above 30%. When studying the impact of hydrological loads using GNSS vertical displacement, it is effective and necessary to deduct the effects of atmospheric loads, non-tidal ocean loads, and other loads.

4.3. The Comparison with the Vertical Displacement Annual Variation from the GNSS and Hydrological Load

This study employs the SSA method to separately extract the annual signals of GNSS vertical displacement and hydrological load deformation [36,37,38], as shown in Figure 8. The red curve represents the annual signal of GNSS vertical displacement, while the blue curve represents the annual signal of hydrological load. The statistical results are presented in

Table 2. Statistical values of the vertical annual variation from the GNSS and hydrological load.

Station	GNSS		Hydrological Load
	Annual Amplitude /mm	Annual Phase /°	Annual Amplitude /mm	Annual Phase /°
XIAG	5.71 ± 0.24	19.30 ± 5.26	3.75 ± 0.16	23.46 ± 5.47
YNCX	5.52 ± 0.29	31.18 ± 5.70	3.57 ± 0.17	22.93 ± 6.98
YNGM	5.89 ± 0.48	28.38 ± 4.64	4.25 ± 0.20	29.25 ± 5.76
YNJD	7.71 ± 0.34	26.25 ± 5.64	3.85 ± 0.19	22.20 ± 6.09
YNLC	6.61 ± 0.44	15.87 ± 6.47	4.09 ± 0.20	21.55 ± 5.78
YNLJ	6.19 ± 0.22	33.08 ± 4.47	3.70 ± 0.18	24.29 ± 5.76
YNRL	7.52 ± 0.30	35.55 ± 5.16	4.63 ± 0.20	22.88 ± 5.29
YNSD	7.51 ± 0.26	21.34 ± 4.37	3.99 ± 0.20	23.51 ± 6.04
YNTC	8.10 ± 0.27	32.69 ± 3.78	4.36 ± 0.19	23.12 ± 5.38
YNYA	7.18 ± 0.23	34.71 ± 4.10	3.55 ± 0.16	23.74 ± 5.67
YNYL	8.01 ± 0.33	20.96 ± 5.14	3.71 ± 0.18	25.16 ± 5.83
YNYS	6.61 ± 0.17	34.11 ± 3.34	3.52 ± 0.16	25.29 ± 5.53

. Both signals exhibit strong consistency in amplitude and period. The annual signal of GNSS vertical displacement fluctuates within ±5.5 mm to ±8 mm, while the annual signal of hydrological load fluctuates within ±3.5 mm to ±5 mm. The difference in annual phase values ranges from 0.87° to 12.67°. There is significant variation in the amplitude of vertical displacement among different CORS, indicating that GNSS signals are more sensitive to vertical displacement changes and can effectively reveal local characteristics of the region.

As shown in

Table 3. R and WRMS (%) between vertical annual variation from the GNSS and hydrological load in the western Yunnan region from 2018 to 2020.

Station	Correlation Coefficient	WRMS (%)	Station	Correlation Coefficient	WRMS (%)
XIAG	0.91	51	YNRL	0.93	51
YNCX	0.97	45	YNSD	0.94	47
YNGM	0.92	43	YNTC	0.95	43
YNJD	0.94	44	YNYA	0.96	46
YNLC	0.92	44	YNYL	0.94	42
YNLJ	0.96	54	YNYS	0.95	48

, the average correlation coefficient between vertical annual variation from the GNSS and hydrological load is 0.94. The lowest correlation coefficient is observed at the XIAG station, which is 0.91, while the highest correlation coefficient is observed at the YNCX station, which is 0.97. There is a strong correlation between GNSS vertical annual variation and hydrological load vertical annual variation. The average WRMS (%) is 46.5%, indicating that hydrological loads contribute to nearly 50% of the non-tectonic annual signals of GNSS vertical displacements.

5. Discussion

In traditional hydrological load impact studies, global hydrological models, or GRACE data, are generally used. The spatial resolution of global hydrological models is typically 0.5° × 0.5° or 0.25° × 0.25°, while the GRACE gravity satellite has a resolution of about 300 km [39,40]. Both have a temporal resolution on the monthly scale. However, neither global hydrological models nor GRACE data can meet the demand for high precision and high resolution in medium- and small-scale regional areas. This article proposes using the remove–restore method to assimilate high-resolution hydrological data in regional areas, thereby improving the applicability of hydrological load impact in medium- and small-scale regions. The remove–restore method not only preserves the long-wave information of the global model, that is, the far-zone impact, but also introduces the mid-short-wave information of the regional model, which is a refined expression of the near-zone impact. Experimental trials were conducted using the western Yunnan region as an example, obtaining hydrological load results with a resolution of 1′ × 1′, thus demonstrating the feasibility of this method. Furthermore, GNSS data were utilized for comparative analysis, validating the reliability of this method. From the perspective of temporal resolution, the hydrological load impact in this study yields results on a weekly scale, representing some improvement. However, the temporal resolution of regional hydrological data is 3 h, and the weekly averaging method does not fully exploit its advantage of high temporal resolution. In future research, the applicability of daily hydrological load results will be further investigated.

The findings of this study indicate that monitoring hydrological load displacement using GNSS is feasible, and GNSS signals are more sensitive and accurate, capable of fully revealing the local characteristics of the region [19,41]. The strong consistency between the hydrological load displacement monitored by GNSS and regional hydrological data is evident from the correlation coefficient results. However, there are still some differences in amplitude between the two. The annual signal of GNSS vertical displacement fluctuates within ±5.5 mm to ±8 mm, while the annual signal of hydrological load fluctuates within ±3.5 mm to ±5 mm. The differences are not only related to model errors in regional hydrological data but also to the processing of GNSS vertical displacement time series. GNSS time series not only include non-tectonic motions such as hydrological load, atmospheric load, and non-tidal ocean load but also incorporate other environmental changes and geodynamic variations, such as earthquakes, volcanoes, and tectonic motions. In particular, the influence of atmospheric load and non-tidal ocean load is significant. When studying hydrological load deformation using GNSS, it is necessary to remove the influence of atmospheric load and non-tidal ocean load. After load correction, there is a noticeable improvement in the correlation between GNSS vertical deformation and hydrological load, demonstrating the feasibility and necessity of the technical approach proposed in this paper.

In addition, thermal expansion effects can induce changes in the vertical displacement of CORS stations. Jiang found that the thermal expansion effect is more pronounced at GNSS stations located in mid- to high latitudes, with a maximum vertical displacement amplitude of 1.8 mm, while the average amplitude at low-latitude stations is only 0.16 mm [42]. This effect is one of the important factors contributing to the differences in annual amplitudes between these locations, but it was overlooked in this study. In future research, it is crucial to fully consider all influencing factors in GNSS signals to accurately capture hydrological load displacement in GNSS vertical displacement.

6. Conclusions

To enhance the applicability of hydrological load results in small to medium-scale regions, this study introduced the concept of the remove–restore method, assimilating regional high-resolution hydrological data to obtain higher spatiotemporal-resolution hydrological load results. Subsequently, this study utilized GNSS vertical displacement from 12 CORS and hydrological load displacement spanning the time period from 2018 to 2020 to investigate the seasonal variations in vertical motion in the western Yunnan region. This study compared and analyzed the quantitative relationship and characteristic variations between GNSS vertical displacement and hydrological load displacement. Additionally, the SSA method was employed to extract the annual variation signals of GNSS vertical displacement and hydrological load displacement, demonstrating strong consistency between the two. The results are as follows:

(1)

The hydrological load displacement calculated in this study based on the remove–restore method has a higher spatiotemporal resolution. This provides scientific support and important references for future research aiming to remove non-tectonic deformations, such as hydrological loads, from GNSS vertical displacement and study surface tectonic deformations.

(2)

The seasonal motion trends of GNSS vertical displacement and hydrological load displacement are consistent. However, the displacement values of hydrological load are generally smaller than those of GNSS, indicating that hydrological load displacement can explain a portion of the seasonal variations in GNSS vertical motion. After removing the effects of atmospheric loads and non-tidal ocean loads from GNSS signals, the average correlation coefficient between the two increased from 0.77 to 0.84, and the average WRMS (%) increased from 29.59% to 37.17%. This suggests that in the western Yunnan region, approximately 30% or more of the non-tectonic deformations in GNSS vertical displacement originate from hydrological load.

(3)

The average correlation coefficient between the annual variations of GNSS and hydrological load reaches 0.94, indicating a strong correlation. The average WRMS (%) is 46.5%, suggesting that hydrological loads contribute to nearly 50% of the non-tectonic annual variations of GNSS vertical displacement.

The research results can effectively reveal the finer deformation characteristics of hydrological loads in the region and help to explore the mechanism of action between water storage migration and surface displacement.