Estimation of Terrestrial Water Storage Variations in Sichuan-Yunnan Region from GPS Observations Using Independent Component Analysis

: GPS can be used to measure land motions induced by mass loading variations on the Earth’s surface. This paper presents an independent component analysis (ICA)-based inversion method that uses vertical GPS coordinate time series to estimate the change of terrestrial water storage (TWS) in the Sichuan-Yunnan region in China. The ICA method was applied to extract the hydrological deformation signals from the vertical coordinate time series of GPS stations in the Sichuan-Yunnan region from the Crustal Movement Observation Network of China (CMONC). These vertical deformation signals were then inverted to TWS variations. Comparative experiments were conducted based on Gravity Recovery and Climate Experiment (GRACE) data and a hydrological model for validation. The results demonstrate that the TWS changes estimated from GPS(ICA) deformations are highly correlated with the water variations derived from the GRACE data and hydrological model in Sichuan-Yunnan region. The TWS variations are overestimated by the vertical GPS observations the northwestern Sichuan-Yunnan region. The anomalies are likely caused by inaccurate atmospheric loading correction models or residual tropospheric errors in the region with high topographic variability and can be reduced by ICA preprocessing.


Introduction
Terrestrial water storage (TWS) variations play an important role in the study of the global water cycle and water resource management [1,2]. Modern space-geodetic satellites provide an effective and efficient approach to monitor and quantify TWS variations. For instance, Gravity Recovery and Climate Experiment (GRACE) data have been widely used to measure the variations of large-scale water storage [3][4][5]. However, GRACE-derived TWS values have low spatial resolution (300-500 km) and low temporal resolution (monthly). Another popular method is the use of hydrological assimilation models, such as Global Land Data Assimilation Systems (GLDAS) [6], to understand the water cycle process and effect of land surface processes on climate [5]. However, GLDAS models neglect either reservoir water or groundwater and lack comprehensive in-situ measurements.
TWS variations can deform the Earth's surfaces, especially in the vertical direction. These elastic responses can be precisely monitored by GPS measurements. Some studies have used continuous GPS measurements to investigate hydrological deformations on the Earth's surface in regions with significant water variations, and the results are in good agreement with the GRACE-derived displacements [7][8][9][10]. These studies indicate that GPS is a useful real-time tool to estimate TWS variations. Chew and Small [11] used vertical GPS observations to assess the TWS anomalies in the High Plains (United States) caused by the 2012 drought. In that study, the vertical GPS time series were assumed to be linearly correlated with TWS variations, but an explanation of their theoretical relationships was not provided. Argus et al. [12] proposed that the TWS inversion method uses vertical GPS time series of land motions based on Green's functions [13] and successfully applied the method to monitor the TWS in California. Fu et al. [14] used GPS as an independent measuring method to estimate the TWS variations in the Washington and Oregon areas and made synthetic tests of the inversion method. Jin and Zhang [15] used continuous GPS observations to estimate the TWS variations in the southwestern U.S., and the inversion results showed good agreement with the GRACE-derived TWS and hydrological models. Zhang et al. [16] estimated the annual EWH (equivalent water height) changes of the lower Three-Rivers headwater region in southwestern China using the vertical displacements of GPS time series obtained from the Crustal Movement Observation Network of China (CMONOC). Zhong et al. [17] analyzed GPS-inverted mass changes under different density distributions of GPS stations and compared the inversion results with GRACE/GFO solutions in spatial and temporal domains.
To improve the accuracy of TWS inversion using GPS observations, Fok and Liu [18] proposed using terrain-corrected GPS and GRACE-derived "virtual GPS stations" to solve the problem of irregularly distributed GPS stations and the neglect of the terrain effect in southwestern China. Liu et al. [19] introduced external TWS-derived displacements and used them as pseudo GPS sites in the inversion process. To improve the TWS inversion results in the boundary regions, Shen et al. [20] developed a boundary-included inversion model that accounts for the mass change effect from the near but exterior region. Lai et al. [21] applied truncated singular value decomposition regularization to more stably solve the inversion problem.
Relatively few studies have discussed the reliability in obtaining hydrological deformation information from GPS observations. Most previous studies considered TWS variations to be the main cause of such elastic deformations on the Earth's surface. The existing TWS inversion methods therefore generally assume that the vertical GPS coordinate changes are identical to the hydrological deformations [12,[14][15][16][18][19][20]. However, other mass loading changes, such as atmospheric loading (ATML) [22,23] cannot be neglected. An estimation of the ATML and corresponding corrections are sometimes necessary, but accurate ATML correction models are typically difficult to obtain due to the limited atmospheric pressure data resolution and terrain change complexities [24]. Continuous GPS time series also include other (seasonal) error sources, such as non-tidal ocean loading (NTOL), higher-order ionospheric effects, bedrock thermal expansion, and errors in the antenna phase center variation models [25,26]. Ignoring these effects is conducive to accurately estimating TWS variations. Independent component analysis (ICA) has been proposed to extract the hydrological signals from regional GPS network observations [27][28][29][30]. Previous studies [31][32][33] have confirmed that ICA is an effective method to extract hydrological deformations from continuous GPS time series.
In this paper, we propose to extract hydrological deformations from regional GPS networks and use the independent GPS displacements to estimate the TWS variations (hereafter referred to as GPS(ICA)). A case study is carried out in the Sichuan-Yunnan region, and the results are compared with those from the GRACE and GLDAS models.

