Inversion of GNSS Vertical Displacements for Terrestrial Water Storage Changes Using Slepian Basis Functions

The surface displacements measured by the Global Navigation Satellite System (GNSS) provide a unique insight for studying terrestrial water storage (TWS) changes. In this study, we recovered the TWS changes from GNSS vertical displacements in Southwest China (SWC) using Slepian basis function (SBF) from January 2011 to December 2020. The performance of the TWS changes estimated by SBF was validated against the Gravity Recovery and Climate Experiment (GRACE)/GRACE Follow‐On (GFO) and the GNSS‐inverted TWS changes estimated by Green's function method. The results showed that the spatial patterns, seasonal, and linear trends of the TWS changes derived from GNSS using SBF agreed with the GRACE/GFO estimates. However, there are still evident differences in the local scope, and the GNSS‐derived TWS changes presented stronger amplitudes and more details in the spatio‐temporal domains than the GRACE/GFO estimates in SWC. The unconstrained GNSS inversion results using SBF also presented stronger signal amplitudes than those estimates with Green's function, and the TWS changes estimated by SBF were more reliable than Green's function method for regions with sparsely distributed GNSS stations in SWC. Additionally, the average distance between the GNSS stations can be considered as a reasonable filtering radius of SBF, and the SBF‐estimated TWS changes with different Gaussian filtering radii had comparable signal amplitudes and spatial patterns with the estimates of Green's function and GRACE/GFO.

For the SBF method, Han and Razeghi (2017) first introduced the SBF into the inversion of surface mass changes using GNSS vertical displacements in the Australian continent, and the results demonstrated that the GNSS-inverted surface mass changes were in good agreement with the GRACE and hydrologic model estimates. Based on the research of Han and Razeghi (2017), Jiang, Hsu, Yuan, Cheng, et al. (2021), and Jiang et al. (2022) inverted the daily GNSS vertical displacements to estimate TWS changes using the SBF in mainland China and Brazil, and the GNSS-inferred TWS changes could clearly reveal the hydrological drought events in these two regions. Cheng et al. (2021) investigated the TWS changes in Sichuan, Yunnan, and Chongqing in China derived from GNSS observations using the SBF, and the GNSS-derived TWS changes presented better consistency with precipitation than the estimated results of GRACE and hydrologic model. Previous studies have focused on validating the performance and application of SBF for recovering TWS changes though external remote sensing and hydrometeorological data (e.g., GRACE/GFO, hydrologic model, and precipitation). However, other features and advantages of SBF for recovering TWS changes from GNSS observations should be further investigated. For instance, how should we determine the maximum degree and filtering radius of SBF corresponding to the specific study region and GNSS stations? What are the differences in the TWS changes from GNSS observations between the SBF and Green's function methods?
In this study, we investigated the performance of SBF for recovering TWS changes in Southwest China (SWC). The vertical displacement time series of 83 GNSS stations from the Crustal Movement Observation Network of China (CMONOC) were used to infer the TWS changes for the period of January 2011 to December 2020. First, we conducted a closed-loop simulation to determine the maximum degree of the SBF and the optimal concentration factor of the study region, and investigate the impacts of different station distributions on GNSS inversion results. Second, the SBF-derived TWS changes from GNSS vertical displacements were validated against the GRACE/GFO mascon solutions in the spatio-temporal domains. Third, we determined the optimal filtering radius of the SBF and analyzed the differences in the TWS changes derived from GNSS vertical displacements between the SBF and Green's function methods in the spatio-temporal domains.
The rest of this study is organized as follows: Section 2 provides an overview of the study region, the data set resources and the corresponding post-processing. Section 3 presents the theory of SBF and introduces an inversion model for inferring TWS changes from GNSS vertical displacements using SBF. Section 4 presents the closed-loop simulation and validates the performance of the SBF method with GRACE/GFO estimates in the spatio-temporal domains. Section 5 compares the TWS changes derived from GNSS using SBF with different Gaussian filtering radii and the Green's function method with a priori constraints. Additionally, the differences between the GNSS-inverted TWS changes based on the SBF and GRACE/GFO estimates are further discussed. Finally, some conclusions and suggestions are given in Section 6. 10.1029/2022EA002608 3 of 18

