Surface Mass Variations from GPS and GRACE/GFO: A Case Study in Southwest China

: Surface mass variations inferred from the Global Positioning System (GPS), and observed by the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GFO) complement each other in terms of spatial and temporal coverage. This paper presents an analysis of regional surface mass variations inverted from GPS vertical displacements under different density distributions of GPS stations, and compares the GPS-derived mass variations with GRACE/GFO inversion results in spatial and temporal domains. To this end, GPS vertical displacement data from a total of 85 permanent GPS stations of the Crustal Movement Observation Network of China (CMONOC), the latest GRACE/GFO RL06 spherical harmonic (SH) solutions and GRACE RL06 mascon solutions are used to investigate surface mass variations in four regions or basins, including the Yunnan Province (YNP), Min River Basin (MRB), Jialing River Basin (JLRB), and Wu River Basin (WRB) in Southwest China. Our results showed that the spatial distributions and seasonal characteristics of GPS-derived mass change time series agree well with those from GRACE/GFO observations, especially in regions with relatively dense distributions of GPS stations (e.g., in the YNP and MRB), but there are still obvious discrepancies between the GPS and GRACE/GFO results. Scale factor methods (both basin-scaled and pixel-scaled) were employed to reduce the amplitude discrepancies between GPS and GRACE/GFO results. The results also showed that the one-year gap between the GRACE and GFO missions can be bridged by scaled GPS-derived mass change time series in the four studied regions, especially in the YNP and MRB regions (with relatively dense distributions of GPS stations). validation, X.L., J.C., and Q.L.; data curation, X.L. and T.L.; funding acquisition, B.Z.; writing—original draft preparation, X.L. and B.Z.; writing—review and editing, B.Z., J.C., and Q.L. All authors read and


Introduction
Surface mass variations mainly include terrestrial water storage (TWS) change, melting of polar ice sheets and mountain glaciers, atmospheric pressure change, ocean mass change, and mass changes associated with solid Earth geophysical phenomena. At present, large scale surface mass variations are primarily quantified by using the Gravity Recovery and Climate Experiment (GRACE) time-variable gravity field solutions and outputs from hydrological, atmospheric, and ocean models [1][2][3][4]. Even though GRACE satellite gravimetry has demonstrated the great potentials of measuring large-scale mass changes associated with land hydrological cycle, glacial melting and sea level changes, earthquakes, and human activities with unprecedented accuracy, it still has some major limitations: (1) The truncation of GRACE Spherical Harmonic (SH) coefficients and needed spatial filtering and smoothing to suppress the longitudinal stripping noises will introduce leakage bias and degrade the spatial resolution of GRACE-determined surface mass variations, which is limited to ~300-500 km [5][6][7]; (2) the uncertainty of GRACE low-degree SH coefficients is relatively larger, especially the degree-2 zonal term C20, which needs to be replaced with independent estimates from Satellite Laser Ranging (SLR) [8,9]; (3) the data quality of GRACE observations degrades after August 2016, which affects the accuracy of GRACE-inferred mass changes; (4) there is an approximately oneyear gap between the GRACE and GRACE Follow-On (GFO) missions, which affects the interpretation and applications of GRACE/GFO data.
Surface deformation due to loading caused by the redistribution of Earth's surface mass can be continuously monitored by the ground-based Global Positioning System (GPS) observation network. Since the surface loading inversion approach was proposed by Blewitt et al. [10], mass variations inverted from GPS surface deformation data have been serving as a complement to the GRACE/GFO observations both at global and regional scales [11][12][13][14][15][16][17]. Compared with GRACE/GFO, GPS has its own advantages: (1) The GPS network is widely distributed globally, and the observation data can be obtained in near real time at millimeter-level accuracy with higher spatio-temporal resolution; (2) the GPS surface deformation, especially the vertical displacement, is more sensitive to changes of local (~10 km) and regional mass loading, and the spatial resolution of mass changes inverted from a dense GPS network can reach ~50 km [16,17]; (3) due to the continuous and real-time monitoring of the surface deformation by GPS networks, GPS-derived mass variations can help bridge the gap between the GRACE and GFO missions [15,16].
Because Earth's surface deformation due to surface mass loading can also be estimated by GRACE SH gravity solutions [13], numerous previous studies have focused on the consistency analysis of surface deformations observed by GPS and GRACE. Davis et al. [18] demonstrated the high correlation of annual variations of GPS-measured and GRACE-derived vertical displacement time series in the Amazon Basin. Since then, several other studies had shown good consistency between GPS-measured and GRACE-derived vertical displacements, like in the Nepal Himalayas [19], West Africa [20], Greenland [21], Southeastern Poland [22], the Southeastern Tibetan Plateau [23], Southwest China [24], and the North China Plain [25]. Moreover, Wahr et al. [26] found that the GPS vertical displacement component combined with the horizontal displacement component could be used to determine whether nearby loading was concentrated in a small region, and where that region was. However, due to the different measurement modes and error sources of GPS and GRACE, poor agreements between GPS and GRACE vertical displacement estimates have also been reported in Central American [27] and Europe [28]. Nevertheless, these previous studies suggested that surface deformations determined by GPS and GRACE are consistent in most of the studied regions, and the surface mass variations can be inverted by GPS deformation data when non-loading effects are carefully removed.
A number of previous studies have also demonstrated that the regional surface mass variations can be independently inverted by using continuous GPS surface deformation data. Argus et al. [16] used GPS vertical displacement data to invert seasonal variations of TWS, and found that the inversion results showed good consistency with those from GRACE and hydrological models. The study further showed that TWS change can be recovered by GPS vertical displacements at a spatial resolution of ~50 km over California. Fu et al. [17] inverted monthly surface mass change time series using GPS vertical displacement data in the Washington and Oregon regions, and showed that GPSinferred water mass variations were highly correlated and consistent with those derived from GRACE and hydrological models, and could help fill in the gaps of GRACE. Jin and Zhang [29] estimated the mass changes using continuous GPS vertical displacement data in Southwestern USA, and the severe drought in 2012 was captured by the GPS-derived mass change time series, indicating that continuous GPS deformation observation has the potential to serve as a real-time drought indicator. Enzminger et al. [30] employed GPS vertical displacements to assess the accuracy of the snow water equivalent over Western USA Mountains. Zhang et al. [31] utilized GPS vertical displacements from 29 GPS stations to study mass variations at a spatial resolution of approximately 100 km over Yunnan Province in China. More recently, Ferreira et al. [32] presented a resolution analysis of daily mass variations inverted by GPS vertical displacement data over South America. Fok et al. [33] discussed a method for improving the inversion accuracy of GPS-derived mass variations in Southwestern China through combined applications of terrain-corrected GPS vertical displacements and supplementary GRACE spatial data constraints.
This study aims to further analyze the performance of regional mass variations inverted by GPS vertical displacement data under different density distributions of GPS stations, and examine the consistency of mass changes inferred from GPS vertical displacements, GRACE, and the latest GFO observations in the spatial and temporal domains, and demonstrate the feasibility of GPS-derived mass variations to fill in the gap between GRACE and GFO missions. To this end, continuous GPS vertical displacement time series from the Crustal Movement Observation Network of China (CMONOC), the latest GRACE/GFO RL06 SH solutions, and GRACE RL06 mascon solutions from the Center for Space Research (CSR) were used to investigate surface mass variations in four regions of Southwest China. The spatial distributions and seasonal characteristics of the mass change time series derived from GPS vertical displacements and GRACE/GFO were analyzed and compared in detail. The scale factor methods (both basin-scaled and pixel-scaled) used in Landerer et al. [34] were employed to reduce the amplitude discrepancies between GPS-derived and GRACE/GFO-inferred mass change time series. The rest of the paper is organized as the following. Overview of the studied area and the datasets used for inferring surface mass variations are described in Section 2. Data postprocessing of GPS vertical displacements, and inversion model for surface mass variations are presented in Section 3. In Section 4, comparisons of different mass change inversion results and discussions of differences among these results are presented, and the potentials of using GPS-derived mass change time series to bridge the GRACE/GFO mission gap are also demonstrated. Finally, conclusions are drawn in Section 5.