Continuous GPS Observation and Inversion Method
The study area is located in the Sichuan-Yunnan region of China, as shown in Figure 1. We selected 47 GPS sites from CMONOC to obtain their six-year vertical coordinate time series (2011)(2012)(2013)(2014)(2015)(2016). GAMIT/GLOBK software was used to generate the daily solutions in ITRF2008 reference frame, with main parameters including the station coordinates, zenith tropospheric delays (one per hour), and gradients (one per day in the N/S and E/W directions). Tropospheric delays were corrected using the Vienna Mapping Function HGMF and DGMF. The ANTEX model was used to correct the antenna phase center variations. Non-tidal surface loading effects and high-order ionospheric delays are not corrected in the time series.

Continuous GPS Observation and Inversion Method
The study area is located in the Sichuan-Yunnan region of China, as shown in Figure  1. We selected 47 GPS sites from CMONOC to obtain their six-year vertical coordinate time series (2011)(2012)(2013)(2014)(2015)(2016). GAMIT/GLOBK software was used to generate the daily solutions in ITRF2008 reference frame, with main parameters including the station coordinates, zenith tropospheric delays (one per hour), and gradients (one per day in the N/S and E/W directions). Tropospheric delays were corrected using the Vienna Mapping Function HGMF and DGMF. The ANTEX model was used to correct the antenna phase center variations. Non-tidal surface loading effects and high-order ionospheric delays are not corrected in the time series. Cycle variations in the continuous GPS time series are mainly caused by mass loading changes on the Earth's surface. The deformations caused by mass loading can be computed based on elastic theory [13]. The deformations observed by GPS can therefore be used to study the mass loading changes on the Earth's surface. The elastic deformation induced by mass loading can be calculated by: where is the angular distance between the point load and location, ∆ is the mass load, and are the radius and mass of the Earth, respectively, and ℎ and are the elastic Love numbers and Legendre polynomials, respectively. Equation (1) can be expressed as ( ) = ( ) • , where ( ) represents the Green's functions and is the surface grid mass load. In this paper, we used the code provided by Chanard et al. [34] to compute the Green's functions [35] based on the Preliminary Reference Earth Model (PREM) [36]. Furthermore, as elastic vertical motion rapidly decreases with distance from a load [14], we thus extended the studied area by 2° in both longitude and latitude when computing the Green's functions to avoid the edge effect [14].
The elastic deformation is contributed from the total mass load on the surface. Suppose that there are grids and GPS sites in the inversion region. The displacement at the th GPS site position can then be expressed as: Cycle variations in the continuous GPS time series are mainly caused by mass loading changes on the Earth's surface. The deformations caused by mass loading can be computed based on elastic theory [13]. The deformations observed by GPS can therefore be used to study the mass loading changes on the Earth's surface. The elastic deformation induced by mass loading can be calculated by: where θ is the angular distance between the point load and location, ∆M is the mass load, R and M e are the radius and mass of the Earth, respectively, and h n and P n are the elastic Love numbers and Legendre polynomials, respectively. Equation (1) can be expressed as u(θ) = G(θ)·x, where G(θ) represents the Green's functions and x is the surface grid mass load. In this paper, we used the code provided by Chanard et al. [34] to compute the Green's functions [35] based on the Preliminary Reference Earth Model (PREM) [36]. Furthermore, as elastic vertical motion rapidly decreases with distance from a load [14], we thus extended the studied area by 2 • in both longitude and latitude when computing the Green's functions to avoid the edge effect [14]. The elastic deformation u is contributed from the total mass load on the surface. Suppose that there are m grids and n GPS sites in the inversion region. The displacement at the ith GPS site position can then be expressed as: where A i is the design matrix, θ i,j is the angular distance between a grid mass in the inversion area and the ith GPS site position, and x is the grid mass load in the inversion area, because the distant mass load is neglected. The grid mass load derived displacements at the GPS sites can be expressed as: Considering that there is usually a limited number of GPS sites in the inversion, a regularization constraint is applied to obtain the TWS variations across the study area. The inversion model in Equation (3) with constraints thus becomes: where A is the design matrix consisting of Green's functions, x is the unknown mass load at each grid cell, u is the vertical displacement at each GPS site, and L is the Laplacian operator [37]. The templates of the Laplacian smoothing constraints on neighboring grids are shown in Figure 2.
where is the design matrix, , is the angular distance between a grid mass in the inversion area and the th GPS site position, and is the grid mass load in the inversion area, because the distant mass load is neglected. The grid mass load derived displacements at the GPS sites can be expressed as: Considering that there is usually a limited number of GPS sites in the inversion, a regularization constraint is applied to obtain the TWS variations across the study area. The inversion model in Equation (3) with constraints thus becomes: where is the design matrix consisting of Green's functions, is the unknown mass load at each grid cell, is the vertical displacement at each GPS site, and is the Laplacian operator [37]. The templates of the Laplacian smoothing constraints on neighboring grids are shown in Figure 2. The regularized solution to estimate the grid mass load in Equation (4) is used to minimize the combination in Equation (5): where is a roughness factor that specifies by how much the neighboring pixels values may differ, and is the weight matrix. The roughness factor balances the relative importance between data misfit ( ( − ) 2 ) and smoothness (‖( ) 2 ‖) and is selected based on the L-curve method [38] in our study. Figure 3 shows the changes of misfit and smoothness with different roughness factors. We evaluated a group of roughness factors and selected = 0.03 in our study. The regularized solution to estimate the grid mass load x in Equation (4) is used to minimize the combination in Equation (5): where λ is a roughness factor that specifies by how much the neighboring pixels values may differ, and P is the weight matrix. The roughness factor λ balances the relative importance between data misfit ( ((Ax − u)) 2 P ) and smoothness ( (Lx) 2 ) and is selected based on the L-curve method [38] in our study. Figure 3 shows the changes of misfit and smoothness with different roughness factors. We evaluated a group of roughness factors and selected λ = 0.03 in our study.  The final TWS variations can be estimated by: The covariance matrix of the final estimation can be expressed as [39]: The final TWS variations can be estimated by: The covariance matrix D λ of the final estimation can be expressed as [39]: where σ 2 0 is the variance of unit weight. In the case of a certain number of grids to be inverted, the more GPS sites number means the more observation equations, and the inversion would be with better precision and stability.