Study Region
SWC (21°08'-34°21'N and 97°21'-110°11'E) is located in the upper reaches of the Yangtze River in the southeastern part of the Qinghai-Tibet Plateau, which includes Yunnan Province (YNP), Sichuan Province (SCP), Guizhou Province (GZP), and Chongqing Municipality (CQM) (see Figure 1) . Figure 1 shows the geographic overview of SWC, and the locations of the major rivers and GNSS stations. The topography of SWC is high in the northwest and low in the southeast. The study region is dominated by plateaus and mountains, and the karst, river valley and basin landforms are widely distributed over this region. SWC has a subtropical monsoon climate, with an average annual temperature of ∼14°C-24°C and precipitation of ∼600-2,300 mm. Additionally, SWC possesses abundant natural resources, for example, water, biological, mineral, and land resources. Therefore, accurately investigating the spatio-temporal characteristic of TWS changes in SWC is of great significance for dealing with water shortages, climate changes, and regional socioeconomic developments.

GNSS Data and Post-Processing
In this study, we inverted the vertical displacement time series of 83 GNSS stations from the Crustal Movement Observation Network of China (CMONOC) to obtain the TWS changes in SWC from January 2011 to December 2020. The geographical locations of these 83 GNSS stations are presented in Figure 1, and the average station distance ( √

GNSS
, where S and n GNSS are the area (unit: km 2 ) and the number of GNSS stations within the study region) between GNSS stations within SWC is ∼150 km. These daily GNSS vertical displacement time series were obtained from the National Earthquake Data Center (NEDC) and processed at the First Monitoring and Application Center, China Earthquake Administration. The GNSS vertical displacement time series published by the NEDC were processed by GAMIT/GLOBK software (version 10.4) (Herring et al., 2010) and the Quasi-observation Combination Analysis (QOCA) software (Dong et al., 1998) based on the double-differenced model. First, the daily free network solutions of the CMONOC and the International GNSS Service core stations were processed using GAMIT software, and all appropriate corrections were applied in this procedure: the GNSS receiver and satellite antennas were calibrated by the absolute antenna phase center model (Schmid et al., 2016); the solid Earth tides and pole tides were corrected according to the International Earth Rotation Service conventions 2003 and the ocean tidal loading was corrected using the Finite Element Solutions 2004 model (Lyard et al., 2006); the hydrostatic and wet zenith delays were estimated by the gridded Vienna Mapping Function and the tropospheric delays were modeled using Global Mapping Function (Böhm et al., 2006). Second, the QOCA software was used for joint adjustment to obtain the final GNSS displacement time series relative to the International Terrestrial Reference Frame 2008 (Altamimi et al., 2011).
The downloaded GNSS time series include linear signals (e.g., Glacial Isostatic Adjustment [GIA], tectonic motion and the porous response of solid Earth to groundwater filling and emptying the aquifers), seasonal signals (e.g., TWS loading, non-tidal atmospheric loading [NTAL], non-tidal oceanic loading [NTOL], draconitic term and thermal expansion) and other components (e.g., outliers, offsets, and postseismic deformation). In order to process the data jumps, outliers and discontinuities of the GNSS data due to the influences of human activities, equipment replacement, and natural environmental changes, the data post-processing was performed on these GNSS observations to obtain accurate hydrological loading displacement time series. First, the least squares fitting method was used to remove the data jumps. Second, the quartile gross error detection method proposed by Nikolaidis (2002) was used to remove the gross errors. Third, the Kriged Kalman filter interpolation software proposed by N. Liu et al. (2018) was used to interpolate the missing data of the GNSS observations. Figure 2 shows the vertical displacement time series of four representative GNSS stations (i.e., GZGY, SCDF, YNSM, and YNYL) before and after post-processing; the geographical locations of these four GNSS stations are marked in Figure 1. As shown in Figure 2, the post-processing could well eliminate data jumps and gross errors and fill in the missing data.
The surface mass changes are dominated by seasonal signals in the study region (see Figure 2), and the long-term trend mainly includes the potential contributions of tectonic motion and GIA effects (Fu et al., 2015). Therefore, to better understand the seasonal characteristics of the surface mass changes in SWC, we detrended the GNSS observations to remove potential tectonic motion or GIA effects. Moreover, to obtain the final TWS changes from the GNSS observations, the NTAL and NTOL effects were corrected by the geophysical fluid loading products published by the Earth System Modeling group at Deutsches GeoForschungsZentrum (ESMGFZ) (Dill & Dobslaw, 2013).

GRACE/GFO Data
To assess the performance of the SBF method for recovering the TWS changes from GNSS vertical displacements, we used the GRACE/GFO RL06 monthly mascon solutions (Save, 2019) to infer the TWS changes in SWC over the period from January 2011 to June 2017 and from June 2018 to December 2020 as a comparison. Given the large uncertainties in the low-degree terms of the GRACE/GFO observations, various corrections have been made to the released GRACE/GFO mascon products. For instance, the terms C 20 and C 30 were replaced by the Satellite Laser Ranging estimates (Loomis et al., 2020), and the degree-1 coefficients were replaced by the estimates from the GRACE project as supplementary datasets (GRACE Technical Note 13) (Landerer, 2019). Moreover, the GIA effects were corrected by the ICE-6G_D model (Peltier et al., 2018).