Study Area
In this paper, Southwest China was chosen as the study area, which includes four regions or basins ( Figure 1): The Yunnan Province (YNP), Min River Basin (MRB), Jialing River Basin (JLRB) and Wu River Basin (WRB). Table 1 lists the areas, number of GPS stations, and station density for each of the regions. Rivers and valleys are widely distributed over this studied area, and the landform mainly consists of plateaus and mountains with great relief. Furthermore, it belongs to subtropical monsoon climate zone, with uneven annual average temperature distribution and abundant rainfall. Hence, it is of great significance for the management and utilization of natural resources to investigate surface mass variations in Southwest China.

GPS Vertical Displacement Data
GPS vertical displacement data during the period of January 2011 to August 2019 from a total of 85 permanent GPS stations of CMONOC (ftp://ftp.cgps.ac.cn/products/position/gamit/) were used to invert the surface mass variations over Southwest China. All the permanent GPS stations of CMONOC are equipped with multi-frequency receivers with a sampling interval of 30 s, and the GPS station coordinates time series are processed by the GAMIT/GLOBK software (release 10.4) based on the double-differenced model [25,35]. Firstly, loosely constrained daily coordinates and covariance were processed using GAMIT software. In this procedure, satellite orbit parameters and satellite antenna offsets were adjusted using the reprocessed orbits from the International GNSS Service (IGS). Solid Earth tides and pole tides were corrected according to the International Earth Rotation Service (IERS) conventions 2003, and ocean tidal loading was corrected using the Finite Element Solutions 2004 (FES2004) model. Receiver antenna phase center was corrected using the absolute antenna phase center model [36]. Tropospheric delays were modeled using Global Mapping Function (GMF) [37] every hour, and the first order ionospheric delays were eliminated through the ionosphere-free linear combination. Atmospheric zenith delay and atmospheric gradients were estimated every hour and every day, respectively. Secondly, the loosely constrained daily solutions with station coordinates, covariance, atmosphere delays, and orbit parameters were combined with global solutions produced by the Scripps Orbital and Position Analysis Center (SOPAC, http://sopac.ucsd.edu/) using the GLOBK software. Finally, the datum transformation was performed to obtain a station coordinate time series with respect to the International Terrestrial Reference Frame (ITRF) 2008 using approximately 50 globally distributed IGS stations through seven-parameter transformation model, and then the daily GPS vertical displacement time series were obtained. More details of GPS data processing can be found from the manual of CMONOC reference stations (ftp://ftp.cgps.ac.cn/doc/processing_manual.pdf).

GRACE/GFO-Inferred Surface Mass Variations
GRACE and GFO monthly SH solutions (up to degree and order 60) were used to infer surface mass variations in accordance with the approach presented by Wahr et al. [38]. All of the SH solutions used in our study are published by the CSR at the University of Texas at Austin (UTCSR), including CSR RL06 GRACE SH solutions (hereinafter called CSR GRACE-SH) covering the period August 2002 to June 2017 (ftp://isdcftp.gfz-potsdam.de/grace/Level-2/CSR/RL06/) and the latest CSR RL06 GFO solutions (hereinafter called CSR GFO-SH) from June 2018 to August 2019 (ftp://isdcftp.gfzpotsdam.de/grace-fo/Level-2/CSR/RL06/). Since the coordinate origin of the GRACE/GFO data system is defined as the center of mass of the total Earth (CM) frame, the degree-1 SH terms are not included in GRACE/GFO SH solutions. However, GPS surface deformations are observed with respect to the center of surface figure (CF) frame, and the degree-1 SH terms reflect the motion of the geocenter with respect to the CM frame. In order to keep the consistency of the reference frame between GPS and GRACE/GFO, the degree-1 SH coefficients had been added back to the GRACE/GFO SH solutions using the degree-1 components provided by the GRACE project as supplementary datasets (GRACE Technical Note 13) [39]. Meanwhile, due to the uncertainty of lowdegree SH coefficients from GRACE/GFO data, the C20 and C30 coefficients had been replaced by SLR solutions provided by the GRACE/GFO Science Data System (SDS) team [40].
Due to SH coefficients at high degrees and orders being contaminated by the longitudinal stripes and other errors, the de-striping method P3M6 [41] and Gaussian filtering with the averaging radius of 300 km were used to suppress to spatial noises. But the above post-processing of SH solutions will reduce the spatial resolution and lead to signal leakage. In this study, the leakage errors were corrected by the basin scale factor method [34], and the scaling factors were derived from the Global Land Data Assimilation System (GLDAS) Noah outputs [42] by calculating the ratio between the original TWS and the truncated and filtered TWS in each studied region. Moreover, in order to be consistent with the inversion results from GPS vertical displacement data, the atmosphere combined with ocean (GAC) fields from Atmosphere and Ocean De-aliasing Level-1B (AOD1B) RL06 products (ftp://isdcftp.gfz-potsdam.de/grace/Level-2/CSR/RL06/) had been added back to the SH solutions. By doing this, both of the inversion results from GPS and GRACE/GFO contained the mass changes from TWS, non-tidal atmosphere, and ocean.
Due to the lack of adequate in situ measurements to validate GRACE/GFO observations, it is challenging to assess uncertainty of GRACE/GFO-inferred mass changes. Generally, the Root Mean Square (RMS) of mass variations over an open ocean area covering the same latitude range of the four studied regions can be approximately regarded as the uncertainty in GRACE/GFO mass change estimates [43,44]. In this study, an ocean area in the Pacific [20°N-35°N,170°E-220°E] was picked up to estimate the uncertainty of GRACE/GFO-inferred mass variations over Southwest China. It is approximately at the same latitudes as Southwest China and far away from land to avoid leakage of mass variation from land hydrology. It should be noted that this uncertainty assessment does not account for leakage errors in the four studied regions, which would lead to underestimate the uncertainty of GRACE/GFO mass change estimates.

GRACE Mascon Solutions
GRACE mascon solutions are typically presented as equivalent surface mass variations, and the longitudinal stripes are suppressed by applying regularization constraints [45][46][47]. Therefore, GRACE mascon solutions do not require de-striping or spatial filtering like in the SH solutions, and they can offer improved spatial resolution for analysis of regional surface mass variations. For these reasons, CSR GRACE RL06 mascon solutions (hereinafter called CSR GRACE-M, downloaded at http://www2.csr.utexas.edu/grace/RL06_mascons.html) were used to compare with the inversion results obtained from GPS vertical displacements and GRACE SH solutions.
The CSR mascon solutions were estimated on an equal area geodesic grid consisting of hexagonal tiles approximately 120 km wide (~1° × 1° at the equator), and mass variations in each mascon cell were computed from GRACE range-rate observations via partials with respect to SH expansion. Constraints on the solution include application of a time-variable regularization matrix from 200 km Gaussian smoothed regularized SH solutions, separation of land and ocean signals to reduce leakage, and forward modeling of annual and trend signals from ice losses in polar regions and mountain glaciers to avoid over-constraining the mascon estimation [47,48]. The CSR mascon solutions were provided on 0.5 × 0.5° grids, but the actual spatial resolution was significantly lower than that of grid resolution. Therefore, direct estimates of mass changes from CSR mascon solutions were also subjective to potential leakage bias [3,49]. Like the GRACE SH solutions, the leakage corrections for CSR mascon solutions were also reduced by the basin scale factor method, but a 200 km Gaussian filtering was applied to calculate the scaling factors based on the GLDAS Noah outputs. The uncertainty assessment of CSR mascon solutions over Southwest China is the same as that mentioned in Subsection 2.3.

Post-processing of GPS Vertical Displacements
Because the effects of solid Earth tides, ocean tides loadings, and pole tides have been removed from the GPS vertical displacement time series of CMONOC, the GPS vertical displacements mainly reflect surface deformation caused by mass changes from TWS, non-tidal atmosphere, and ocean. However, the draconitic errors (i.e., ~352 days [50]) that are contained in the GPS vertical displacements must be removed before estimating surface mass variations. Unlike Glaser et al. [51], who used the combination of GPS and SLR data to remove the draconitic errors, we followed the method used in Fok et al. [33] to reduce the influences of draconitic errors on the desired seasonal variation (annual and semi-annual) signals by applying spectral filtering in the frequency domain. In addition, due to the influences of anthropogenic activities and natural environments, there would be data jumps and discontinuities in the GPS vertical displacement time series, e.g., the HBZG station shown in Figure 2. It is essential to carry out outlier detection and data interpolation for these vertical displacement time series. Here, the gross error detection method from Nikolaidis et al. [52] was used to remove the outliers, and the cubic spline interpolation and Principal Component Analysis (PCA) methods [53] were used to complement the missing data. Then the GPS vertical displacement time series of each station was fitted with a position, velocity, and sinusoids with semi-annual and annual periods. Because the elastic long-term trend (due to glacial isostatic adjustment or tectonic effects) can hardly be separated from the GPS vertical displacement time series, and mass changes in the studied regions are dominated by seasonal variability, we deducted it from GPS vertical displacement time series and focused on annual and semi-annual variations. Furthermore, in order to be consistent with the temporal resolution of GRACE/GFO monthly gravity field solutions, the daily vertical displacement observations were converted into monthly averages before inverting surface mass variations. Figure 2 shows the detrended GPS vertical displacement time series (blue line) and their fitted seasonal time series (yellow line) for three selected permanent GPS stations (YNSM, HBZG, and GZGY) in Southwest China, and the monthly variations of vertical displacements derived from CSR GRACE-SH (red line) and CSR GFO-SH (green line) are also presented. Atmospheric and ocean no-tidal loading effects were added back to GRACE/GFO estimates using the GAC fields from AOD1B RL06 products. As shown in Figure 2, for some stations (e.g., YNSM and HBZG), GPS-measured and GRACE/GFO-derived vertical displacement time series show good consistency with each other, while for some other stations (e.g., GZGY), their consistency is poor, perhaps due to technical errors in GPS data processing.  Figure 3 shows the correlation coefficient between GPS-measured and GRACE-derived vertical displacement time series for each GPS station over the period January 2011 to June 2017. Here, the GRACE-derived vertical displacements were inferred from CSR GRACE-SH. The dots represent the distribution of permanent GPS stations, and colors of the dots represent the correlation coefficient values. In addition, the locations of YNSM, HBZG, GZGY, and KMIN station are also marked. Among the 85 permanent GPS stations over Southwest China, 71% of them had a correlation coefficient greater than 0.5, and 50% had a correlation coefficient greater than 0.6. In addition, most of the GPS stations with good correlations are located in the YNP and MRB regions. However, the correlations of individual GPS stations are obviously poor, like GZGY station in WRB and KMIN station in YNP. In addition to technical errors in GPS data processing, another important reason that accounts for the poor correlation between GPS-measured and GRACE/GFO-inferred vertical displacements is that GPS deformations are more sensitive to local mass loading. One of the criteria for identifying GPS abnormal stations is to check whether there is any obvious difference compared with the nearby stations. However, due to the relatively sparse distributions of GPS stations in the four studied regions, it is difficult to distinguish whether or not the abnormal station is affected by local effects. In our study, the GPS station data checked with poor correlation and consistency would be discreetly used based on the influence analysis of a single GPS station (e.g., GZGY and KMIN) on the inversion of mass changes over Southwest China.

Inversion Model for Surface Mass Variations
On the basis of elastic load theory, the Earth's deformation caused by the redistribution of surface mass loading can be approximately expressed as elastic deformation, and elastic deformation induced by the mass variation can be expressed as the Green's Function [54,55]. It is assumed that a uniform mass disk with an angle radius of α is placed on the surface of the Earth with spherical symmetry, self-gravitation, elasticity, and isotropy. The total mass of the disk is as follows [26]: in which w h and w ρ are the equivalent water height (EWH) and density of water, respectively, and R represents the radius of the Earth. The vertical displacement u induced by a mass loading can be expressed as [29]: where e M is the mass of the Earth, θ is the angular distance, l h are the load Love numbers from Wang et al. [56], and l P are the Legendre polynomials of order l.
In order to better investigate surface mass variations, the whole studied area of Southwest China was divided into uniform square grid cells with a spatial resolution of 0.25 × 0.25°. Then, by substituting Equation (1) into Equation (2), the unknown EWH vector for the above grid cells can be solved based on least-squares estimation. However, the correlation of the GPS vertical displacement observations will result in a discrete ill-posed problem for inverting regional surface mass variations. To deal with the ill-posed problem, the Tikhonov regularization method was applied as follows [16]: where A is the design matrix decided by Green's Function, X is the unknown EWH vector that needs to be solved, u is the observation vector of GPS vertical displacements, σ is the noise standard deviation of the observations, L is the regularized Laplacian matrix [57] which can constrain the mass changes more stably between the adjacent grid cells, and λ is the regularization parameter which can be estimated by the generalized cross validation (GCV) method [58,59]. Finally, the EWH vector can be estimated as follows: ( ) The uncertainty of surface mass variations inverted from GPS vertical displacements can be acquired from covariance analysis of the estimated EWH vector. The covariance matrix can be obtained by inverting the normal equation matrix, and the diagonal elements of covariance matrix are the error variances of estimated EWH.

Spatial Distributions of Surface Mass Variations Derived from GPS and GRACE/GFO
Both GPS-measured and GRACE/GFO-derived vertical displacement time series showed apparent seasonal variations over Southwest China, and these two types of vertical displacement time series showed good correlation and consistency except for some individual stations, e.g., GZGY station (see Figures 2 and 3), and KMIN station (see Figure 3). As mentioned earlier, the poor correlation between GPS-measured and GRACE/GFO-inferred vertical displacements may have been caused by technical errors in GPS data processing or local effects on GPS station. In order to analyze the influence of a single GPS station (i.e., GZGY and KMIN) on the inversion of mass variations over Southwest China, a closed-loop simulation was conducted based on the simulated GPS vertical displacements. Using the locations of 85 permanent GPS stations of CMONOC, the simulated GPS vertical displacement data were generated from CSR GRACE-SH plus GAC fields in March and November 2011, respectively. For these two months, surface mass variations derived from the above input signals were different from each other, and the signals in March 2011 were much weaker than those in November 2011. For instance, the maximum, minimum, mean, and standard deviation of  Firstly, compared with KMIN station, there are no other stations around GZGY within ~200 km, and there is no constraint information available nearby. Secondly, there are only two GPS stations located in WRB region, and GZGY station may provide a representative sample of WRB, so whether GZGY station data was used or not would make an obvious influence on the inversion results, especially in the months with relatively larger signal variations (e.g., in November 2011). However, these differences caused by GZGY or KMIN station are very small compared to the mass change signals themselves (see Figure 4). Although the poor correlations between GPS-measured and GRACE/GFOinferred vertical displacements were found at GZGY and KMIN stations, we still included GZGY and KMIN station data to invert surface mass variations over Southwest China based on the above analysis. In order to examine the performance of GPS and GRACE/GFO for monitoring surface mass variations, we employed the GPS vertical displacement data from a total of 85 permanent GPS stations of CMONOC to invert mass variations over Southwest China, and compared them with GRACE/GFO observations. Figure 6 shows the spatial distributions of surface mass variations obtained from GPS vertical displacements and CSR RL06 GRACE/GFO SH solutions in February and July 2011, and January and July 2019, respectively, and the red dots in Figure 6a,d represent the location of GPS stations.
As presented in Figure 6, in February 2011, surface mass variations derived from CSR GRACE-SH show mainly negative anomalies in the YNP region, the northwestern of MRB, and JLRB, while indicating positive anomalies in WRB, the southeast of MRB, and JLRB. The spatial distribution of mass variations inverted by GPS vertical displacements is comparable to that observed by GRACE, but the positive anomalies in WRB, the southeast of MRB, and JLRB were smaller than those of CSR GRACE-SH estimates. In July 2011, the spatial distributions of the mass variations obtained by GPS vertical displacements and CSR GRACE-SH were almost the same. Among them, positive anomalies were mainly presented in the YNP and MRB regions, while negative anomalies appeared mostly in the JLRB and WRB regions. In January 2019, mass variations derived from CSR GFO-SH showed positive anomalies in nearly the entire YNP and WRB and southeast of JLRB and MRB, while the northwest of JLRB and MRB were dominated by negative anomalies. The spatial distribution of GPSinferred mass variations was comparable to that of the CSR GFO-SH result, but the GPS magnitudes of negative anomalies in the northwest of JLRB and MRB were larger than those from CSR GFO-SH. The distribution of mass variations in July 2019 almost showed an inverse pattern compared with that of January 2019. In general, for most months (not all of them are shown in this paper), spatial distributions of mass variations obtained by GPS and GRACE/GFO showed good consistency. These results indicate that GPS can help independently monitor surface mass variations in our studied regions, and mass changes retrieved by GPS vertical displacements can be verified by GRACE/GFO results.

Seasonal Surface Mass Variations Revealed by GPS and GRACE/GFO
In order to better understand surface mass variations over Southwest China, seasonal variations of the inverted results from GPS and GRACE/GFO were further analyzed in the four regions, i.e., the YNP, MRB, JLRB, and WRB, respectively. The annual amplitudes of the mass variations estimated by GPS vertical displacements, CSR GRACE-SH, and CSR GRACE-M from January 2011 to June 2017 are presented in the 3 panels of Figure 7, respectively. As shown in Figure 7, the largest annual amplitudes are located in southwest YNP, while the northeast YNP and the other three regions show significantly smaller annual amplitudes. Furthermore, the distribution of annual amplitude obtained by GPS was in good agreement with those of CSR GRACE-SH and CSR GRACE-M. This suggests that the GPS-inferred mass changes are fairly reliable at seasonal time scales, based on comparisons with different GRACE results.  As shown in Figure 8a, YNP has the largest area with the densest distribution of GPS stations (Table 1) For the JLRB region shown in Figure 8c, GPS-derived surface mass changes agreed well with GRACE/GFO estimates. However, due to the sparse and uneven distribution of GPS stations in the JLRB region, there was a small annual phase lag between GPS and GRACE/GFO results. Figure 8d shows that, from 2011 to 2014, there were fairly obvious differences between GPS and GRACE/GFO derived mass changes in the WRB region. One reason is to think that the consistency of GPSmeasured and GRACE-derived vertical displacements for the GZGY station located in the WRB was poor during this period, as shown in Figure 2. Because there are no other GPS stations around GZGY station and no constraints available nearby, it is difficult to determine whether the disagreement was caused by GPS technical errors or local effects. From the simulation in Subsection 4.1, we found that a single GPS station (i.e., GZGY) has little influence on the inversion of mass variations over Southwest China. We think the main reason is that WRB has the smallest area which is close to the limit of basin size that GRACE can observe, and the lowest station density (see Table 1) in the WRB region leads to larger uncertainty in the inverted surface mass variations.  Table 2 lists the amplitudes and phases of annual and semi-annual mass change time series for the four regions based on GPS and GRACE inversion results from January 2011 to June 2017, and the RMS differences of the mass change time series between GPS and GRACE (CSR GRACE-SH, CSR GRACE-M) are also presented. The semi-annual amplitudes were much smaller compared to the annual amplitudes in the four studied regions, and uncertainties at semi-annual time scale were expected to be relatively larger. For the YNP region, both annual amplitudes and phases agreed well between GPS and GRACE inversion results. The annual amplitudes from GRACE were slightly larger than that from GPS, and the annual amplitude from CSR GRACE-M agreed better with GPS estimates. However, the semi-annual amplitude of GPS estimates was notably larger than the two GRACE results (CSR GRACE-SH, CSR GRACE-M). Moreover, the RMS differences of mass change time series between GPS and CSR GRACE-SH, and GPS and CSR GRACE-M were ~24.61 and 26.11 Gt, respectively, which indicates that there are still obvious discrepancies between GPS and GRACE inversion results.  For the MRB region listed in Table 2, annual amplitudes from GRACE (CSR GRACE-SH, CSR GRACE-M) were smaller than that from GPS, while the GPS semi-annual amplitude was smaller than those of GRACE. Meanwhile, the RMS differences of mass change time series between GPS and CSR GRACE-SH, GPS and CSR GRACE-M were ~8.74 and 7.48 Gt, respectively. For the JLRB, annual amplitudes of GPS and GRACE estimates were in good agreement, but the CSR GRACE-M estimates showed slightly larger annual amplitude. Additionally, the semi-annual amplitudes of the three estimates were almost the same, while the GPS series showed a 3-4 month phase lag compared with GRACE series. The RMS differences of the mass change time series between GPS and CSR GRACE-SH, GPS, and CSR GRACE-M were ~7.82 and 6.33 Gt, respectively. As for the WRB, the annual amplitude of mass changes inferred from GPS was underestimated, while the semi-annual amplitude from GPS was larger than those from the two GRACE inversion results. The RMS differences of the mass change time series between GPS and CSR GRACE-SH, GPS, and CSR GRACE-M were ~4.42 and 3.96 Gt, respectively.
In order to verify the reliability of mass variations derived from GPS and GRACE in the four studied regions, we further estimated the uncertainties from these two types of data. As mentioned earlier, the uncertainty of surface mass variations inverted from GPS vertical displacements was acquired from covariance analysis of the estimated EWH. The error RMS of mass changes derived from monthly covariance matrices for the YNP, MRB, JLRB, and WRB regions were 8.58, 3.61, 3.50, and 1.91 Gt, respectively. The uncertainty in GRACE estimates were obtained from the RMS of mass variations over an ocean area that covered the same latitudinal range of the four studied regions. The error RMS derived from CSR GRACE-SH and CSR GRACE-M for the YNP, MRB, JLRB, and WRB regions were 3.58 vs. 4 Table 2, all the above error RMS in the four studied regions are smaller. However, if the area of the studied region is smaller, the uncertainty of the corresponding estimation results will be larger. Especially for the WRB, the error RMS of differences of the mass changes between GPS and GRACE (2.07 or 2.14 Gt) were close to their RMS differences (4.42 or 3.96 Gt). As mentioned before, one reason is that WRB has the smallest area which is close to the limit of basin size that GRACE can observe (although the leakage was corrected in this study). Another reason is that WRB has the lowest station density compared to the other three regions, which leads to larger uncertainty in GPS estimates.
To sum up, the good consistencies between GPS-and GRACE/GFO-derived surface mass change time series show that GPS could well detect the characteristics of seasonal variations over Southwest China (especially in YNP and MRB) from January 2011 to August 2019, and the GPS-derived mass change time series provided an alternative way to fill in the gap between the GRACE and GFO missions. In addition, there are still discrepancies between GPS and GRACE/GFO inversion results. The reasons for this are twofold. Firstly, GRACE/GFO estimates mainly reflect large-scale mass variations with spatial resolution greater than 300km and strongly constrain the total mass variations, while the GPS principally detects regional surface mass variation at a local point (thus constraining the distribution of surface mass), and relatively sparse GPS stations cannot fully infer the details of surface mass variations in the whole studied area. Secondly, differences of error sources contained in GPS vertical displacements and GRACE/GFO data, and different processing strategies (e.g., smooth filtering and constraint solving techniques) can also induce discrepancies between the inversion results of these two types of data.

GPS-Derived Surface Mass Variations to Bridge GRACE/GFO Mission Gap
As discussed above, surface mass change time series inverted by GPS and GRACE/GFO over Southwest China agree well, especially in the YNP and MRB regions. Because GPS can continuously monitor surface mass changes on regional scales, GPS-derived mass variations can serve as a supplementary data to bridge the GRACE/GFO mission gap. To minimize potential biases between GPS-and GRACE/GFO-derived mass changes due to differences in data processing procedures, we employed scale factor methods (both basin-scaled and pixel-scaled) adopted in Landerer et al. [34] to reduce the amplitude discrepancies between GPS-and GRACE/GFO-derived mass change time series. Specifically, we took CSR GRACE-SH-derived mass change time series from January 2011 to June 2017 as a reference, and estimated basin scale and pixel scale factors by minimizing the misfits between GPS-derived and GRACE-inferred surface mass change time series based on drainage basin (basin-scaled) and each grid cell (pixel-scaled), respectively. The feasibility of these scaling factors was tested by comparison of the scaled GPS-derived and GFO-inferred mass change time series during the period from June 2018 to August 2019.
The four panels of Figure 9 show  As shown in Figure 9, in the YNP and MBR regions, all of the surface mass change time series inverted by GPS and GFO show good consistency, and GPS-derived mass change time series scaled by basin-scaled and pixel-scaled methods are almost the same. The annual amplitudes of mass change time series obtained by basin-scaled and pixel-scaled methods were slightly smaller than those of the raw surface mass change time series inverted by GPS, but were much closer to the CSR GFO-SH results. As for JLRB and WRB, surface mass change time series inverted by GPS and GFO showed some phase differences, and the annual amplitudes of basin-scaled and pixel-scaled GPS mass changes were smaller than those of the CSR GFO-SH results. The discrepancies shown in the latter two regions demonstrate that the relatively sparse and uneven distributions of GPS stations limit the accuracy of GPS-inferred regional mass variations, and the limited spatial resolution of GFO also leads to larger uncertainty in mass change estimates. Table 3 lists the annual amplitudes and phases for the four regions (i.e., the YNP, MRB, JLRB, and WRB) based on the raw GPS-derived (Raw GPS) mass change time series, basin-scaled (basinscaled GPS), and pixel-scaled (pixel-scaled GPS) mass change time series, and GFO-observed (CSR GFO-SH) mass change time series from June 2018 to August 2019, and the RMS of mass change time series minus seasonal variations for CSR GFO-SH and those three GPS inversion results are also given for comparisons.  As shown in Table 3, in the YNP region, the annual amplitude of the raw GPS-derived mass change time series is much larger than that of GFO observations (72.69 vs. 59.19 Gt), but becomes substantially smaller after the scaling factor corrections. Moreover, the annual amplitude of basinscaled GPS mass change time series agreed remarkably well with GFO observations (59.19 vs. 59.84 Gt), and the RMS residual also reduced correspondingly after the basin-scaled factor corrections (16.71 vs. 13.61 Gt). However, the RMS residual increased by using the pixel-scaled factor (16.71 vs. 18.04 Gt). For the MRB, the annual amplitude of raw GPS mass change time series was also significantly larger than that of GFO observations (22.42 vs. 12.53 Gt), and after scaling factor corrections, GPS annual amplitudes agreed notably better with GFO observations (10.03 or 14.71 vs. 12.53 Gt). The RMS residuals were significantly reduced (5.50 vs. 2.49 or 3.46 Gt). However, in the JLRB and WRB regions, neither basin-scaled nor pixel-scaled GPS mass change time series showed improved agreements with GFO observations, indicating the challenges of inverting regional mass changes when inadequate numbers of GPS stations are available. Certainly, due to the limited spatial resolution of GFO, the larger uncertainty of GFO estimates would fairly affect the comparisons. This may be evident in the JLRB region, where GFO observations did not show a clear and typical seasonal signal as in the other three regions. Nevertheless, the RMS residuals did show obvious decreases in the WRB region after the two scaled factors corrections (2.58 vs. 1.51 or 1.87 Gt). In addition, the RMS residuals of basin-scaled GPS time series were smaller than those of the pixel-scaled GPS time series in the four studied regions, which shows that basin-scaled factor method is more suitable to scale the GPS-derived mass change time series than the pixel-scaled factor method in our study.
In order to demonstrate the scaled GPS-derived surface mass variations can help bridge the GRACE/GFO mission gap, Figure 10 Figure 10, either basin-scaled or pixel-scaled mass change time series appear to well bridge the GRACE/GFO mission gap in these four regions, especially in YNP and MRB. Besides, due to the relatively dense distributions of GPS stations distributed in the YNP and MRB regions, seasonal variations during the one-year mission gap were also well captured. Seasonal variations in JLRB and WRB were significantly smaller than those in YNP and MRB, which also played a role in the poorer agreements between GPS and GFO results in the JLRB and WRB regions. As can be seen from the above analysis, the annual amplitudes of GPS-derived mass change time series obtained from relatively dense distributions of GPS stations in the YNP and MRB are usually greater than that of GFO-inferred mass change time series. For sparsely uneven distributions of GPS stations in the JLRB and WRB, the annual amplitudes of GPS-derived mass change time series were smaller than that of GFO-inferred mass change time series. However, the scaled GPS-derived mass change time series can bridge the one-year mission gap between GRACE and GFO in these four regions, especially in the YNP and MRB where GPS stations are relatively densely distributed. From the RMS residuals listed in Table 3, we believe the basin-scaled factor method is more suitable to scale the GPS-derived mass change time series than the pixel-scaled factor method in our study.

Conclusions
In this study, GPS vertical displacement data from a total of 85 permanent GPS stations of CMONOC, the latest GRACE/GFO RL06 SH solutions, and GRACE RL06 mascon solutions were used to investigate surface mass variations over Southwest China. Spatial distributions and seasonal characteristics of surface mass variations inverted by GPS and GRACE/GFO were analyzed and compared in detail in the four studied regions (i.e., the YNP, MRB, JLRB, and WRB) in Southwest China. Moreover, the potential of using GPS-derived mass change time series to bridge the GRACE/GFO mission gap was also demonstrated.
For most months, spatial distributions of GPS and GRACE/GFO derived surface mass variations showed good consistency during the period of January 2011 to August 2019. Annual amplitudes and phases of GPS-derived surface mass changes agreed well with those of GRACE estimates (e.g., CSR GRACE-SH and CSR GRACE-M). Seasonal characteristics of surface mass changes over the period January 2011 to August 2019 were well captured by GPS inversions in each of the four regions, and GPS-derived mass change time series showed good agreements with GRACE/GFO observations, especially in the YNP and MRB regions where GPS stations are relatively densely distributed. However, evident discrepancies between GPS and GRACE/GFO derived mass changes still exist (e.g., in JLRB and WRB). One particular reason is that the relatively sparse GPS stations in JLRB and WRB are not adequate for inverting the details of regional mass variations. Another major reason is that both JLRB and WRB have smaller areas which are close to the limit of basin size that GRACE/GFO can observe (although the leakage was corrected), and the uncertainty of GRACE/GFO estimates would largely affect the comparisons in these cases. In addition, the error sources contained in GPS vertical displacement measurements and GRACE/GFO data, as well as different data processing strategies, will also induce discrepancies between GPS and GRACE/GFO derived mass change results.
Scale factor methods (both basin-scaled and pixel-scaled) were employed to reduce the amplitude discrepancies between the inversion results from GPS and GRACE/GFO. We showed that the scaled GPS-derived mass change time series can help bridge the gap between the GRACE and GFO missions, and the seasonal variations were accurately captured in the YNP and MRB regions by the scaled GPS results during the one-year mission gap. Therefore, GPS-derived surface mass variations can complement the data gaps of GRACE/GFO to obtain a long-term and continuous mass change time series at regional scales when adequate numbers of GPS stations are available. Moreover, the basin-scaled factor method was verified to be more suitable than the pixel-scaled factor method to scale the GPS-derived mass changes through the RMS residuals analysis. Future work will benefit from reprocessed GPS vertical displacement data and increased spatial coverage of GPS stations. In addition, as two independent geodetic techniques, joint inversions of surface mass variations from GPS and GRACE/GFO data may help increase the spatial and temporal resolution of GRACE/GFOderived surface mass changes.