Inversion Method Based on ICA
The TWS inversion method using GPS measurements is based on the assumption that the deformations monitored by GPS are caused by TWS variations. ICA has been proven to be a useful tool for extracting geophysical signals from observed geophysical series, such as GPS time series [27][28][29][30][31][32]40], GRACE time series [41][42][43], and InSAR monitoring data [44,45]. In this paper, ICA is used to extract hydrological deformations from GPS time series, and the FastICA algorithm [46,47] is applied in the separation. Regional GPS time series are used as mixed signals and treated as ICA inputs. Some temporally independent components (ICs) can be separated from the GPS coordinate time series by Y(t) = BX(t), where X(t) represents the GPS coordinate time series and B is the separating matrix. Then the inputs can be expressed as: where the row vectors in Y are ICs estimated by ICA processing and the column vectors in the invert matrix of B are the spatial responses of ICs. Principal component analysis (PCA) was used to preprocess the inputs before we performed ICA and an appropriate number of principal components (PCs) need to remain for the following ICA processing. We used a trial and error approach to select the number of PCs [48], and found that in most cases the choice of the number of remaining PCs had little effect on the separation of the top several ICs. The contribution of each IC to the GPS observations can be characterized using ratio values defined as: where the displacement GPS ICk is the product of ICk and the corresponding spatial response. A smaller ratio value implies a more significant contribution for the corresponding IC [28], and the ICs are reordered in ascending order [40]. The top six ICs are shown in Figure 4a and the ratio values of the GPS ICs are shown in Figure 4b. We found IC1 contributes most for the GPS coordinate time series. The separated ICs were compared with external data to trace their potential geophysical causes. We computed the daily hydrological loadings at each GPS site and compared the results with GPS IC1 , as shown in Figure 5. The averaged hydrological loading displacement is consistent with GPS IC1 in terms of its amplitude and phase, and the correlation coefficient is 0.78. The term GPS IC1 is therefore thought to represent hydrological deformation, and GPS IC1 is then used in the subsequent inversion step.
where the displacement GPSICk is the product of ICk and the corresponding spatial response. A smaller ratio value implies a more significant contribution for the corresponding IC [28], and the ICs are reordered in ascending order [40]. The top six ICs are shown in Figure 4a and the ratio values of the GPSICs are shown in Figure 4b. We found IC1 contributes most for the GPS coordinate time series.  The separated ICs were compared with external data to trace their potential geophysical causes. We computed the daily hydrological loadings at each GPS site and compared the results with GPS IC1, as shown in Figure 5. The averaged hydrological loading displacement is consistent with GPS IC1 in terms of its amplitude and phase, and the correlation coefficient is 0.78. The term GPSIC1 is therefore thought to represent hydrological deformation, and GPSIC1 is then used in the subsequent inversion step.    The separated ICs were compared with external data to trace their potential geophysical causes. We computed the daily hydrological loadings at each GPS site and compared the results with GPS IC1, as shown in Figure 5. The averaged hydrological loading displacement is consistent with GPS IC1 in terms of its amplitude and phase, and the correlation coefficient is 0.78. The term GPSIC1 is therefore thought to represent hydrological deformation, and GPSIC1 is then used in the subsequent inversion step. An inversion result using mass loading-corrected displacements is also applied for comparison. The ATML and NTOL were computed using the International Mass Loading Service [49] in the center-of-figure (CF) frame using the model of Modern Era Retrospectiveanalysis for Research and Applications (MERRA) and the Max Planck Institute Ocean Model (MPIOM) 06, respectively. We averaged the ATML and NTOL series into daily results and removed them from the GPS time series. The corrected GPS time series was then used to invert the TWS variation in the Sichuan-Yunnan Region. Previous studies on the thermal effect on the vertical displacements of GPS stations in China [50,51] indicate that the effect of thermal expansion on the vertical GPS coordinate time series is limited in low-latitude areas. The thermal expansion effect is thus not considered in our case study of the Sichuan-Yunnan region.