Global Land Data Assimilation System Data
The Global Land Data Assimilation System (GLDAS) hydrologic model is jointly established by the National Aeronautics and Space Administration's Goddard Space Flight Center and National Oceanic and Atmospheric Administration's National Centers of Environmental Prediction (Rodell et al., 2004). The GLDAS derives four models (i.e., NOAH, VIC, CLM, and MOSAIC), which integrates a variety of land surface flux information in the form of grids (soil moisture, soil temperature, evaporation, rainfall, runoff and snow, etc.). In order to investigate the performance of SBF, the TWS changes (the sum of canopy water, soil moisture (0-2 m) and snow water equivalent) from GLDAS_NOAH_V10 products with the spatial resolution of 1° × 1° from January 2011 to December 2020 were used to generate the synthetic TWS changes for closedloop simulation.

Theory of Slepian Basis Function
Geophysical fields (e.g., gravitational, magnetic, and electric fields) on or outside the Earth's sphere can be represented by a series of surface spherical harmonic coefficients (SHCs), and determination of the surface SHCs depends on observational data with global coverage. Nevertheless, for a space-limited region, the surface spherical harmonic functions (SHFs) do not satisfy orthogonality due to the limited observation data. Therefore, other appropriate and reliable methods should be introduced to represent these geophysical signals. Slepian (1983) first proposed the Slepian problem, mainly to solve the energy concentration problem in the time and frequency domain in one-dimensional continuity space. Subsequently, the Slepian problem was introduced into three-dimensional space, and had been used to solve related scientific problems in the field of geophysics and cosmology (Dahlen & Simons, 2008;. The SBF ( ( , ) ) on a sphere can be represented by the surface SHFs  as follows: where and are the colatitude and longitude of the grid cells within the study region, and are the degree and order of the SHCs, L is the maximum expansion degree, , represents the SHCs of the expansion of ( , ) , and ( , ) represents the surface SHFs. To ensure that the energy of ( , ) is concentrated in the local region R (i.e., the spectral energy of the ( , ) is maximized within the region R), ( , ) should satisfy ): where < < represents the concentration factor, which indicates the concentration of ( , ) within the R, such that, the closer is to 1, the more concentrated ( , ) is within R; represent the entire sphere. The maximization problem of in Equation 2 is called the Slepian problem. The solutions of the Slepian problem can be expressed as a system of eigenvalue equations in the frequency domain as follows: The concentration factor and the Slepian harmonic expansion coefficient , are the solutions of the eigenvalues and eigenvectors of Equation 3, respectively.
According to the size of the study region and the distributions of the GNSS stations, the maximum expansion degree of the SBFs was expanded to degree-65 in SWC (see Section 4.1). Therefore, there are (65 + 1) 2 = 4,356 eigenvalues (i.e., the concentration factor ) and SBFs for our study region. Figure 3 shows the concentration factor corresponding to different values of the index . There are 30 concentration factors greater than 0.1, and the optimal index for the study region is 15 (see Section 4.1) and the corresponding is 0.7214. Additionally, the total contribution of the first 15 concentration factors used in this study is 70.10%. Furthermore, Figure 4 shows the spatial patterns of examples for the first 15 SBFs concentrated in SWC. We can see that the SBFs signals gradually attenuate outside the extended region with an increase in the index , and the first 15 SBFs can

Inversion Model for Estimating TWS Changes Using GNSS
According to the characteristics of the local concentration of SBF, the land motions observed by GNSS in a specific study region (e.g., SWC) can be represented by a series of SBFs as follows (Han & Razeghi, 2017): where are the GNSS vertical displacement observations, a (6,378.137 km) is the average radius of the Earth and is the Slepian coefficient corresponding to the α-th SBF ( , ) . The linear model of Equation 4 can be expressed as follows: where is the GNSS vertical displacement observation vector, is the coefficient matrix corresponding to the SBF ( , ) determined by the study region and geographical locations of the GNSS stations, is the Slepian coefficient vector (i.e., ) to be estimated, and is the GNSS observation error vector. Hence, the estimates of can be expressed through least squares adjustment as follows: According to the elastic load theory (Farrell, 1972), the Slepian coefficient corresponding to the vertical displacements can be converted into the Slepian coefficient corresponding to the mass loads by load Love numbers as follows (Han & Razeghi, 2017): where (5,517 km/m 3 ) is the average density of the Earth, (1,000 km/m 3 ) is the density of fresh water and ′ represent the load Love numbers considered from Wang et al. (2012). The eigenvalue-weighted regional water mass loads, in terms of equivalent water height, can be expressed as follows (Han & Razeghi, 2017): where J is the optimal concentration factor.

