*Article* **Recovering Regional Groundwater Storage Anomalies by Combining GNSS and Surface Mass Load Data: A Case Study in Western Yunnan**

**Pengfei Xu 1,2, Tao Jiang 2,\*, Chuanyin Zhang 2, Ke Shi <sup>3</sup> and Wanqiu Li <sup>4</sup>**


**Abstract:** The redistribution of surface mass (e.g., atmosphere, soil water, oceans, and groundwater) can cause load responses, resulting in vertical deformations of the crust. Indeed, the global navigation satellite system (GNSS)-based continuously operating reference stations (CORS) are able to accurately measure the vertical deformation caused by surface mass loads. In this study, the CORS was used to invert groundwater storage anomalies (GWSA), represented by the equivalent water height (EWH), after removing the effect of the non-groundwater surface mass load (atmospheric, groundwater, and non-tidal oceanic loads) from the vertical deformation monitored by CORS. In addition, the global and regional high-resolution surface mass models were combined to calculate the high-precision load deformation field in in western Yunnan using the remove–restore method, thereby obtaining more accurate surface mass load data and improving the accuracy of the inverted GWSA results. In order to assess the feasibility of the CORS inversion for the GWSA used, 66 CORS stations in western Yunnan Province were considered, presenting weekly GWSA data from 10 January 2018 to 31 December 2020. The results revealed significant seasonal variation in GWSA in the study area, showing an amplitude range of −200–200 mm. This approach is based on the already-established CORS network without requiring additional set-up costs. In addition, the reliability of CORS inverse results was assessed using Gravity Recovery and Climate Experiment (GRACE) inverse results and actual groundwater monitoring data. According to the obtained results, GWSA can be monitored by both CORS and GRACE data; however, CORS provided a more effective spatiotemporal resolution of GWSA. Therefore, the CORS network combined with surface mass load data is able to effectively monitor the spatiotemporal dynamics of GWSA in small-scale areas and provides important references for the study of hydrology.

**Keywords:** ground water storage; surface mass load; groundwater monitoring station; GNSS; remove-restore method

#### **1. Introduction**

As a densely populated country, China is relatively poor in water resources per capita [1]. Moreover, China is characterized by an uneven spatiotemporal distribution of water resources due to the impacts of climatic and topographic characteristics. The GWSA (groundwater storage anomalies) aims to estimate groundwater storage, thereby maintaining the ecological balance of groundwater systems. The overexploitation of groundwater resources has seriously restrained sustainable development in many regions [2,3]. Therefore, the assessment of the spatiotemporal distribution of regional GWS (groundwater storage) is of great scientific importance for studying water circulation and groundwater allocation as well as for preventing groundwater droughts [4,5].

**Citation:** Xu, P.; Jiang, T.; Zhang, C.; Shi, K.; Li, W. Recovering Regional Groundwater Storage Anomalies by Combining GNSS and Surface Mass Load Data: A Case Study in Western Yunnan. *Remote Sens.* **2022**, *14*, 4032. https://doi.org/10.3390/rs14164032

Academic Editors: Gino Dardanelli and Shuanggen Jin

Received: 27 June 2022 Accepted: 16 August 2022 Published: 18 August 2022

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

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

At present, the two main GWS monitoring approaches consist of in situ monitoring of groundwater level fluctuations and using Gravity Recovery and Climate Experiment (GRACE) and Global Land Data Assimilation System (GLDAS) data to invert the spatiotemporal characteristics of GWS [6,7]. In fact, consecutive groundwater level data, with relatively high accuracies, can be obtained at monitoring well stations using the first approach. However, this approach is costly since it requires several monitoring stations with high density to ensure accurate regional monitoring. Moreover, most well stations can be located in plains regions, making it challenging to monitor groundwater storage in mountainous areas. Meanwhile, due to the complex regional geological structures, some parameters (e.g., specific yield) are difficult to determine, making it difficult to effectively assess GWSA using discretely distributed well stations. On the other hand, the GRACE and GLDAS data-based approach is often applied to large-scale regional GWSA monitoring [8,9]. Although the quality of gravity satellite data is relatively high, there is considerable interference upon high-order spherical harmonic coefficients, resulting in low spatial resolutions. Furthermore, these satellite data have been lacking since 2011, with a data gap of nearly one year between GRACE/GRACE Follow-on (GRACE-FO), explaining the inability of GRACE data to provide high-resolution continuous monitoring the GWSA data in small-scale regions [10,11]. Therefore, given the complex physical conditions, high cost of monitoring groundwater storage through well stations, and the low-resolution GRACE monitoring data, the continuously operating reference system (CORS) inversion method for GWSA in small-scale areas requires further study and discussion.