GRACE Measurements and GLDAS Hydrological Models
The GRACE-derived TWS variations and products of the GLDAS hydrological models were used to compare with our inversion results. For the GRACE TWS data, we used the surface mass change data based on the spherical harmonics from JPL (maximum degree/order: n = 60) [52,53] (available at https://grace.jpl.nasa.gov/, accessed on 16 November 2021). For the products, the non-tidal variability was corrected using the time variable background model known as Atmosphere and Ocean De-aliasing Level-1B (AOD1B). The degree 1 coefficients were estimated from Sun et al. [54]. The C20 (degree 2 order 0) and C30 coefficients were replaced with satellite laser ranging solutions [55]. We multiplied the GRACE data by the provided scaling grid product. A detailed description of the data processing, scaling factor derivation, and caveats are available in [52]. The missing month was then interpolated from the adjacent two months. The TWS was also estimated using the sum of soil moisture in all layers, accumulated snow, and plant canopy surface water from the GLDAS model, and further used to compare with the GRACE and GPS-derived TWS variation results. We used the products of the monthly data from the GLDAS-2.1 Noah model with a spatial resolution of 1 • × 1 • .

Estimated TWS Variations in the Sichuan-Yunnan Region Using Different Methods
We estimated the monthly TWS variations of year 2011 based on the inversion theory and hydrological deformation extracted from the vertical GPS coordinate time series in Section 2.2. The results are shown in Figure 6. We also adopted the conventional approach to invert the TWS variations using GPS observations, in which the ATML and NTOL were removed. The inversion results are shown in Figure 7. The GRACE and GLDAS-derived TWS variations in the Sichuan-Yunnan are presented in Figures 8 and 9, respectively. Note that there are differences in the color bar scale in Figures 6-9.   A comparison of the monthly TWS variations obtained by the four different methods in Figures 6-9 shows that the monthly changes and spatial variations are similar overall. In general, the TWS gradually increases from northeast to southwest, and the variations in the southwest of Yunnan Province are the most apparent. However, the TWS variation results obtained by the GPS inversion exhibit abnormal changes in the northwestern region of Sichuan-Yunnan, where there is very complex terrain and some GPS sites are located above 4000 m. Additional comparisons of the TWS variations in the northwestern Sichuan-Yunnan region are analyzed and presented in the next section. We spatially averaged the TWS grid data in the Sichuan-Yunnan region to analyze the TWS change rules of the four  Figure 10a. The temporal change rules of the TWS obtained using the different methods are in good agreement; all reach a minimum in February-March (winter), and a maximum in August-October (late summer). By computing the difference between the TWS values in September and March in each grid, we obtained the distribution of the maximum TWS for each method. The results are shown in Figure 10b.     A comparison of the monthly TWS variations obtained by the four different methods in Figures 6-9 shows that the monthly changes and spatial variations are similar overall. In general, the TWS gradually increases from northeast to southwest, and the variations

Northwestern Sichuan-Yunnan Region
Another clear difference among the four inverted results is the abnormal TWS change of the GPS inverted results in northwestern Sichuan Yunnan, as shown in Figure 10. To further investigate this anomaly, we compare the TWS change in the northwestern part of the Sichuan-Yunnan region and the entire region. We average the monthly TWS results in

Northwestern Sichuan-Yunnan Region
Another clear difference among the four inverted results is the abnormal TWS change of the GPS inverted results in northwestern Sichuan Yunnan, as shown in Figure 10. To further investigate this anomaly, we compare the TWS change in the northwestern part of the Sichuan-Yunnan region and the entire region. We average the monthly TWS results in northwestern Sichuan-Yunnan (longitude: 99.5 • E-101.5 • E, latitude: 29.5 • N-31.5 • N) and the entire region for all of the investigated years. The results are shown in Figure 11, in which the legends of northwestern Sichuan-Yunnan results are marked with an asterisk.
Remote Sens. 2022, 14, x FOR PEER REVIEW 11 of 15 GPS(ICA), GRACE, and GLDAS. In contrast, the variations of GPS-inverted TWS in northwestern Sichuan-Yunnan are higher than those in the entire region from April to August for all six years. The comparisons indicate that the GPS(ICA) inversion results are more consistent with those of GRACE and GLDAS than that of the GPS (ATML and NTOL corrected) observations. To quantitatively compare the estimated TWS variations in northwestern Sichuan-Yunnan and the entire region using the different methods, referring to the fitting model for GPS coordinate time series [56], we fit the series in Figure 11 according to: ( ) = + + sin(2 ) + cos(2 ) + sin(4 ) + cos(4 ) + where − are the coefficients and can be estimated using the least squares method. We can obtain the annual amplitude for each series by = √ + . The annual amplitudes of the estimated TWS variations for the different methods in northwestern Sichuan-Yunnan and the entire region are shown in Table 1. The in Table  1 is defined as = ( − )/ . The seasonal amplitude of the GPS-derived TWS in northwestern Sichuan-Yunnan region shows significant difference, nearly 50% larger than the results of the entire region, whereas no such notable discrepancy is found in the results obtained using the other three methods.  There is no significant difference between the TWS variations in northwestern Sichuan-Yunnan and the entire Sichuan-Yunnan region for the results obtained using GPS(ICA), GRACE, and GLDAS. In contrast, the variations of GPS-inverted TWS in northwestern Sichuan-Yunnan are higher than those in the entire region from April to August for all six years. The comparisons indicate that the GPS(ICA) inversion results are more consistent with those of GRACE and GLDAS than that of the GPS (ATML and NTOL corrected) observations.
To quantitatively compare the estimated TWS variations in northwestern Sichuan-Yunnan and the entire region using the different methods, referring to the fitting model for GPS coordinate time series [56], we fit the series in Figure 11 according to: where a − f are the coefficients and can be estimated using the least squares method. We can obtain the annual amplitude A for each series by A = √ c 2 + d 2 . The annual amplitudes of the estimated TWS variations for the different methods in northwestern Sichuan-Yunnan and the entire region are shown in Table 1. The r in Table 1 is defined as r = (A 2 − A 1 )/A 1 . The seasonal amplitude of the GPS-derived TWS in northwestern Sichuan-Yunnan region shows significant difference, nearly 50% larger than the results of the entire region, whereas no such notable discrepancy is found in the results obtained using the other three methods.

Discussion
The anomalies of the GPS-derived TWS variations in northwestern Sichuan-Yunnan indicate that the coordinate time series used in the inversion contains deformation information or errors other than hydrological deformation. Considering the fact that northwestern Sichuan-Yunnan is near the Qinghai-Tibet Plateau and the terrain is extremely complex (Figure 1), the inaccurate ATML correction in this area may be the cause of anomalies in the GPS inversion result. The height of the GPS sites in northwestern Sichuan-Yunnan is generally around 3000 m, and the discrepancies between these sites and others can reach 2000-3000 m. In the ATML calculation, the atmospheric pressure data usually do not reflect the real atmospheric pressure in the Sichuan-Yunnan region. The ATML is therefore often inaccurately estimated in areas with complex terrain changes such as the Sichuan-Yunnan region [24].
We computed the averaged NTOL and ATML displacements and their RMS values at each site, as shown in Figure 12. The effect of the NTOL is found to be limited in the Sichuan-Yunnan region, the amplitude of which is only approximately 1 mm. The period from April to August is found to have large ATML variations, which is consistent with the time period of the TWS anomalies that occurred in the northwestern Sichuan-Yunnan region. If the ATML during this period in the northwestern Sichuan-Yunnan region is overestimated, the elastic responses observed by GPS would also be overestimated, which would be thought to be hydrological deformations. The inverted TWS will therefore be overestimated, which reasonably explains why there are anomalies in the northwestern Sichuan-Yunnan region from April to August. : annual amplitude for each series in the entire region; : annual amplitude for each series in northwestern Sichuan-Yunnan

Discussion
The anomalies of the GPS-derived TWS variations in northwestern Sichuan-Yunnan indicate that the coordinate time series used in the inversion contains deformation information or errors other than hydrological deformation. Considering the fact that northwestern Sichuan-Yunnan is near the Qinghai-Tibet Plateau and the terrain is extremely complex (Figure 1), the inaccurate ATML correction in this area may be the cause of anomalies in the GPS inversion result. The height of the GPS sites in northwestern Sichuan-Yunnan is generally around 3000 m, and the discrepancies between these sites and others can reach 2000-3000 m. In the ATML calculation, the atmospheric pressure data usually do not reflect the real atmospheric pressure in the Sichuan-Yunnan region. The ATML is therefore often inaccurately estimated in areas with complex terrain changes such as the Sichuan-Yunnan region [24].
We computed the averaged NTOL and ATML displacements and their RMS values at each site, as shown in Figure 12. The effect of the NTOL is found to be limited in the Sichuan-Yunnan region, the amplitude of which is only approximately 1 mm. The period from April to August is found to have large ATML variations, which is consistent with the time period of the TWS anomalies that occurred in the northwestern Sichuan-Yunnan region. If the ATML during this period in the northwestern Sichuan-Yunnan region is overestimated, the elastic responses observed by GPS would also be overestimated, which would be thought to be hydrological deformations. The inverted TWS will therefore be overestimated, which reasonably explains why there are anomalies in the northwestern Sichuan-Yunnan region from April to August. Another possible reason for the anomaly in northwestern Sichuan-Yunnan may be the effect of high topographic variability that can cause more residual tropospheric errors Another possible reason for the anomaly in northwestern Sichuan-Yunnan may be the effect of high topographic variability that can cause more residual tropospheric errors in processing GPS observations.

Conclusions
GPS can be used to measure elastic deformation caused by mass loading variations on the Earth's surface. Accurate hydrological deformation information from the GPS coordinate time series is the most critical step in TWS inversion. We propose a TWS inversion method using GPS observations based on independent component analysis (ICA). The ICA method is first used to separate the hydrological deformations from GPS coordinate time series. The TWS variations are then quantitatively estimated using the separated hydrological deformation displacements.
We estimate the TWS variations in the Sichuan-Yunnan region in China using different approaches, including GPS, GPS(ICA), GRACE, and GLDAS methods. The overall results of the four methods show similar change patterns in the spatial and temporal domain. Some differences are found in the results due to the completely different observation methods and implied sampling and process-representations. The TWS variations inverted by GPS and GPS(ICA) show more notable variations compared with the GRACE and GLDAS results. We find that anomalies occur in the northwestern Sichuan-Yunnan region for the GPS inversion results, where GPS sites are located within extreme and complex terrain changes.
We quantitatively evaluate the differences between the amplitude of the four different method-derived TWS variations in northwestern Sichuan-Yunnan and the entire region and find that the annual amplitude of the GPS-inverted TWS in the northwestern region is nearly twice that of the entire region. The anomalies in the GPS inversion results are further explored, and we suggest that the extreme and complex terrain changes may be the main cause of the GPS inverted TWS anomalies in the northwestern Sichuan-Yunnan region. The proposed GPS(ICA) method attenuates this problem and more accurately estimates the TWS variation results in the Sichuan-Yunnan region, which provides a new method for estimating regional TWS variations using GPS observations.