Green's Function Method
Considering that the inversion algorithm of the Green's function method is relatively mature enough, we employed the estimated results of Green's function method to assess the performance of SBF-estimated TWS changes in this study. Unlike the SBF method, the Green's function method represents the surface deformation by the convolution of the mass loads and Green's function in the space domain (Farrell, 1972) as follows: where ( ) is the surface vertical displacement, is the angular distance; ∆ is the mass load; (5.965 × 10 24 kg) is the mass of the Earth; and are the Legendre polynomials.
By establishing the relationship between the mass loads and vertical displacements through Green's function, the surface mass changes (e.g., TWS changes) can be recovered using the least squares adjustment. It should be noted that inversion of TWS changes from GNSS observations using Green's function method is a serious ill-posed problem, and the regularization method like Tikhonov and Truncated Singular Value Decomposition methods (Argus et al., 2014;Lai et al., 2020) should be introduced to stabilize the inversion results. The objective function of Green's function method based on the Tikhonov regularization is as follows (Argus et al., 2014;Fu et al., 2015): where is the design matrix associated with Green's function; is the unknown TWS change vector that needs to be estimated; is the GNSS observation vector, is the noise standard deviation (STD) of the GNSS observations; is the Laplacian matrix; and is the regularization parameter. Then the unknown TWS change vector can be determined as follows:

Closed-Loop Simulation
Determining the optimal maximum degree of SBF and the corresponding concentration factor of the study region helps us to obtain more reliable TWS changes derived from GNSS observations based on SBF. For this purpose, we conducted a closed-loop simulation study. We fitted and obtained the annual amplitudes of TWS change time series derived from GLDAS hydrologic model from January 2011 to December 2020 at each grid cell in SWC, and used the GLDAS-based annual amplitudes to generate the synthetic TWS changes (synthetic signals). To reduce the edge effects and make use of more GNSS station data as much as possible, we extended the boundaries of the study region as 2° (see the red curve in Figure 1). Then we used the synthetic signals to simulate the vertical displacements at 83 GNSS stations by Equation 9. We determined the maximum degree and optimal concentration factor by comparing the correlation coefficient (CC) and STD between the synthetic signals and inversion results from the simulated 83 GNSS observations, the maximum degree and concentration factor corresponding to the minimum STD and highest CC were considered as the optimal values. Figure 5 shows the CCs and STDs and the optimal concentration factors based on different maximum degrees (i.e., from degree-40 to degree-90 with 5° interval). As shown in Figure 5, the CCs (black curve) between the synthetic signals and inversion results corresponding to these different maximum degrees generally keep a high level (i.e., between 0.74 and 0.85), and the CCs of degree-60, -65, and -70 are the largest (i.e., 0.85). The STDs (red curve) corresponding from degree-40 to degree-75 are 15-20 mm, but with the maximum degrees greater than degree-75, the STDs increase and the corresponding CCs decrease rapidly. This is because the using of higher maximum degree makes more parameters to be estimated in Equation 5, which leads to the instability of estimated Slepian coefficients. For the optimal concentration factors (blue curve), they vary between 14 and 15. In general, the using of maximum degrees between degree-45 and degree-70 can well recover the synthetic signals derived from the simulated GNSS observations based on SBF, the corresponding differences of STDs and CCs between them are less than 0.96 and 0.02 mm, respectively. We chose the degree-65 and 15, whose CC and STD (i.e., 0.85 and 16.61 mm) are better than other degrees, as the maximum degree and optimal concentration factor of our measured GNSS inversion below. It is noteworthy that the minimum and maximum distances between adjacent stations within the study region are ∼30 and ∼260 km, respectively, and the average station distance is ∼150 km. The degree 65 is generally close to the ratio of 20,000/260 km (20,000 is a half of the circumference of Earth). Additionally, the first 15 concentration factors can well express the TWS signal in SWC (see Figure 4).
Additionally, we further investigated the impacts of GNSS station distributions (i.e., the evenly/unevenly and sparsely/densely distributed stations) on the inversion results based on SBF in SWC. For the evenly and unevenly distributed stations test, given that the GNSS stations within SWC (i.e., 57 stations in this study) are the main contributors to the final inversion results, we only set these 57 stations as evenly distributed stations (e.g., corresponding 1.5° interval) and compared them with the unevenly distributed stations in the actual situation. For the sparsely and densely distributed stations test, we designed the GNSS stations within SWC with the intervals of 2°, 1°, and 0.5° and the corresponding numbers of GNSS stations within SWC are 26, 105, and 421, respectively. Note that the number of GNSS stations between the SWC and the extended boundary is 26. In addition, to be consistent with the real GNSS vertical displacements, which are contaminated by observation noises, and considering that the surface mass loads and thermal expansion effects account for ∼60% of the observed amplitude of GNSS vertical displacements (Yan et al., 2009), Gaussian white noise with a STD of 0.46 mm (i.e., occupying 40% of the simulated GNSS observations) was added to the simulated GNSS observations. Figure 6 shows the patterns of the synthetic signals and the inversion results from simulated GNSS observations using SBF with different station distributions, and the corresponding STDs and CCs between the synthetic signals and different GNSS inversion results are also presented. As shown in Figures 6b and 6c, for the comparison of evenly and unevenly distributed stations, the spatial patterns of TWS changes derived from the actual station distributions and even station distributions within SWC present good consistency with the synthetic signals. However, the estimates of even station distributions are closer to the synthetic signals than those of actual station In general, the 83 GNSS stations can be well used to recover the TWS changes in SWC.