The CORS is a ground observation system based on the global navigation satellite system (GNSS) observation system. These data are derived from long consecutive monitoring data of satellite navigation signals, providing real-time and periodic data through communication facilities. The CORS can monitor the dynamics of regional geodetic heights and their consecutive long-time series [12,13]. Indeed, several researchers have demonstrated the ability of the CORS network to obtain real-time information from selected monitoring stations in certain regions, providing references for monitoring spatial dynamics as well as comprehensive continuous regional observation data [14–16]. In addition, the loaddeformation theory demonstrated that changes in the surface environmental mass (e.g., atmosphere, surface water, groundwater, and ocean) lead to load vertical deformation, thereby influencing the seasonal periodic signals in the geodesic height variations of CORS stations [17–19]. Argus [20] used the seasonal signals of GNSS vertical data to invert the terrestrial water storage anomalies (TWSA) in California, showing consistent spatial distribution data with that inverted using GRACE. In addition, He [21] inverted the changes in the TWSA in Yunnan Province from 2010 to 2014 and compared the GNSS data with those of GRACE and GLDAS, discussing the possibility of GNSS to separately operate and monitor the TWSA. Several studies have revealed consistent findings with the hypothesis that the seasonal vertical displacement of GNSS is strongly related to the load-deformation caused by TWSA and have assessed the reliability of GNSS data to invert the TWSA [22–28]. However, studies on groundwater inversion have only assessed the correlation between seasonal GNSS and GWSA data due to the complexity of seasonal signals of GNSS vertical displacements, while only a few studies have inverted GWSA data using the CORS.

To address this challenge, this study aims to verify the feasibility and reliability of the combined CORS high-resolution surface mass load inversion for GWSA based on an already-established CORS network, thus avoiding additional establishment costs. Moreover, the density of CORS stations and their continuous and immediate high-precision observation make it possible to perform high-resolution GWSA monitoring in a smallscale region. Yunnan Province in China has experienced a severe water shortage and continuous drought. Indeed, the GWSAs have directly affected the regional economy and local ecological environment in Yunnan Province. Therefore, assessing the spatiotemporal distribution characteristics of the GWSA can provide important reference significance for human production activities as well as for relevant decision-makers in the regions. In order to improve the precision and stability of the CORS inversion, the global and regional

high-resolution surface mass models were combined in this study to calculate the highprecision load deformation field in Yunnan Province based on the remove-restore method. Indeed, researchers have eliminated three surface mass loads (atmospheric, soil water, and non-tidal oceanic loads) from the non-linear time series of CORS geodetic heights through data processing, and then the CORS data were used to invert the GWSAs based on load deformation theory and inversion models [29,30]. To assess the applicability of the method, researchers have used data from 66 CORS stations in Western Yunnan Province and high-resolution surface mass load results and performed a weekly GWSA grid in the study area from 10 January 2018 to 31 December 2020 to assess the effectiveness and reliability of CORS inverse results using GRACE inverse results and groundwater monitoring station data.

#### **2. Materials and Methods**

#### *2.1. Data Used*

#### 2.1.1. CORS Network Data

As is shown in Figure 1, the study region is circled by red lines. The data were collected from 66 CORS stations located in the study region, covering the 2018–2020 period. In addition, 15 international GNSS service (IGS) stations were selected in this study to obtain high-precision CORS coordinate time series and to resolve the solution of CORS data. The longitude and latitude coordinates of the IGS stations are reported in Table 1. GAMIT and GLOBK were also used in this study to process the collected data [31]. On the other hand, by correcting the daily errors of GNSS data for each station, a solution was achieved in some zones of the study region through GAMIT (Table 2). GLOBK was used for a global adjustment of the network and to present the time series of station coordinates in the international terrestrial reference frame (ITRF) [32]. As is shown in Table 2, the CORS data solution removed the influences of the solid, sea, and atmospheric tides, while the influences of non-tidal elements (e.g., atmosphere pressure, soil water, and sea level) were maintained.

**Figure 1.** Geographic location of the study region.


**Table 1.** The 15 IGS sites used for the calculation of the CORS network in Western Yunnan Province.

**Table 2.** The main parameters used in the GAMIT calculation.


#### 2.1.2. Atmospheric Pressure Data

In order to calculate the effect of atmospheric load, the global atmospheric pressure data were derived in this study from ECMWF's (European Centre for Medium-Range Weather Forecasts) re-analysis of the 0.25◦ × 0.25◦ ERA-interim surface pressure product data (https://www.ecmwf.int/, accessed on 15 May 2021), covering the 3 January 2018–30 December 2020 period [33]. The original data were averaged by week, which is reported by Wednesday of each week. Similar to the solutions of non-tidal oceanic and soil water loadings, the weekly grid time series of atmospheric loading were obtained using Equations (1) and (2).