Spatial Patterns of TWS Changes
To further assess the performance of SBF for recovering TWS changes in SWC, we inverted the measured GNSS vertical displacement time series of 83 GNSS stations for TWS changes from January 2011 to December 2020. In addition, the GRACE/GFO-derived TWS changes are presented for comparison. It should be noted that the TWS changes in SWC derived from GNSS observations using SBF do not add any external information as constraints.
The SBF-estimated TWS changes include the comprehensive contributions of water loads and other noises. Therefore, to suppress the possible errors in the GNSS inversion results, we applied a 150-km (corresponding to the average distance between GNSS stations within SWC) Gaussian filtering  to the GNSS inversion results. Figure 7 illustrates the spatial patterns of the TWS changes derived from GNSS and GRACE/ GFO in August 2011, April 2012, September 2019, and March 2020. As illustrated in Figure 7, in August 2011, the GNSS-and GRACE/GFO-derived TWS changes mainly presented negative and positive anomalies in the eastern and western parts of SWC, respectively. However, the GNSS-derived TWS changes were more significant and stronger than those of GRACE estimates. In April 2012, both GNSS and GRACE well detected the negative anomalies in the southern part of SWC, but the positive anomalies in the northwest and northeast parts of SWC were not well detected by GRACE. In September 2019, the TWS changes in all of SWC generally presented positive anomalies, but the amplitudes of the GNSS estimates were stronger than those of the GFO estimates. In March 2020, the spatial patterns of TWS changes from GNSS and GFO were similar to those in April 2012, and the GNSS-derived TWS changes had stronger amplitudes in the southwest part of SWC than those of GFO, while the GFO-derived TWS changes had stronger amplitudes in the northeast part of SWC than those of GNSS. In general, the spatial patterns of TWS changes derived from GNSS and GRACE/GFO are generally consistent with each other, but there are still evident differences in their local scope; in particular, the amplitudes of GNSS estimates are stronger than those of GRACE/GFO estimates. Figure 8 further presents the spatial patterns of annual amplitude for the TWS changes derived from GNSS based on SBF and GRACE/GFO from January 2011 to December 2020. In particular, the annual amplitudes of the TWS changes derived from GRACE/GFO were simultaneously fitted using GRACE and GFO estimates. As presented in Figure 8, the annual amplitudes of the TWS changes from GNSS and GRACE/GFO generally keep good consistency with each other, and the TWS changes gradually increase from northeast to southwest of SWC; the corresponding average annual amplitudes are 146.11 ± 3.86 and 72.61 ± 12.56 mm, respectively. This demonstrates that the annual amplitudes of GNSS estimates are significantly stronger than those of GRACE/ GFO estimates, although the GNSS inversion results have been smoothed by a 150-km Gaussian filtering. This is because the ground-based GNSS observations are more sensitive to changes in regional mass loads within 100-km around each GNSS station (Bevis et al., 2005), while the satellite-based GRACE/GFO observations are more sensitive to large-scale mass loads.

Time-Variant TWS Changes
Investigating the time-variant TWS changes helps us better understand the regional water cycle of SWC. We calculated the region-averaged TWS change time series and fitted them with a position, a linear trend and a seasonal term (i.e., = + + ( π + ) , where and t are the time series and time, a and b are the position and linear trend, and and are annual amplitude and phase, respectively.). Figure 9 presents the time series of raw TWS changes, and the seasonal and linear trends derived from daily GNSS and monthly GRACE/GFO observations over the period of January 2011-December 2020. Table 1 lists the corresponding annual amplitudes, annual phases and linear trends of the TWS change time series derived from GNSS and  GRACE/GFO during the study period. As presented in Figure 9a, both the GNSS-and GRACE/GFO-derived TWS changes presented obvious seasonal characteristics over the study region, and the CC and STD between the GNSS and GRACE/GFO estimates were 0.68 and 87.96 mm, respectively. It should be noted that the daily GNSS-derived TWS changes presented a stronger signal amplitude and more details than those of GRACE/ GFO estimates because GNSS is more sensitive to regional TWS changes than GRACE/GFO, and the amplitude of GNSS is almost twice that of GRACE/GFO (see Table 1). For the annual phases, the GRACE/GFO-derived TWS change time series lagged behind the GNSS estimates by more than 1 month, which demonstrates that GNSS is also more responsive to TWS changes than GRACE/GFO in the time domain (Jiang, Hsu, Yuan, Cheng, et al., 2021). For the linear trends, the TWS changes mainly showed an increasing trend from January 2011 to June 2017 (corresponding to the GRACE operation period) and a decreasing trend from July 2017 (June 2018 for GFO) to December 2020 (corresponding to the GFO operation period). For the GRACE operation period, the increasing rate of TWS changes derived from GRACE solutions (i.e., 7.62 ± 2.73 mm/yr) is greater than that of GNSS estimates (i.e., 3.69 ± 1.21 mm/yr); while for the GFO operation period, the decreasing rate of TWS changes derived from GNSS solutions (i.e., −18.78 ± 2.81 mm/yr) is greater than that of the GFO estimates (i.e., −14.10 ± 14.16 mm/yr). In general, although there are still some differences in the annual amplitudes, annual phases and linear trend rates of the TWS changes derived from GNSS and GRACE/GFO solutions, all of them can clearly reveal the seasonal and long-term variation characteristics in SWC during the study period.
As mentioned above, GNSS observations are more sensitive to regional TWS changes than GRACE/GFO observations in the spatial domain. Another advantage of GNSS over GRACE/GFO observations is that GNSS has a daily or even sub-daily time resolution. To better compare the differences in the TWS change time series from GNSS and GRACE observations (considering the large gap between the GRACE and GFO, we only used GRACE data here), we performed a power spectrum analysis on the TWS change time series derived from these two geodetic technologies from January 2011 to June 2017 using Fast Fourier Transformation (FFT). The missing months' observations during the operation period of the GRACE mission were interpolated using cubic splines before FFT. Figure 10 illustrates the power spectrum of the TWS change time series derived from GNSS and GRACE. We can see that both the GNSS and GRACE well reveal the significant seasonal changes (mainly corresponding to the annual and semi-annual changes), but the daily GNSS observations reveal more high-frequency signals than those of GRACE observations. The power of the GNSS estimates was significantly larger than that of the GRACE estimates, this is because the ground-based GNSS in more sensitive to the TWS changes than satellite-based GRACE and the corresponding signal amplitude of TWS change derived from GNSS is stronger (see Figure 9 and Table 1). Additionally, the GNSS is more sensitive to the semi-annual variations of TWS than GRACE due to the higher spatio-temporal resolution.