On the other hand, the regional high-resolution atmospheric pressure data were downloaded from the CLDAS-V2.0 data (CLDAS Atmospheric Driving Field) product of the China Meteorological Administration Land Data Assimilation System (CLDAS) from the China Meteorological Data Network (http://data.cma.cn/, accessed on 15 May 2021). The product data consist of temperature, pressure, humidity, wind speed, precipitation, and other meteorological parameter data. The hourly product data cover the Asian region (0◦~65◦N, 60◦~160◦E), with a spatial resolution of 0.0625◦ × 0.0625◦. The dataset includes data observed at over 2400 and 40,000 national and regional automatic meteorological stations, respectively, for operational assessment. By combining global and regional atmospheric pressure data, a high-resolution atmospheric loading deformation field of the study area was determined using the remove–restore method. The soil water loading was calculated using the same method.

#### 2.1.3. Sea Level Anomaly Data

In order to calculate the effect of Non-tidal ocean load, the daily global sea-level anomaly (SLA) data were derived from the 0.25◦ × 0.25◦ archiving validation and interpretation of satellite oceanography (AVISO) data (https://www.aviso.altimetry.fr/, accessed on 15 May 2021), covering the 2 January 2018–30 December 2020 period. The AVISO data

integrate not only multi-satellite-derived sea-level data (e.g., TOPEX/Poseidon, Jason-1/2, and Envisat) but also the associated geophysical correction, including tidal and inverse barometer corrections, thereby representing the effects of non-tidal ocean loading [34,35].

#### 2.1.4. Soil Water Data

Calculation of soil water load requires soil water data. In this study, the V2.1 GLDAS/ Noah model, provided by the National Aeronautics and Space Administration (NASA) and National Centers for Environmental Prediction (NCEP), was used to obtain the soil water data, with a spatial resolution of 0.25◦ × 0.25◦ (https://mirador.gsfc.nasa.gov, accessed on 15 May 2021) [36]. The GLDAS model takes into account soil water, canopy water, and snow water (ranging from 0 to 200 cm), while groundwater is excluded. The temporal resolutions of the data are 3-hourly and monthly. In this study, the 3-hourly resolution GLDAS data, from 3 January 2018 to 30 December 2020, were first downloaded, then averaged by week to obtain weekly soil water data. These data were used to determine the influence of soil water loading, whereas the monthly soil water data covered the January 2018–December 2020 period. Indeed, the monthly GLDAS data combined with the monthly GRACE data can be used to determine GWSA as well as to assess the validity and feasibility of the CORS inversion for GWSA.

On the other hand, the regional high-resolution surface water height data were downloaded in this study from the CLDAS data product of the China Meteorological Data Network, which consists of hourly soil moisture data covering the Asian region (0◦–65◦N, 60◦–160◦E), with a spatial resolution of 0.0625◦ × 0.0625◦.

#### *2.2. GRACE-FO Mascon Solutions*

To compare and analyze the CORS and GRACE inverse results for GWSA, GRACE-FO RL06 mascon solutions, provided by the Center for Space Research (CSR) with a spatial resolution of 0.25◦ × 0.25◦, were used to reflect the monthly changes in terrestrial water storage (http://www2.csr.utexas.edu/grace/RL06\_mascons.html, accessed on 15 May 2021) [37]. The data covered the October 2018–December 2020 period. The collected data were corrected for degree-1, C20 (degree 2 order 0), C30 (degree 3 order 0), and GIA to reflect the dynamics of terrestrial water storage (Table 3). It should be noted that in comparison with GRACE, the GRACE-FO data processing added the C30 correction into its procedure. Therefore, the regional GWSA could be calculated through the combination of monthly data of terrestrial water storage and the monthly data of soil water provided by GLDAS.

**Table 3.** The main corrections used in Mascon solutions.


#### *2.3. Groundwater Monitoring Station Data*

The daily data observed at 15 groundwater monitoring stations (Figure 1) from 26 September 2018 to 25 November 2020 were used in this study to assess the accuracy of the CORS inverse results for the GWSA. The collected data were averaged by week. The applicability and reliability of the integrated solution method can be assessed by comparing the observed data with those of the CORS inversion for the GWSA.

#### *2.4. Precipitation Data from Weather Stations*

To analyze the influences of precipitation on the spatiotemporal distribution of GWSA in the study region, daily precipitation data observed at eight weather stations in Yunnan Province from (Figure 1) 1 January 2018 to 31 December 2020 were used in this study (http://data.cma.cn, accessed on 15 May 2021).

#### **3. Method**

#### *3.1. High-Resolution Surface Mass Load Based on the Remove–Restore Method*

The spherical harmonic analysis method was expressed using Equation (1), reflecting the medium- and long-wave components at global scales, while Equation (2) is the load Green's function method, which mainly reflects the short-wave components at small scales. Since equivalent water height (EWH) can reflect the changes in the mass of the surface environment (e.g., atmosphere, soil water, oceans, and groundwater) based on the load theory of spherical harmonic coefficients, the changes in geodetic heights caused by loads can be calculated using the following equation [17,38]:

$$H(\boldsymbol{\varrho}, \boldsymbol{\lambda}, t) = 3 \frac{\rho\_{\rm w}}{\rho\_{\rm \varepsilon}} \frac{GM}{\gamma R} \sum\_{l=2}^{L} \frac{h\_{l}^{\prime}}{2l+1} \sum\_{m=0}^{l} \left[ \Delta C\_{lm}^{\boldsymbol{\varrho}} \cos m\boldsymbol{\lambda} + \Delta S\_{lm}^{\boldsymbol{\varrho}} \sin m\boldsymbol{\lambda} \right] \overline{\mathcal{P}}\_{lm}(\sin \boldsymbol{\varrho}) \tag{1}$$

where *<sup>ρ</sup><sup>e</sup>* ≈ 5.5 × 103 kg·m−<sup>3</sup> denotes the average density of the solid earth; *<sup>G</sup>* is the gravitational constant; and *γ* is the average gravity of the ground.

On the other hand, the changes in the geodetic heights of CORS stations resulting from the changes in surface mass load can be expressed using the load Green's function equation [17,38]:

$$H(\varphi,\lambda,t) = \int\_0^{2\pi} d\lambda' \int\_0^{\pi} \rho\_W \Delta h\_W G(\psi) a^2 \sin \lambda' d\varphi' \tag{2}$$

where *H*(*ϕ*, *λ*, *t*) denotes the geodetic height; *t* represents time; (*ϕ*, *λ*) denotes the longitude and latitude of the calculated point; and (*ϕ* , *λ* ) denotes the load point on the ground. The density of water was considered as *<sup>ρ</sup><sup>w</sup>* ≈ <sup>10</sup><sup>3</sup> kg·m<sup>−</sup>3; *<sup>ψ</sup>* denotes the spherical angular distance between the calculated and load points; *G*(*ψ*) denotes the green function of the radial loads; and *R* denotes the average earth radius.

The green function of the radial loads can be expressed as follows [17,38]:

$$G(\psi) = \frac{Rh\_{\infty}'}{2M\sin(\psi/2)} + \frac{a}{M} \sum\_{n=0}^{N} \left(h\_l' - h\_{\infty}'\right) P\_n(\cos\psi) \tag{3}$$

where *h <sup>l</sup>* is radial load Love number; *M* denotes the earth mass; and *Pn* denotes the Legendre function.

In this study, the global and regional high-resolution models were combined to calculate the high-precision load deformation field in the study area using the remove–restore method. Figure 2 shows the procedure used in this study [12,19].

First, the global EWH model was used to estimate the reference EWH and load vertical deformation in the study area using the spherical harmonic analysis method, thus obtaining the spatial long-wave component information. Second, the high-resolution regional EWH model was encrypted into a 1 × 1 grid before removing the reference EWH grid from the regional high-resolution EWH grid to obtain the residual EWH using the remove method. Third, the load impact of the residual EWH change on the geodetic height was estimated using the load Green's function method, introducing the short wave components of the study area. Finally, the load impacts obtained in the first and third steps were added using the restore method. By using these steps, the high-resolution and high-precision load deformation field of the study area was obtained. This approach not only introduces the short-wave component of the study area and improves the applicability of the model in the small-scale study area but also effectively removes the truncation error generated during the spherical harmonic expansion and reduces the calculation uncertainties.

**Figure 2.** Flowchart of the high-precision load deformation field solution.

#### *3.2. GWSA Inversion Using the Combined CORS and Surface Mass Load Method*

EWH was used in this study to reflect the mass changes in the surface environment (atmospheric pressure, sea level, soil water, and groundwater levels) [17]. Based on load deformation theory, the normalized expansion formula of the spherical harmonic load was determined using the load-deformation theory, according to the following equations [17,38,39]:

$$
\Delta h\_{\text{w}}(\varphi,\lambda) = R \sum\_{d=1}^{L} \sum\_{m=0}^{l} \left[ \Delta C\_{lm}^{q} \cos m\lambda + \Delta S\_{lm}^{q} \sin m\lambda \right] \overline{P}\_{lm}(\sin q) \tag{4}
$$

where (*ϕ*, *λ*) denotes the geocentric latitude and longitude of the ground calculation point; Δ*C<sup>q</sup> lm*, <sup>Δ</sup>*S<sup>q</sup> lm* are the load sphere harmonic coefficients, with *l* degree and *m* order; and *Plm*(sin *ϕ*) denotes the associated Legendre function, with *l* degree and *m* order.

Different thicknesses and radius disks of the same mass were placed on the surface. Figure 3 shows the related load deformation and its relationship with distance, indicating clear load responses near fields. The load-deformation was 1/10 of that of the center at a distance from the center of the disk equal to three times its radius, while no loaddeformation was observed at a distance from the center of disk equation to 10 times its radius. Therefore, the vertical displacement and load-deformation were mainly influenced by the load points within a limited distance range. The use of GNSS in the inversion of the regional GWS showed that stations on the edges of the study area tended to be simultaneously influenced by the inside and outside loads of the study area. By neglecting the influence of the outside mass changes, the outside loads can be restrained to the border region, thus leading to false inversion on the edges of the study area, explaining the larger inverted region than the study area. The expansion should consider the sensitive range of load responses. Therefore, the 66 CORS stations covered the study region and its peripheral areas.

**Figure 3.** Vertical deformation caused by disks of different masses.

The inversion model of Argus [20] used in this paper:

$$((Ax - b)/\sigma)^2 + \beta^2 (L(x))^2 \to \min \tag{5}$$

where *σ* denotes the standard deviation of the observation value vector of geodetic height changes; *A* denotes the coefficient matrix of the green function, which can be calculated using Equation (3); *x* is the EWH of the corresponding grid point; *b* denotes the observation value of the GNSS geodetic height changes of CORS stations; *β* denotes the smoothing factor; and *L* denotes the Laplace operator.

Due to the limited number of GPS observation stations in the study area, the number of equations is smaller than the number of unknowns, resulting in a rank loss of the coefficient matrix of the normal equation. The inversion of the equivalent water height changes using GNSS observation data is, therefore, an ill-posed equation problem. To address this challenge, the ridge estimation reported by Hoerl [40] was used in this study, which is a classical regularized approach, to determine the solution of *β* in Equation (5). Indeed, the L-curve method was used in previous studies to select the best-regularized parameters [41].

The influence of surface mass changes was calculated using the load-deformation theory in Equation (1). Besides the impact of groundwater changes, atmospheric pressure, sea level, and soil water changes were considered in the present study, while other less influential factors (e.g., rivers, lakes, and reservoirs) were neglected. Based on the weekly geodetic height values, the non-linear time series of each CORS station were obtained. After removing the influences of the non-groundwater data, the residual time series were obtained, including load and non-load vertical deformations of groundwater. Afterward, Equations (2) and (5) were used to invert GWSA from the residual time series. The method used is illustrated in Figure 4. In this study, the integral radius of the load Green's function and the smoothing factor (*β*) were set to 2◦ and 0.01, respectively.

**Figure 4.** Flowchart of CORS inversion for GWSA.

#### *3.3. GWSA Inversion Using GRACE Product Data*

The terrestrial water storages (*TWS*) are the summation of soi l water, groundwater, surface runoff, snow, and canopy water. Contrary to soil water and groundwater, the contributions of surface runoff and canopy water changes to TWS are minor, and thus they can be neglected. Indeed, TWS in the study region is mainly affected by soil water, snow water, and groundwater [10]. Δ*GWS* was calculated in this study using the following formula:

$$
\Delta GWS = \Delta TWS - \left(\Delta SM + \Delta SWE\right) \tag{6}
$$

where Δ*GWS* is the groundwater storage; Δ*TWS* denotes the changes in TWS obtained from GRACE data; Δ*SM* and Δ*SWE* denote the changes in soil water and snow water, respectively, obtained from GLDAS product data.

#### **4. Results**

#### *4.1. GWSA Inversion Results Using CORS Data*

The time series of non-linear geodetic height changes were first processed using gross error detection and linear item removal, as shown by the red points in Figure 5. The black curve in Figure 5 shows the time series of geodetic heights after the periodic fast Fourier transformation (FFT) reconstruction. According to the magnitude of the power spectral density, this paper uses the first eight periodic signals to reconstruct [42,43]. Indeed, the periodic FFT reconstruction is a low-pass filtering method used to decrease or inhibit the high-frequency noise of the time series, thus improving the stability of the inversion results. The geodetic heights of CORS stations corresponded to the temporal resolution of surface mass loads. In order to improve the accuracy of the CORS inversion, the influences of high-resolution surface mass loads, including atmospheric, soil water, and non-tidal oceanic loads, were removed from the reconstructed time series of non-linear geodetic heights, as shown by the purple, blue and green curves in Figure 5. The residual time series obtained after the removal process, according to Equations (2) and (5), can therefore be used to invert GWSA within the coverage of the CORS network, which is reflected by EWH.

**Figure 5.** Processed geodetic height changes of different CORS stations and surface mass loads.

Figure 6 shows the annual differences between the maximum and minimum values of the reconstructed geodetic height time series at the 66 CORS stations, ranging from 15 to 35 mm. According to Figure 5, the atmospheric and soil water loads revealed the highest effects among the three types of loads. The results revealed that the vertical deformation, caused by the soil water loads, ranged from −8 to 8 mm, while non-tidal oceanic load showed the lowest effect, with a vertical deformation range of −2–2 mm.

**Figure 6.** Range of the reconstructed geodetic heights at each CORS station in the study area.

Surface mass loads, namely atmospheric, non-tidal oceanic, and soil water loads, showed different correlations with the geodetic height time series of CORS stations, demonstrating the influences of these surface mass loads. Therefore, in order to compare the obtained results, the influences of loads were first removed, then the WRMS ratio of the GNSS time series was calculated according to the following formula [44]:

$$WRMS(\%) = \frac{WRMS\_{GNSS} - WRMS\_{GNSS-load}}{WRMS\_{GNSS}} \tag{7}$$

where *WRMSGNSS* is the WRMS of the reconstructed GNSS geodetic height time series. Positive and negative values of *WRMS*(%) indicate a decrease and increase in the WRMS of the GNSS time series, respectively. It should be noted that the absolute value of *WRMS*(%) can reflect the load influences on the GNSS's non-linear geodetic height time series.

On the other hand, the Pearson correlation coefficient (*R*) was computed in this study according to the following equation. It is widely used to measure the degree of correlation between two variables [45,46]:

$$R = \frac{\mathbb{C}ov(X, Y)}{\sqrt{Var(X)\_{\prime}Var(Y)}}\tag{8}$$

where *X* = (*x*1, *x*2, ··· , *xN*) denotes the reconstructed time series of the observed GNSS geodetic heights; *Y* = (*y*1, *y*2, ··· , *yN*) denotes the time series of the vertical deformation caused by surface mass loads. The R values range from −1 to 1, indicating strong negative and positive correlations, respectively, between the periodic phases of the two series.

In total, 12 stations were selected to assess the influences of the three surface mass loads on the reconstructed time series of CORS geodetic heights using *R* and *WRMS* values (Table 4). Unlike the non-tidal oceanic load, the *R* and *WRMS* values of atmospheric and soil water loads were all positive. According to the obtained results, the *R* and *WRMS* values of the atmospheric load ranged from 0.50 and 0.66 and 7.61 to 12.80%, respectively, whereas the *R* and *WRMS* values of the soil water load ranged from 0.64 to 0.81 and 20.54 to 34.80%, respectively. On the other hand, the non-tidal oceanic load was negatively correlated with the reconstructed time series of the CORS geodetic heights, which is consistent with the reported by Munekane and Nordman [47,48]. However, the absolute *WRMS* value of the non-tidal oceanic load ranged from 3.24 to 6.39%, suggesting a low influence on the CORS geodetic heights. This result was due to the offsetting effects of the atmospheric and soil water loads, making the influence of the non-tidal oceanic load difficult to reflect in the time series of the CORS geodetic heights. In other words, atmospheric and soil water loads exhibited stronger influences than that of the non-tidal oceanic load. The *R* and *WRMS* values of total loads ranged from 0.81 to 0.89 and 33.84 and 43.93%, respectively (Table 4). On the other hand, by removing the surface mass loads, decreases in the *WRMSGNSS*−*load* values were observed. This finding demonstrates not only the reliability of the reconstructed CORS geodetic height time series and surface mass load deformations but also the effectiveness of the integrated solving process used in the present study.


**Table 4.** Correlation coefficient and WRMS values between different loads and the reconstructed geodetic height time series.

In order to ensure the reliability of the inversion and reduce the influence of the lowquality CORS data, the quality of each CORS dataset was first evaluated before inversion, and then different weights were assigned to CORS stations. The initial weight of every station was set at 1, while the weight of the low-quality stations was decreased based on the entire time series data and their RMS values.

It was demonstrated earlier that the integration solution of the CORS data revealed a weekly GWSA grid over the 3 January 2018–30 December 2020 period was established, with the GWSA expressed in EWH. As is shown in Figure 7, the GWSA showed significant spatial and seasonal variations. The results revealed decreases in GWS from February to June each year, which might be due to the significant seasonal and spatial decreases in groundwater recharge (e.g., rainfall infiltration). From August to the next January, however, increases in GWS were observed in the study area from August to January due to the increase in rainfall amounts. In addition, the inverted GWS using the CORS network showed a relatively high spatial resolution (Figure 7).

**Figure 7.** Weekly GWSA over the January 2018-December 2020 period.

The errors generated from the GNSS data processing, the surface mass model, and the uncertainty in the CORS inversion model may affect the accuracy of the GWSA inversion. The inverse distance weight interpolation method was used in this study to map the time series of the GWSA of CORS stations. Table 5 shows the statistical results of the interpolated GWSA time series of 12 CORS stations. The maximum and minimum values of the inverted GWSA were −200 and 200 mm, respectively. In addition, there were some differences in the

amplitude variation of GWSA at different CORS stations, indicating that the CORS inverse results had a high resolution. The mean and standard deviation (SD) values of GWSA at all CORS stations ranged from 3.16 to 15.03 mm and 55.87 to 110.39 mm, respectively. The reliability of the CORS inverse results was assessed in the subsequent part using the GRACE inverse results and groundwater monitoring station data.


**Table 5.** Statistics of CORS inverse results for GWSA (mm).

In order to further analyze the temporal features of GWSA based on the sequence data of the GWSA grids, the study area was first divided into four sub-study areas, namely the eastern (BCHU, SYUN, and NJIA), northern (BAIS, YNYL, and JCHU), western (SUDI, YINJ, and TOBG), and southern (CHA3, MENT, and YNSD) sub-study areas, and then four CORS stations were selected randomly as examples. In Figure 8, the curves show the time series of GWSA at the corresponding stations, while the bar charts indicate the precipitation amounts in certain sub-study areas observed at nearby weather stations (Figure 1). In comparison with Figure 7, Figure 8 more directly reflects the temporal features of GWSA. The GWSA dynamics in different sub-study areas were relatively similar. Indeed, since groundwater recharge is derived mainly from precipitation in the study region, GWSA was strongly correlated with precipitation and exhibited significant seasonal differences. According to the obtained results, GWS exhibited main wave crest shapes. From February and March of each year, GWS showed the lowest values from February/March to June/July, then increased significantly from July to October due to the significant increase in the precipitation amounts, followed by a decrease in GWS in December and January. Therefore, precipitation was the main influencing factor on the GWSA in the study region.

In order to quantitatively analyze the groundwater drought in the study area within 3 years, the groundwater severity index (DSI) [49] was used in this study. This index was calculated using the following formula:

$$CORS-DSI\_{i,j} = \frac{GWSA\_{i,j} - \overline{GWSA\_j}}{\sigma\_j} \tag{9}$$

where *i* denotes the year from 2018 to 2020; *j* represents the month from January to December; *CORS* − *DSIi*,*<sup>j</sup>* represents the groundwater drought index DSI; *GWSAj* and *σ<sup>j</sup>* represent the mean and standard deviation of GWSA, respectively. The classification results of groundwater drought, derived from DSI, are reported in Table 6. As reported above, the monthly GWSA data were obtained from the monthly average of the weekly GWSA time series.

**Figure 8.** GWSA at Different CORS stations and precipitation amounts in different sub-study areas.

Figure 9 shows the temporal variation of CORS-DSI in the study area from 2018 to 2020. The results showed a significant seasonal variation in CORS-DSI. In addition, the CORS-DSI values in the eastern and northern sub-regions of the study area were above 0.8 in most months, while only a few months showed CORS-DSI values slightly below 0.8, indicating a mild groundwater drought (Figure 9a,b). On the other hand, the western and southern subregions of the study area showed several months with CORS-DSI values below −0.8, while the southern sub-region exhibited relatively significant groundwater drought, with CORS-DSI close to −1.3 over the November 2018–February 2019 and April 2020–August 2020 periods (Figure 9c,d). It is worth noting that the groundwater drought analysis requires more than 10 years of observation data, while only 3 years of CORS station data were used in this study to invert the GWSA. Therefore, to obtain more detailed and accurate drought analysis results, it is suggested to consider long-term observation data in future groundwater drought analyses.

**Figure 9.** CORS-DSI in different sub-study areas.


**Table 6.** Correlation coefficients and WRMS (%) between different mass surface loads and the reconstructed geodetic height time series.

#### *4.2. Comparison with GRACE Data*

In order to assess the reliability of the CORS inverse results for GWSA, the GRACE inverse results for GWSA and the observed groundwater level data were used in this study. EWH is considered in the GRACE data to reflect the TWS. GWSA can be computed by removing soil water storage provided by the GLDAS data product. In this study, the period between October 2018 and December 2020 was considered due to the loss of data between GRACE and GRACE-FO. The monthly GRACE data were compared in this study with the monthly inverted CORS-based GWSA data during the year 2019. The obtained results are shown in Figure 10.

**Figure 10.** Comparison between CORS and GRACE inverse results.

Feng [50] showed a delay in the GRACE-based land water storage changes by 2 months compared to precipitation changes. Therefore, a 2-month phase delay was performed for the GRACE time series. As can be seen from Figure 10, relatively consistent GRACE- and CORS-based GWSA trends were obtained following the phase delay correction.

In total, 15 CORS were selected randomly as examples. The values in Table 7 indicate the correlation coefficients of the CORS- and GRACE-based GWSA results. Due to the missing data in 2018, the correlation coefficients were computed in this study based on the 2019 data. Compared with the inversion CORS-based data, the GRACE data had a 2-month phase delay. After phase delay correction, the coefficients were significantly improved. Except for FGON, LJGC, and YNRL, showing correlation coefficients of 0.68, 0.69, and 0.65, respectively, the correlation coefficients of all other stations were above 0.7. The highest correlation coefficient was 0.85 in YNJD, indicating a strong correlation. The computed correlation coefficients showed a spatial variation due to the low spatial resolution of the inversion GRACE-based on data. As is shown in Figure 10, the interpolation method used in this study resulted in relatively consistent trends of GRACE monitoring results at different sites without showing significant spatial differences. In addition, the CORS-based

GWSA results were based on long-term monitoring GNSS data. Indeed, the density of the CORS stations in the study region allowed us to obtain a higher spatial resolution (Figure 7), showing clear differences in the GWSA time series that are also clearer (Figure 10).


**Table 7.** Correlation coefficients between CORS and GRACE inverse results.

Figure 11 shows the annual variation in CORS- and GRACE-based GWSA from 2018 to 2020. The spatial distribution characteristics of the two monitoring results were slightly consistent, showing relatively large annual changes in GWSA in the eastern and southern parts of the study area. The main spatial differences between the two monitoring results were observed in the south-central part of the study area. Overall, the amplitude of the annual variation in CORS-based GWSA was greater than that of GRACE-based GWSA. This finding may be due to the different monitoring methods of the two data products. Indeed, GRACE measures the integrated regional effect of mass redistribution. Its low spatial resolution makes it difficult to comprehensively reflect the effective information of the small-scale region, whereas the CORS data product is based on GNSS measurement methods, which provides real-time geometric deformation information at different locations in the study area, thereby fully reflecting the signal changes in the region.

**Figure 11.** Annual variation in CORS- and GRACE-based GWSA results.

In summary, both inverse results reflected the changes in GWSA. The results indicate that the GNSS-based CORS data are more sensitive to GWSA, showing obvious local spatial characteristics compared to those of the GRACE data.

#### *4.3. Comparison with Groundwater Monitoring Data*

In order to further assess the reliability of CORS-based GWSA results, the inverted GWSA data were interpolated and compared to groundwater levels observed at groundwater monitoring stations. The daily groundwater level data were collected from September 2018 to December 2020. These data were processed for gross error detection and weekly averaging. In order to quantitatively compare the groundwater level data with GWSA, it is necessary to convert the groundwater level data to equivalent water height by considering the specific yield of the area [10]. The specific yield value for the study area was set to 0.05 [51,52].

The blue curves in Figure 12 indicate the CORS-based GWSA results observed in four stations, while the red line indicates the result of multiplying the groundwater level by the specific yield value. Both results are the average values of data observed over the year 2019. The results revealed consistent trends of CORS-based GWSA and groundwater levels, showing the same order of amplitude changes. On the other hand, the correlation coefficients between the two data ranged from 0.62 to 0.82, indicating strong positive correlations (Table 8). Therefore, the results demonstrated that the combined CORS and high-resolution surface mass load data effectively invert GWSA in small-scale areas, providing accurate seasonal trends for GWSA.

**Figure 12.** Comparison between CORS-based GWSA and groundwater monitoring results.

**Table 8.** Correlation coefficients between CORS-based GWSA and observed GWSA at groundwater monitoring stations.


#### **5. Discussion**

Considering the difficult physical conditions and high cost of constructing groundwater monitoring stations, as well as the low-resolution problem of the GRACE inversion results, the present study aims to assess the feasibility and reliability of CORS for GWSA inversion. Numerous studies have analyzed the correlation of TWSA or GWSA with GNSS vertical displacement [24,26,27]. In addition, He [21] inverted TWSA using GNSS vertical displacement. However, few studies have directly inverted GWSA using the CORS data product. Therefore, to improve the accuracy and stability of the inversion results, global and regional high-resolution surface mass models were combined to determine the high-precision load deformation field in western Yunnan using the remove–restore method. In addition, the reliability of CORS-based GWSA results was assessed using the observed groundwater levels. Although the obtained results showed accurate GWSA and demonstrated the validity and reliability of the methods used in this study, there are still some challenges to address. Data covering a period longer than 2.5 years were observed only at 15 stations among the 66 CORS stations, thereby affecting the accuracy of the inverse results for some months. Since the data processing is complicated, the study period considered in this study was only 3 years, resulting in non-obvious drought characteristics and trends. Therefore, this study focused mainly on the seasonal variation in GWSA in the

study area. In addition, the specific yield in this study was set to 0.05 for the entire study area. It should be noted that the specific yield can vary significantly depending on the depth of groundwater wells and the lithological characteristic of the aquifers. Therefore, different specific yield values need to be provided for different groundwater wells to obtain more accurate GWSA results using groundwater monitoring data.

The correlation coefficients between GRACE- and CORS-based GWSA were above 0.65 following the 2-month delay phase correction. The variation in the EWH of GRACEand CORS-based GWSA showed a consistent magnitude of variation. In addition, the annual variation in CORS-based GWSA results was slightly larger than that of GRACEbased GWSA results. Although the principles of GRACE and CORS methods used in the GWSA inversion are different, the results suggested that both methods can effectively reflect the seasonal trend of GWSA in the study area. Indeed, GRACE product data are based on remote sensing techniques to comprehensively reflect the regional effects with redistributed mass, resulting in a low spatial resolution and difficulties in effectively transmitting small-scale regional information. Moreover, the filtering process used in the GRACE data processing might affect the real signal. On the other hand, the CORS data product is derived from GNSS-based monitoring data, allowing us to obtain in real time the geodetic height changes at a certain location, thus assuring faster response and better performance in reflecting the inter-region signal dynamics. Both methods can reflect the GWSA trends. However, the signal intensity and GRACE trend monitoring results at different locations can be basically the same in small-scale areas, while the differences in CORS network monitoring results are more obvious. GRACE is ineffective in determining features of small-scale study areas due to its limited spatial resolution, while CORS network monitoring results exhibit higher spatiotemporal resolution, with continuous and highprecision observation data, making the CORS method advantageous. Moreover, the CORS method can capture local area signals that are difficult to monitor by remote sensing techniques (e.g., GRACE). In summary, the CORS-based GWSA results revealed higher spatiotemporal resolution than that of GRACE-based GWSA results and are more sensitive to groundwater storage changes.

#### **6. Conclusions**

The present study aims to invert GWSA in western Yunnan using CORS and highprecision surface mass load data. This approach is based on the already-established CORS network without requiring extra construction costs. This approach is, indeed, able to provide a high-resolution GWSA grid sequence within the coverage of the CORS network due to the high density of CORS stations and consecutive high-precision GNSS monitoring data. Western Yunnan Province was considered in the present study to assess the validity of the approach used. The results demonstrated that the approach is able to independently and effectively invert a high-resolution GWSA grid within the CORS network. In addition, the spatiotemporal distribution characteristics of GWSA were analyzed in this study using the precipitation data observed at meteorological stations in western Yunnan. Moreover, to test the reliability of the obtained results, the GRACE and groundwater monitoring data were used in the present study. The main conclusions are as follows:

1. The correlation coefficients between the CORS geodetic height time series and the vertical deformation of the surface mass loads were all above 0.8, indicating strong positive correlations. In addition, the percentage of WRMS decreased from 33.84 to 43.93% following the load removal, demonstrating the effectiveness of the data processing used in the present study and the feasibility of CORS inversion for GWSA. The vertical deformations caused by surface mass loads contributed significantly to the seasonal signals of CORS geodetic height changes. Among the three surface mass loads, atmospheric and soil water loads were more influential, with an amplitude ranging from −8 to 8 mm, while the non-tidal oceanic load showed the lowest influence with an amplitude range of −2–2 mm.


**Author Contributions:** P.X., T.J. and C.Z. conceived and designed the experiments; P.X. performed the experiments, analyzed the data, and wrote the paper; T.J., W.L. and C.Z. helped in the discussion and revision. K.S. provided the data from the CORS network. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (No. 42074020), the Open Fund Project of the Key Laboratory of Marine Environmental Survey Technology and Application of the Ministry of Natural Resources (MESTA-2020-A001), and the Comprehensive CORS network and multi-source data monitoring of surface stability changes and pilot analysis of geological disaster precursor capture (AR2114).

**Data Availability Statement:** The raw/processed CORS network and groundwater level data are available from the author upon reasonable request. Correspondence and requests for data should be addressed to T.J.

**Acknowledgments:** Thanks to the CORS network data provided by the Yunnan Technical Center of Base Surveying and Mapping and the data from groundwater monitoring stations provided by the China Geological Environment Monitoring Institute.

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

#### **References**