Discussion
To further investigate the performance of SBF for recovering TWS changes in SWC, we compared the SBF-estimated TWS changes with different Gaussian filtering radii and the Green's function-derived TWS changes with a priori constraints (i.e., the Laplacian matrix constraint) (Argus et al., 2014). First, based on the closed-loop simulation results of Section 4.1, we further used the GLDAS hydrologic model from January 2011 to December 2020 to generate the synthetic signals and simulate the vertical displacement time  , and then the TWS changes derived from the simulated GNSS observations with errors using the Green's function and SBF method with different filtering radii were used for comparison. Figure 11 displays the comparison of TWS changes from the synthetic signals and the GNSS-inverted TWS changes using Green's function and SBF method with no, 150 and 220-km Gaussian filtering. As shown in Figure 11, the unconstrained SBF-estimated TWS changes derived from the simulated GNSS observations with no filtering contain the comprehensive contributions of input synthetic signals and errors, which presents a larger amplitude than that of synthetic signals. However, the TWS changes derived from the simulated GNSS observations estimated by the Green's function method presents a smaller amplitude than that of the synthetic signals owing to the regularization constraint. It should be noted that the SBF-derived TWS changes with 150 and 220-km Gaussian filtering were close to the synthetic signals and estimated results from Green's function method in SWC, respectively. It demonstrated that using the average distance between the GNSS stations within the study region as the Gaussian filtering radius is reasonable. In addition, the TWS changes estimated by Green's function method with Laplacian matrix constraint is comparable to the SBF estimates with 220-km Gaussian filtering in SWC.
Additionally, we further explored the performance of different Gaussian filtering radii for recovering TWS changes from measured GNSS observations using SBF method. Figure 12 shows the spatial patterns of the annual amplitudes for TWS changes derived from measured GNSS observations using SBF with different Gaussian filtering radii and Green's function with a Laplacian smoothing constraint; the GRACE/GFO-derived annual amplitudes are also presented for comparison. As shown in Figure 12, for the SBF-estimated TWS changes with no filtering (see Figure 12a), the annual amplitudes presented in the western and southern parts of SWC were very significant, and the corresponding amplitudes are larger than 400 mm. This is because the TWS changes estimated from GNSS vertical displacements using SBF without any constraints contain the comprehensive contributions of actual TWS changes and other errors. Furthermore, Gaussian filtering with different filtering radii was used to reduce the high-frequency noise. We can see that the SBF-estimated annual amplitudes with different filtering radii (see Figures 12b, 12c, 12e, and 12f) have analogous spatial patterns to the results of no filtering (see Figure 12a). The Gaussian filtering suppresses the high-frequency noise and simultaneously weakens the signal amplitudes of the actual TWS changes. The annual amplitudes of SBF-derived TWS changes with a 220-km Gaussian filtering (see Figure 12c) are comparable in terms of the magnitude and spatial patterns of the Green's function-derived TWS changes with a priori constraints (i.e., Figure 12c vs. Figure 12d) and the corresponding average annual amplitudes are 104.02 ± 2.60 and 101.40 ± 3.90 mm, respectively, but the SBF estimates are smoother than the Green's function estimates. This is consistent with the simulation results presented in Figure 11. It should be noted that the inversion of TWS changes from the GNSS vertical displacements using the Green's function method is a typical discrete ill-posed problem, so an external priori constraint should be introduced to stabilize the inversion results. The second-order difference Laplacian operator (i.e., the Laplacian matrix) is commonly used as the a priori constraints, which can also attenuate the signal amplitude along with the noise, similar to Gaussian filtering. The a priori constraints in the regularization process is equivalent to the post-processing filter (e.g., Gaussian filtering) (Swenson & Wahr, 2011). However, the inversion algorithm for inverting the GNSS vertical displacements for TWS changes using SBF is relatively stable and the external constraint is not necessary. In general, compared with the Green's function method, the unconstrained SBF-estimated GNSS inversion results without filtering present stronger signal amplitudes for the TWS changes in SWC, but the contributions of potential errors should be considered. Additionally, the SBF-estimated GNSS inversion results with a 150-km Gaussian filtering are regarded as reasonable in SWC, and the SBF-derived TWS changes with a 220-km Gaussian filtering present comparable signal amplitudes to the estimates of Green's function in SWC. It should be noted that the proper filtering radius depends on the station coverage and distribution in different study regions.
As mentioned above, SBF-estimated TWS changes with a 220-km Gaussian filtering and Green's function-derived TWS changes with a priori constraints have comparable annual amplitudes. To further investigate the differences in the TWS changes derived from GNSS observations using SBF with a 220-km Gaussian filtering and those derived by Green's function with a priori constraints in the time domain, Figure 13 depicts the TWS change time series in YNP, SCP, GZP, and CQM derived by these two inversion methods. As presented in Figure 13, the TWS changes estimated by SBF with a 220-km Gaussian filtering and those of Green's function with a priori constraints are in good agreement in these four provinces (especially for the YNP and SCP), the corresponding CCs are larger than 0.85 and the STDs are between 22.90 and 36.56 mm (see Table 2). Nevertheless, there are still differences in the TWS change time series of these two inversion methods: the inversion results based on the Green's function method for GZP and CQM have some outliers from 2014 to 2016, and the corresponding CCs and STDs of the TWS changes in these two provinces from SBF and Green's function are 0.89 and 36.56 mm, 0.85 and 32.40 mm, respectively. This is because the GNSS stations of GZP and CQM are sparser than those in the other two provinces. It should be pointed out that the Green's function method is more applicable to the regions where the GNSS stations are densely distributed (Han & Razeghi, 2017), and the sparse GNSS stations can introduce great uncertainty to Green's function estimates. Additionally, accurately estimating the optimal regularization parameter of the Green's function inversion method is still a challenging work. However, the unconstrained SBF-estimated TWS with 150-km Gaussian filtering can well weaken these potential effects. It can be concluded that the SBF method is a more suitable method for recovering the TWS changes from GNSS observations than the Green's function method in the regions with sparsely distributed GNSS stations in SWC.
As mentioned in Sections 4.2 and 4.3, although the SBF-derived TWS changes were smoothed by a 150-km Gaussian filtering, the annual amplitude of SBF-estimated TWS changes is still approximately twice that of GRACE/ GFO estimates. To obtain the amplitude of SBF-estimated TWS changes close to that of the GRACE/ GFO estimates, we further smoothed the SBF-estimated TWS changes by a Gaussian filtering with a larger radius (see Figures 12e and 12f). As shown in Figures 12e and 12f, the signal amplitude of the TWS changes was severely suppressed with an increase in the filtering radius. Among them, the SBF-estimated TWS changes with a 300-km filtering radius had a comparable amplitude to the GRACE/GFO estimates. Figure 14 shows the TWS change time series estimated by SBF with a 300-km Gaussian filtering and GRACE/GFO estimates for SWC. We can see that the amplitude of the SBF-estimated TWS changes with a 300-km Gaussian filtering is relatively close to that of the GRACE/GFO estimates, and the corresponding annual amplitudes are 68.11 ± 1.30 and 69.54 ± 7.29 mm, respectively. This demonstrates that although the mascon products of GRACE/GFO have greater spatial resolution for recovering TWS changes in small-scale regions than the spherical harmonic solutions (Watkins et al., 2015), the signal amplitude and spatio-temporal resolution of the TWS changes from GRACE/GFO observations are still limited by the native accuracy of GRACE/GFO. This further verifies the advantages of GNSS observation technology for recovering TWS changes in small-scale regions. In general, the unconstrained TWS changes from the GNSS vertical displacements using the SBF method present stronger signal amplitudes than those of GNSS-inverted TWS changes  using the Green's function method and GRACE/GFO estimates, and the SBF-estimated TWS changes smoothed by different Gaussian filtering radii can present comparable spatial patterns and signal amplitudes to those estimated by Green's function and GRACE/GFO. However, the best way to objectively validate the performance of SBF and extract the actual water mass loads of SBF-derived TWS changes still needs to be further studied.

Conclusions
In this study, we investigated the performance of SBF for inferring TWS changes from GNSS vertical displacements in SWC. We used the vertical displacement time series of 83 GNSS stations from CMONOC to infer the TWS changes in SWC, for the period of January 2011-December 2020. The performance of TWS changes derived from GNSS using SBF was validated against the GRACE/GFO mascon products and GNSS-inverted TWS changes using the Green's function method in the spatio-temporal domains.
The closed-loop simulation study showed that the maximum degree and optimal concentration factor of SBF was degree-65 and 15 in SWC, respectively. Additionally, the evenly distributed stations within SWC can recover more reliable TWS changes than those of the actual distribution stations, and dense GNSS stations can obtain more stable inversion results to a certain extent. Validations with the GRACE/GFO estimates indicated that the spatial patterns and the seasonal and linear trends of TWS changes derived from GNSS observations using SBF and GRACE/GFO were in good agreement. However, SBF-estimated TWS changes from GNSS presented stronger amplitudes and higher time resolution (i.e., daily) in the spatio-temporal domain than those of GRACE/ GFO estimates. Additionally, the TWS changes derived from GNSS using SBF smoothed by a 300-km Gaussian filtering had comparable amplitudes with the GRACE/GFO estimates, which further indicates the advantages of ground-based GNSS observations for recovering TWS changes over satellite-based GRACE/GFO observations in the small-scale regions. Furthermore, it demonstrated that the SBF-estimated TWS changes without any priori constraints also presented stronger amplitudes of the TWS changes than the results estimated by the Green's function method, and the TWS changes derived from GNSS using SBF smoothed by a 220-km Gaussian filtering had comparable amplitudes with those of the Green's function method with a priori constraints (i.e., the Laplacian matrix) in this study. In addition, the average distance between the GNSS stations (i.e., 150-km in this study) can be considered as a reasonable filtering radius in SWC, and the TWS changes derived from GNSS observations using SBF presented better stability than those of the Green's function method for the regions with sparsely distributed GNSS stations.
The unconstrained TWS changes derived from GNSS using SBF contain the comprehensive contributions of actual water mass loads and errors. We employed the averaged distance between the GNSS stations as the Gaussian filtering radius in this work. To obtain more reliable TWS changes from GNSS observations using SBF, more robust filtering algorithms should be introduced to separate the actual TWS changes from GNSS vertical displacements in future work.

Data Availability Statement
The