Improved Estimation of Regional Surface Mass Variations from GRACE Intersatellite Geopotential Differences Using a Priori Constraints

We presented an improved method for estimation of regional surface mass variations from the Gravity Recovery and Climate Experiment (GRACE)-derived precise intersatellite geopotential differences using a priori constraints. An alternative analytic formula was proposed to incorporate the K-band ranging (KBR) range rate into the improved energy balance equation, and precise geopotential differences were estimated from GRACE Level-1B data based on the remove-compute-restore (RCR) technique, which avoids the long-wavelength gravity signals being absorbed by empirical parameters. To reduce the ill condition for inversion of regional mass variations from geopotential differences, a priori information from hydrological models was used to construct the constraint equations, and the optimal regularization parameters were adaptively determined based on iterative least-squares estimation. To assess our improved method, a case study of regional mass variations’ inversion was carried out over South America on 2◦ × 2◦ grids at monthly intervals from January 2005 to December 2010. The results show that regional mascon solutions inverted from geopotential differences estimated by the RCR technique using hydrological models as a priori constraints can retain more signal energy and enhance regional mass variation inversion. The spatial distributions and annual amplitudes of geopotential difference-based regional mascon solutions agree well with the official GRACE mascon solutions, although notable differences exist in spatial patterns and trends, especially in small basins. In addition, our improved method can robustly estimate the mascon solutions, which are less affected by the a priori information. The results from the case study have clearly demonstrated the feasibility and effectiveness of the proposed method.


Introduction
Over the past 16 plus years, surface mass variations derived from the Gravity Recovery and Climate Experiment (GRACE) satellite gravimetry has significantly improved our ability to study the terrestrial water cycle, ice sheets and mountain glacier mass balance, sea level change, and ocean bottom pressure variations, as well as to understand responses to changes in the global climate system [1,2]. The traditional approach for modeling surface mass variations is based on the time-variable gravity field, which is represented by spherical harmonic (SH) basis functions [3,4], and surface mass variations obtained from GRACE Level-2 monthly SH solutions, e.g., the Center for Space Research (CSR), Geoforschungszentrum (GFZ), and Jet Propulsion Laboratory (JPL) published RL05/06 SH accurate than using the approximate correction terms presented by Jekeli [40]. However, intersatellite geopotential differences estimated from the above procedures may still contain systematic errors from the KBR measurements. The most common method to reduce the systematic errors is to absorb these errors by using empirical parameters (e.g., two cycle-per-revolution (CPR) parameters used in Han et al. [41] and Tangdamrongsub et al. [38], four CPR parameters used in Ramillien et al. [37]), but these empirical parameters would also weaken the temporal gravity signals, especially the long-wavelength components [37,39,43,44]. For example, several studies (e.g., Ramillien et al. [37], Frédéric Frappart et al. [43], and Ramillien et al. [44]) used the low degrees of SH solutions to complement the missing long-wavelength information of surface mass variations. Unlike the previous studies, we will employ the remove-compute-restore (RCR) technique to retain the temporal gravity signals to the full extent, and extract precise intersatellite geopotential differences based on an improved energy balance equation.
In the second step, the intersatellite geopotential differences are expressed as a function of surface mass variations, and are conventionally downward-continued from satellite altitude to the Earth's surface by Newton's formula of gravitational potential [34][35][36][37][38]. However, the noises in intersatellite geopotential differences will be amplified, and the estimation of mascon parameters is a well-known ill-posed problem. Therefore, regularization strategies have been employed to find numerically stable mascon solutions from intersatellite geopotential differences, either based on singular value decomposition (SVD) and L-curve analysis [36], or by introducing spatial constraints [37], and these regularization methods were verified from the inversion of surface mass variations over South America at 10-to 30-day intervals using different sizes of blocks. Besides, Han et al. [34] introduced a priori information equation to estimate the mascon parameters by an iterative least-squares estimation. Tangdamrongsub et al. [38] constructed space-domain and space-time covariance functions using a priori hydrological knowledge to obtain stable mascon solutions, in which the optimal regularization parameters were determined based on the L-curve criterion. Because the L-curve criterion is an empirical method to find the optimal regularization parameters, it is often a tedious job to solve for the mascon solutions with certain time intervals from a long period of GRACE data. In order to adaptively determine the optimal regularization parameters according to the data themselves, inspired by the work of Han et al. [34], we will make full use of a priori information to obtain the stable mascon solutions based on the iterative least-squares estimation.
Mascon solutions do not suffer from spectral truncation effects like in GRACE SH solutions, and can offer surface mass variations with higher spatio-temporal resolution at the regional scale [25,34,36,38]. Compared with the inversion of global mascon solutions, estimation of regional surface mass variations in the form of regional mascon solutions are computationally convenient because they can be estimated from a small subset of the global GRACE KBR data. In this study, we focused on the inversion of regional mascon solutions from GRACE intersatellite geopotential differences using a priori constraints. Firstly, we present an alternative analytic formula to incorporate the KBR range-rate into the energy balance equation, and estimate the precise in situ intersatellite geopotential differences from GRACE Level-1B data based on the RCR technique using the energy integral. Secondly, we employ an adaptive algorithm to determine the optimal regularization parameters based on the iterative least-squares method using hydrological model outputs as a priori constraints. Finally, a case study of regional mass variations' inversion over South America on 2 • × 2 • geographic grids at monthly intervals from January 2005 to December 2010 is presented, and our solutions are compared to the official GRACE mascon solutions to demonstrate the feasibility of the proposed methods.

Improved Energy Balance Equation
Due to the approximation of the potential rotation term in the energy integral developed by Jekeli [40], we follow a more accurate formulation of the potential rotation term presented by Remote Sens. 2020, 12, 2553 4 of 28 Guo et al. [42], and an improved energy balance equation can be formulated in the geocentric inertial frame for a single satellite as follows [39]: where V E is the Earth's gravitational potential; r and . r are the orbit position and velocity in the inertial frame; w is the angular velocity vector of the Earth-fixed frame with respect to the inertial frame; a = ∇V T + f includes the non-conservative force f, as well as the forces derived from high-frequency time-variable gravitational potential V T (due to N-body perturbation, solid Earth tides, ocean tides, pole tides, atmosphere and ocean signals, etc.); the integration in Equation (1) is from the initial time t 0 to t; and E 0 is a constant of integration.
The in situ intersatellite geopotential difference between GRACE's two satellites can be derived according to Equation (1) as: where the subscripts represent two GRACE satellites ('1', '2') and their difference ('12'). Because the energy integral in Equation (2) does not explicitly contain the KBR range-rate measurements, Shang et al. [39] introduced an alignment equation to update the relative velocity vector . r 12 by using the precise range-rate . ρ. Unlike the alignment equation, we simply decompose the relative velocity vector . r 12 into the line-of-sight (LOS) direction and its orthogonal direction, which is in the plane containing the relative position vector r 12 and relative velocity vector . r 12 , and the range-rate . ρ can be incorporated into the first term on the right side of Equation (2) (called kinetic energy term) as follows: where E k 12 is the difference of the kinetic energy between the two satellites, the range-rate . ρ = . r 12 ·e 12 , and e 12 = r 12 |r 12 | is the LOS unit vector. Our numerical test based on one month of GRACE data showed that the root-mean-square (RMS) difference between Equation (3) and the alignment equation [39] for estimation of E k 12 is less than 5.0 × 10 −7 m 2 /s 2 , which is negligible compared with the GRACE measurement accuracy.

Inversion Method of Regional Surface Mass Variations
The intersatellite geopotential difference V E 12 obtained from Equation (2) mainly includes static gravitational potential due to the solid Earth and time-variable part due to mass redistributions (including continental hydrology, postglacial rebound, earthquake, etc.). Here, we introduce a static reference geopotential U, and intersatellite geopotential difference due to mass changes, which is estimated from T 12 = V E 12 − U 12 , and U 12 is the geopotential difference computed by the reference geopotential U. Then, the residual geopotential difference T 12 can be used to invert regional surface mass variations. The study area can be divided into regular geographical grids, and the center of the jth block is located at the colatitudes θ j and the longitude λ j , and its surface area is δS j = R 2 ∆θ∆λ sin θ j , in which R is the mean Earth's radius, and ∆θ and ∆λ are the sampling angle intervals along the latitude and the longitude, respectively. For a given period ∆t = t − t 0 , the surface mass change of the j th block δm j can be expressed as: where P n is the Legendre polynomial of order n, and cosψ j 1,2 = cosθ j cosθ 1,2 + sinθ j sinθ 1,2 cos(λ j − λ 1,2 ), (r 1 , θ 1 , λ 1 ), and (r 2 , θ 2 , λ 2 ) are the geocentric radius, geocentric co-latitude, and longitude for each GRACE satellite ('1', '2'). Substituting Equation (4) and Equation (6) into Equation (5), and considering the load effects caused by the surface mass variations, we can obtain the observation equations to estimate regional mass variations from geopotential difference observables as follows: where k n are the load Love numbers that account for the elastic deformation potential for a viscoelastic Earth, and the values of load Love numbers were adopted from Wahr et al. [3]. It should be pointed out that in order to avoid the influence of edge effects on the regional mascon solutions, the boundary of calculation region should be expanded appropriately.

Adaptive Estimation Method Using a Priori Constraints
As it is known, recovering surface mass variations from geopotential differences based on Equation (7) is an ill-posed problem, because the noises will be amplified when geopotential difference data are downward continued from the satellite altitude to the Earth's surface. Several studies employed regularization techniques (e.g., SVD, L-curve) [36], spatial constraints [37] and a priori constraints [38] to reduce the ill condition in the regional mass variation solution from GRACE. However, it is a tedious job to find the optimal regularization parameters for massive GRACE data processing. In this study, we used a priori constraints based on geophysical information (e.g., hydrological models) to stabilize the ill-posed problem, and adaptively determined the optimal regularization parameters from the data themselves based on the iterative least-squares estimation method.
The linear observation equation from Equation (7) for estimating surface mass variations can be expressed as follows: where y is the observation vector of residual geopotential differences, A is the design matrix, x is the EWH parameter vector to be estimated, e is the residual vector with zero expectation and identity cofactor matrix, and σ 2 y is the error variance of geopotential differences. In order to stabilize the solutions, a stochastic model for the unknown parameter vector x is introduced based on the use of a priori information. Here, assuming that we have an a priori expectation x 0 and covariance matrix C x for the parameter vector x, the a priori information equations can be expressed as follows: where e 0 is the residual vector with zero expectation and covariance matrix C x , I x is the identity matrix associated to the parameter vector x.
The idea of regularization is to find the balance between the minimum norm of the residuals and that of the vector of unknown parameters relative to the a priori values. By minimizing Ax − y 2 + α x − x 0 2 = min, the optimal solution can be determined as: in which α is a regularization parameter to balance the weights between the observations and the a priori information. Because both the regularization parameter α and error variance σ 2 y are unknown in Equation (10), we can estimate them together through variance component estimation [45] as follows: in which n obs is the number of observations. As a result, we need to solve forx andασ 2 y iteratively using Equation (10) and Equation (11), respectively, starting with an initial value of ασ 2 y . In general, it is impossible to obtain the accurate a priori covariance matrix C x in Equation (10) and Equation (11), but it can be iterated starting from an initial covariance matrix. Here, we assume a stationary stochastic process for x, and the covariance matrix C x can be obtained from the empirical covariance function as follows: where σ 2 x and β are the variance and correlation distance of the covariance function, and d ij is the distance between the centers of the ith block and the jth block. The empirical values of the covariance function are computed from the gridded values (e.g., a priori hydrological models) in the studied area as follows [38]: where f i and f j are two surface mass variation values on two grid blocks separated by the distance d ij , N is the number of such pairs, and ∆d is an interval range (e.g., with a 2 • grid interval). Then,σ 2 x andβ can be estimated after linearization of Equation (12) based on the least-squares estimation.
Therefore, during the iterative process of solving forx andασ 2 y , we need to add one more step for C x in Equation (12), and the empirical values of covariance in Equation (13) for each iteration are computed from the intermediate estimates ofx. When the covariance matrix C x is updated using the intermediate estimates ofx, the next iteration of solving forx andασ 2 y is followed, and the final estimatedx can be obtained until the solutions converge.

GRACE Level-1B Data and Background Models
The GRACE Level-1B RL03 data (ftp://isdcftp.gfz-potsdam.de/grace/Level-1B/JPL/) from January 2005 to December 2010 were used to invert regional surface mass variations in this paper. The GRACE Level-1B data used in our study include KBR range-rate measurements (KBR1B) with a sampling rate of 5 s, satellite positions, and velocities (GNV1B) with a sampling rate of 5 s, accelerometer measurements (ACC1B) with a sampling rate of 1 s, and satellite attitude measurements (SCA1B) with a sampling rate of 1 s. To keep the sampling rate consistent with other data, ACC1B and SCA1B data were re-sampled to 5 s before the inversion of surface mass variations. It should be noted that the GRACE Level-1B RL03 data are identical to the data used for solving the GRACE RL06 SH solutions (e.g., GFZ, CSR, and JPL RL06 SH solutions) and the GRACE RL06 mascon solutions (e.g., JPL RL06M and CSR RL06M).
The background models used for estimating intersatellite geopotential differences are summarized in Table 1. These background models are basically consistent with GRACE Level-2 Processing Standards Document (for Level-2 Product Release 06) [6][7][8]. For instance, the static earth gravity filed was modeled by the GGM05C model [46], the third-body perturbations were computed based on JPL planetary and lunar ephemeris DE430 [47], and the solid earth tide, solid earth pole tide, and the contribution of general relativity were calculated according to International Earth Rotation Service (IERS) Conventions (2010) [48]. The EOT11 model [49] was employed to calculate the dynamic effects of ocean tide, and the Ray/Ponte model [50] was used to compute S 1 and S 2 air tidal contributions. The ocean pole tide was evaluated by the self-consistent equilibrium model of Desai [51]. In addition, the short period non-tidal variability in the atmosphere and oceans was reduced using the Atmosphere and Ocean De-aliasing Level-1B (AOD1B) RL06 products (ftp://isdcftp.gfz-potsdam.de/grace/Level-1B/GFZ/AOD/RL06/). The Global Land Data Assimilation System (GLDAS) is a global hydrological model [52], jointly developed by the National Aeronautics and Space Administration's (NASA) GSFC and the National Oceanic and Atmospheric Administration's (NOAA) National Centers of Environmental Prediction (NCEP). The GLDAS model uses land surface modeling and data assimilation technology, and the published data mainly include the inputs and outputs of land surface parameters (soil moisture, soil temperature, evaporation, rainfall, runoff, and snow), and then obtain the near real-time information of land surface variations. In this study, soil moisture (0-200 cm), snow melt water, and plant canopy surface water estimates from the GLDAS Noah model (https://disc.gsfc.nasa.gov/datasets/GLDAS_ NOAH10_M_2.1/summary?keywords=GLDAS) at monthly intervals with a spatial resolution of 1 • × 1 • were used to obtain the changes of terrestrial water storage (TWS). GLDAS-derived TWS changes were used as a priori information to reduce the ill condition for solving regional surface mass variations.
However, the data processing strategies of these three GRACE mascon solutions are different from each other. CSR RL06M are provided on 0.25 • × 0.25 • grids but computed on an equal area geodesic grid comprised of hexagonal tiles approximately 120 km wide (~1 • × 1 • at the equator). Constraints on the CSR mascons include the application of a time-variable regularization matrix from 200-km Gaussian smoothed regularized SH solutions, and separation of land and ocean signals to reduce leakage [26,28]. JPL RL06M are provided on 0.5 • × 0.5 • grids but estimated on a series of equal area 3 • × 3 • spherical caps, and the scale factors are provided to infer spatial variations among the 0.5 • × 0.5 • sub-blocks within each 3 • × 3 • block. The JPL mascons are constrained by a priori information from geophysical models, such as global hydrological model, ocean models, and altimetry data [24], and a coastline resolution improvement (CRI) filter is applied to reduce spatial leakage from land to oceans [25]. The GSFC mascon solutions are estimated on a global set of 1 • arcdeg equal-area blocks, which are the same size of blocks as the final solutions provided. GSFC mascons employed the regularization constraint in the form of the signal auto-covariance matrix, which is constructed by spatio-temporal constraint equations [11,23].
It should be noted that the effective resolution of these mascon solutions is still limited by the native GRACE resolution (i.e.,~300 km), although high-resolution interpolated products are provided, spatial leakage will also be present in mascons [53]. In addition, the above three mascon solutions were applied several corrections, e.g., the C 20 term replacement with Satellite Laser Ranging (SLR) estimates, geocenter correction with the degree-1 coefficients, and the glacial isostatic adjustment (GIA) correction.

GRACE Spherical Harmonic (SH) Solutions
The GRACE SH solutions were also used for comparisons. The CSR RL06 SH solutions (CSR RL06SH) up to the degree and order 96 (ftp://isdcftp.gfz-potsdam.de/grace/Level-2/CSR/RL06/) were used to infer regional surface mass variations, in which Gaussian filtering [3] with the averaging radius of 300 km and the decorrelation filter with P3M6 [14] were employed to mitigate the longitudinal stripes. In order to be consistent with the GRACE mascon solutions for later comparisons, the C 20 term replacement, geocenter correction, and GIA correction were applied to GRACE SH solutions. Specifically, the C 20 coefficients were replaced by the SLR solutions provided by the GRACE Science Data System (SDS) team [54], the degree-1 SH coefficients were added back to the GRACE SH solutions using the degree-1 components provided by the GRACE project as supplementary datasets (GRACE Technical Note 13) [55], and the GIA correction was applied based on the ICE6G-D model from Peltier et al. [56].

Geopotential Differences Estimation Based on the RCR Technique
The geopotential difference T 12 in Equation (7) due to surface mass changes can be obtained after the static geopotential and other high-frequency time-variable gravitational potential (listed in Table 1) are removed based on the energy balance Equation (2). However, the derived residual geopotential differences are contaminated by systematic errors, which are mainly caused by the KBR instrument and orbit errors [41,57,58], as well as the background model errors [39,59].
Because the orbit and background models used in Equation (2) are critical inputs, how to reduce the influences of orbit and background model errors is a key issue for precise estimation of geopotential differences. In this study, we treated the GNV1B orbit data as pseudo observations, and estimated purely dynamic orbits based on the dynamic integral approach. With the background models listed in Table 1, the dynamic orbits of two GRACE satellites were independently integrated every day. Meanwhile, the accelerometer data were simultaneously calibrated with respect to purely dynamic orbits. For each satellite, the estimated parameters include the initial state vectors (three positions and three velocities) and the accelerometer parameters (three biases and three scales in the three directions). When the purely dynamic orbits computed from the known background models were used as the input for energy integral Equation (2), the output geopotential differences would only contain time-variable gravity information contributed from the KBR range-rate measurements. The above similar technique had been used in previous studies on GRACE, such as Shang et al. [39], Zhao et al. [59], and Liu et al. [60].
In order to reduce systematic errors, especially at the low-frequency band (mainly due to KBR instrument errors), empirical parameters (i.e., bias, drift, and periodic components) are widely employed as follows [38,41]: where ∆T 12 are the systematic errors; T rev is the revolution period (~5400 s); p 0 , p 1 , p 2 , p 3 , p 4 , and p 5 are empirical parameters to be estimated, in which p 0 contains integration constant and system bias error, p 1 accounts for the drift error, p 2 and p 3 absorb the 1 CPR systematic errors, while p 4 and p 5 absorb the 2 CPR systematic errors. The parameters p 0 , p 1 , p 2 , and p 3 are estimated for every revolution, while p 4 and p 5 are estimated for every half revolution. Unfortunately, the above empirical parameters used for mitigating systematic errors may absorb part of the temporal gravity signals as well, especially for the long-wavelength components [37,39,41,43]. In order to show the weakening effects of empirical parameters on temporal gravity signals, we did a simulation experiment as follows: First, 1 day of GRACE Level-1B data (e.g., 1 September 2005) were used to compute geopotential differences, and then the systematic errors were estimated through the empirical Equation (14); second, simulated geopotential differences (as the true values) for the day of 1 September 2005 were calculated using CSR RL06 SH solutions (truncated to degree and order 96 with the GGM05C static Earth gravity field removed); third, simulated geopotential differences were added with systematic errors that were estimated in the first step, then the above sum values were used to estimate the systematic errors again by using the empirical Equation (14), and the final estimated geopotential differences were obtained by subtracting the re-estimated systematic errors from the sum values. The two panels of Figure 1 show the (a) time series and (b) power spectral density (PSD) of the simulated geopotential difference using the CSR RL06 SH solutions (black solid curves), estimated geopotential difference using empirical parameters (blue solid curves), and geopotential difference residuals (red solid curves) for the day of 1 September 2005.
As shown in Figure 1a, compared with simulated geopotential differences, the variation ranges of the estimated geopotential difference time series are notably decreased, and the RMS of geopotential differences is reduced from 0.0015 to 0.0013 m 2 /s 2 , which means that about 13% of the geopotential difference signals were absorbed by the empirical parameters. In addition, the residuals time series is characterized by mostly low-frequency variations (due to the empirical parameters), and the influence of geopotential difference residuals cannot be ignored. From Figure 1b, we can see that the PSD of the estimated geopotential differences is decreased in the low-frequency band (especially at 1 CPR and 2 CPR), suggesting that the weakening effects of empirical parameters are mainly on long-wavelength temporal signals. It should be noted that, compared with the PSD of estimated geopotential differences, obvious decreases are observed in the simulated geopotential differences at the high-frequency band (up to 96 CPR). This is because the truncated CSR RL06 SH solutions up to the degree and order 96 were used for computing the simulated geopotential differences, which do not contain high-frequency temporal gravity signals. Remote Sens. 2020, 11, x FOR PEER REVIEW 10 of 29 In order to extract precise geopotential differences, and with the aim to preserve both the lowand high-frequency temporal gravity signals, the RCR technique was employed during the estimation of geopotential differences. In the present study, geopotential differences were computed in the following three steps. Firstly, the reference geopotential difference 12 was calculated by using CSR RL06 SH solutions up to the degree and order 50 (corresponding to a 400-km spatial resolution), in which the C20 coefficients were replaced by SLR solutions, and then 12 was removed from the geopotential difference 12 that was computed based on the energy balance Equation (2). Secondly, through the empirical parameters fitting (EPF) using Equation (14), the systematic errors were estimated by using 12 = EPF[ 12 − 12 ], and the corrected geopotential difference was obtained by ( 12 − 12 ) − 12 . Thirdly, the reference geopotential difference 12 was added back to the corrected geopotential difference ( 12 − 12 ) − 12 , and the final estimated geopotential difference was 12 − 12 .
Compared with the RCR technique, the traditional method for estimating the systematic errors is 12 = EPF [ 12 ], and the estimated geopotential difference is 12 − 12 . In order to illustrate the improvement of the RCR technique for estimation of the geopotential differences, Figures 2(a) and 2b show the global map of geopotential differences estimated by the traditional method and RCR technique for the month of September 2005 with the GGM05C static earth gravity field removed. Figure 2c shows the geopotential difference residuals obtained by subtracting estimates of the traditional method from those of the RCR technique. It can be seen from Figures 2a,b that geopotential differences are the direct reflections of surface mass changes, and signals estimated based on the RCR technique are obviously stronger than those of the traditional method in South America, Antarctic Peninsula, southern Africa, and Australia. However, there are also some regions with weakened signals, such as eastern Asia and North Atlantic. The RMS of geopotential differences (outliers are removed according to a 3-sigma rule) estimated by the traditional method and RCR technique are In order to extract precise geopotential differences, and with the aim to preserve both the low-and high-frequency temporal gravity signals, the RCR technique was employed during the estimation of geopotential differences. In the present study, geopotential differences were computed in the following three steps. Firstly, the reference geopotential difference T re f 12 was calculated by using CSR RL06 SH solutions up to the degree and order 50 (corresponding to a 400-km spatial resolution), in which the C 20 coefficients were replaced by SLR solutions, and then T re f 12 was removed from the geopotential difference T 12 that was computed based on the energy balance Equation (2). Secondly, through the empirical parameters fitting (EPF) using Equation (14), the systematic errors were estimated by using ∆T RCR 12 = EPF T 12 − T re f 12 , and the corrected geopotential difference was obtained by ( Thirdly, the reference geopotential difference T re f 12 was added back to the corrected geopotential difference (T 12 − T re f 12 ) − ∆T RCR 12 , and the final estimated geopotential difference was T 12 − ∆T RCR 12 . Compared with the RCR technique, the traditional method for estimating the systematic errors is ∆T 12 = EPF[T 12 ], and the estimated geopotential difference is T 12 − ∆T 12 . In order to illustrate the improvement of the RCR technique for estimation of the geopotential differences, Figure 2a,b show the global map of geopotential differences estimated by the traditional method and RCR technique for the month of September 2005 with the GGM05C static earth gravity field removed. Figure 2c shows the geopotential difference residuals obtained by subtracting estimates of the traditional method from those of the RCR technique. It can be seen from Figure 2a,b that geopotential differences are the direct reflections of surface mass changes, and signals estimated based on the RCR technique are obviously stronger than those of the traditional method in South America, Antarctic Peninsula, southern Africa, and Australia. However, there are also some regions with weakened signals, such as eastern Asia and North Atlantic. The RMS of geopotential differences (outliers are removed according to a 3-sigma rule) estimated by the traditional method and RCR technique are~0.0022 and 0.0025 m 2 /s 2 , respectively. Figure 2c shows that geopotential difference residuals are mainly characterized by low-frequency gravity signals, which are absorbed by the empirical parameters used in the traditional method. Therefore, the RCR technique is helpful to retain the low-frequency gravity signals, and the inversion of surface mass variations at a full spectrum (both low and high frequency) can be obtained.
Remote Sens. 2020, 11, x FOR PEER REVIEW 11 of 29 ~0.0022 and 0.0025 m 2 /s 2 , respectively. Figure 2c shows that geopotential difference residuals are mainly characterized by low-frequency gravity signals, which are absorbed by the empirical parameters used in the traditional method. Therefore, the RCR technique is helpful to retain the lowfrequency gravity signals, and the inversion of surface mass variations at a full spectrum (both low and high frequency) can be obtained.

Preliminary Tests over South America
In order to verify the proposed method (and software), regional surface mass variations on 2° × 2° geographic grids over South America ([60°S-15°N, 90°W -30°W], see Figure3a) were inverted from geopotential difference estimates in September 2005, and GLDAS-derived TWS changes were employed as a priori constraints to reduce the ill condition in regional mass estimations from GRACE. In this study, the optimal regularization parameters were determined based on the iterative leastsquares estimation method. Specifically, we set initial values = 1 and 2 = 0.002 m 2 /s 2 in Equation (10), and took the maximum value of EWH differences between two adjacent iterations to be less than 10 −6 mm as the convergence condition. The final estimated mass variations can be obtained after about 10 iterations. Figures 3b,c show the regional surface mass variations over South America inverted

Preliminary Tests over South America
In order to verify the proposed method (and software), regional surface mass variations on 2 Figure 3a) were inverted from geopotential difference estimates in September 2005, and GLDAS-derived TWS changes were employed as a priori constraints to reduce the ill condition in regional mass estimations from GRACE. In this study, the optimal regularization parameters were determined based on the iterative least-squares estimation method. Specifically, we set initial values α = 1 and σ 2 y = 0.002 m 2 /s 2 in Equation (10), and took the maximum value of EWH differences between two adjacent iterations to be less than 10 −6 mm as the convergence condition. The final estimated mass variations can be obtained after about 10 iterations. Figure 3b,c show the regional surface mass variations over South America inverted from geopotential differences estimated by the traditional method and RCR technique in September 2005, respectively. Figure 3d presents the corresponding surface mass variations over South America derived from the CSR RL06 SH solutions with 300-km Gaussian smoothing and P3M6 decorrelation filtering for comparisons.
Remote Sens. 2020, 11, x FOR PEER REVIEW 12 of 29 from geopotential differences estimated by the traditional method and RCR technique in September 2005, respectively. Figure 3d presents the corresponding surface mass variations over South America derived from the CSR RL06 SH solutions with 300-km Gaussian smoothing and P3M6 decorrelation filtering for comparisons. As shown in Figure 3, there are obvious differences in regional surface mass variations inverted from geopotential differences estimated by the traditional method and RCR technique. The mass variations over Amazon basin show substantially stronger negative anomalies in Figure 3c than in Figure 3b, while the mass variations over Orinoco basin (in the north), western part of Parana basin, and parts of southern South America (e.g., in Uruguay) show stronger positive anomalies in Figure  3b than in Figure 3c. However, compared with surface mass variations inverted from geopotential differences estimated by the traditional method, the results from geopotential differences estimated by the RCR technique are more consistent with CSR RL06 SH solutions. Moreover, mass variation signals in Figures 3c,d are stronger than those in Figure 3b, and the RMS of mass variations in Figures 3b-d are 13.7, 18.4, and 16.7 cm, respectively. The reason, as mentioned before, is that the empirical parameters used in the traditional method absorbed the long-wavelength gravity signals, and thus weakened the mass variation estimates. In addition, the signals of mass variations in Figure 3c are stronger than those in Figure 3d, which strongly suggests that the regional mascon solutions can retain more signal energy and enhance regional mass variation inversion. As shown in Figure 3, there are obvious differences in regional surface mass variations inverted from geopotential differences estimated by the traditional method and RCR technique. The mass variations over Amazon basin show substantially stronger negative anomalies in Figure 3c than in Figure 3b, while the mass variations over Orinoco basin (in the north), western part of Parana basin, and parts of southern South America (e.g., in Uruguay) show stronger positive anomalies in Figure 3b than in Figure 3c. However, compared with surface mass variations inverted from geopotential differences estimated by the traditional method, the results from geopotential differences estimated by the RCR technique are more consistent with CSR RL06 SH solutions. Moreover, mass variation signals in Figure 3c,d are stronger than those in Figure 3b, and the RMS of mass variations in Figure 3b-d are 13.7, 18.4, and 16.7 cm, respectively. The reason, as mentioned before, is that the empirical parameters used in the traditional method absorbed the long-wavelength gravity signals, and thus weakened the mass variation estimates. In addition, the signals of mass variations in Figure 3c are stronger than those in Figure 3d, which strongly suggests that the regional mascon solutions can retain more signal energy and enhance regional mass variation inversion.
To assess the effects of different block sizes on the inversion results, Figure 4 shows regional mass variations estimated on 1 • × 1 • grids, 2 • × 2 • grids, and 4 • × 4 • grids from geopotential differences computed by the RCR technique in September 2005 using GLDAS-derived TWS changes as a priori constraints, and the corresponding residuals between geopotential difference-estimated mass variations and GLDAS-derived TWS changes are also presented. Figure 5 shows the residuals between regional mass variations estimated on different mascon block sizes.  To assess the effects of different block sizes on the inversion results, Figure 4 shows regional mass variations estimated on 1° × 1° grids, 2° × 2° grids, and 4° × 4° grids from geopotential differences computed by the RCR technique in September 2005 using GLDAS-derived TWS changes as a priori constraints, and the corresponding residuals between geopotential difference − estimated mass variations and GLDAS-derived TWS changes are also presented. Figure 5 shows the residuals between regional mass variations estimated on different mascon block sizes.
As shown in Figure 4, although the a priori information from GLDAS TWS changes can contribute to the final estimated mass variations, there are obvious differences between regional mass variation solutions and GLDAS − derived TWS changes, no matter what sizes of grids were used for the estimation. The residuals can be considered as the new contributions from GRACE observations to the final estimated regional mass variations, and the added value from GRACE data can obviously enhance the signals of mass variations over South America, especially in basins of Amazon, Orinoco, between regional mass variations estimated on different mascon block sizes. The RMS differences of the inversion results between 1° × 1° and 2° × 2° grids, 1° × 1° and 4° × 4° grids, 2 × 2° and 4° × 4° grids are ~1.6, 3.3, and 3.1 cm, respectively. Considering the intrinsic optimal spatial resolution of the GRACE data is limited at ~200-300 km, and the 1° × 1° regional mascon solution will have leakage errors, while the 2° × 2° grids computation seem enough to represent the GRACE resolution.

Spatial Distribution of Mass Changes over South America
To examine the performance of regional mass variations inverted from GRACE intersatellite geopotential differences, we employed the GRACE Level − 1B RL03 data from January 2005 to December 2010 to calculate the precise geopotential differences based on the RCR technique, and estimated regional surface mass variations on 2° × 2°grids at monthly intervals over South America. We then compared the results with the official GRACE RL06 mascon solutions and SH solutions. It should be noted that all the mascon solutions and SH solutions had been applied with the corrections as mentioned in Section 2.4.4 (i.e., the C20 term replacement, the geocenter correction, and the GIA correction, unless otherwise mentioned). Figure 6 shows the spatial distributions of regional mass variations over South America in January (the first row), April (the second row), July (the third row), and October (the fourth row) 2005 from "the new" geopotential difference − based regional mascon solutions (hereinafter called GPD Mascon, the first column), CSR RL06 mascon solutions (CSR RL06M, the second column), JPL RL06 mascon solutions (JPL RL06M, the third column), GSFC mascon solutions (GSFC Mascon, the fourth column), and CSR RL06 SH solutions (CSR RL06SH, the fifth column), respectively.
As shown in Figure 6, both of the mascon solutions and the SH solutions can reflect the same spatial distribution characteristics of surface mass variations over South America in the four selected months. However, due to the spatial filtering in post − processing of GRACE SH solutions, the amplitudes of surface mass variations from CSR RL06SH are obviously smaller than those of the four mascon solutions. Moreover, there are apparent residual longitudinal stripes in CSR RL06SH (especially in January and April 2005), although the post − processing was applied to these SH solutions. The reason is that, unlike the mascon solutions, the GRACE SH solutions are solved without using any constraints or a priori information. As shown in Figure 4, although the a priori information from GLDAS TWS changes can contribute to the final estimated mass variations, there are obvious differences between regional mass variation solutions and GLDAS-derived TWS changes, no matter what sizes of grids were used for the estimation. The residuals can be considered as the new contributions from GRACE observations to the final estimated regional mass variations, and the added value from GRACE data can obviously enhance the signals of mass variations over South America, especially in basins of Amazon, Orinoco, and Colorado, etc. Small mascon block sizes are expected to help to improve the level of details in estimated mass variations (e.g., from 4 • × 4 • to 2 • × 2 • grids), but there are no notable changes of signal amplitudes by decreasing the size of the mascon blocks, and no gains of details in our experiments (see Figure 4a-c). In addition, we can draw a similar conclusion from Figure 5 based on the residuals between regional mass variations estimated on different mascon block sizes. The RMS differences of the inversion results between 1 • × 1 • and 2 • × 2 • grids, 1 • × 1 • and 4 • × 4 • grids, 2 × 2 • and 4 • × 4 • grids are~1.6, 3.3, and 3.1 cm, respectively. Considering the intrinsic optimal spatial resolution of the GRACE data is limited at~200-300 km, and the 1 • × 1 • regional mascon solution will have leakage errors, while the 2 • × 2 • grids computation seem enough to represent the GRACE resolution.

Spatial Distribution of Mass Changes over South America
To examine the performance of regional mass variations inverted from GRACE intersatellite geopotential differences, we employed the GRACE Level-1B RL03 data from January 2005 to December 2010 to calculate the precise geopotential differences based on the RCR technique, and estimated regional surface mass variations on 2 • × 2 • grids at monthly intervals over South America. We then compared the results with the official GRACE RL06 mascon solutions and SH solutions. It should be noted that all the mascon solutions and SH solutions had been applied with the corrections as mentioned in Section 2.4.4 (i.e., the C 20 term replacement, the geocenter correction, and the GIA correction, unless otherwise mentioned). Figure 6 shows the spatial distributions of regional mass variations over South America in January (the first row), April (the second row), July (the third row), and October (the fourth row) 2005 from "the new" geopotential difference-based regional mascon solutions (hereinafter called GPD Mascon, the first column), CSR RL06 mascon solutions (CSR RL06M, the second column), JPL RL06 mascon solutions (JPL RL06M, the third column), GSFC mascon solutions (GSFC Mascon, the fourth column), and CSR RL06 SH solutions (CSR RL06SH, the fifth column), respectively.  Table 2, and the ranges of RMS differences between the three official mascon solutions are also presented for comparisons. These results show that GPD Mascon appears closer to GSCF Mascon but differs much more from JPL RL06M. For example, the RMS differences between the GPD Mascon and CSR RL06M, Figure 6. Spatial distributions of regional mass variations over South America in January (the first row), April (the second row), July (the third row), and October (the fourth row) 2005 from geopotential difference-based regional mascon solutions (GPD Mascon, the first column), CSR RL06 mascon solutions (CSR RL06M, the second column), JPL RL06 mascon solutions (JPL RL06M, the third column), GSFC mascon solutions (GSFC Mascon, the fourth column), and CSR RL06 SH solutions (CSR RL06SH, the fifth column), respectively.
As shown in Figure 6, both of the mascon solutions and the SH solutions can reflect the same spatial distribution characteristics of surface mass variations over South America in the four selected months. However, due to the spatial filtering in post-processing of GRACE SH solutions, the amplitudes of surface mass variations from CSR RL06SH are obviously smaller than those of the four mascon solutions. Moreover, there are apparent residual longitudinal stripes in CSR RL06SH (especially in January and April 2005), although the post-processing was applied to these SH solutions. The reason is that, unlike the mascon solutions, the GRACE SH solutions are solved without using any constraints or a priori information.
Comparing the four mascon solutions in Figure 6 Table 2, and the ranges of RMS differences between the three official mascon solutions are also presented for comparisons. These results show that GPD Mascon appears closer to GSCF Mascon but differs much more from JPL RL06M. For example, the RMS differences between the GPD Mascon and CSR RL06M, GPD Mascon and JPL RL06M, and GPD Mascon and GSFC Mascon in April and October 2005 are 6.18 vs. 6.72 cm, 7.25 vs. 7.93 cm, and 5.42 vs. 5.46 cm, respectively. Meanwhile, differences among CSR RL06M, JPL RL06M, and GSCF Mascon are also obvious, e.g., the ranges of RMS differences between the three official mascon solutions are 5.87-6.75 cm in April 2005, and 5.67-6.70 cm in October 2005, respectively. Table 2. Statistics of the root-mean-square (RMS) differences between the geopotential difference-  In order to better understand the surface mass variations inferred from different solutions over South America, a simultaneous fit of annual, semi-annual, trend, and bias was made to the time series for each mascon block of the above five solutions (i.e., GPD Mascon, CSR RL06M, JPL RL06M, GSFC Mascon, and CSR RL06SH). Figure 7 shows the spatial distributions of annual amplitudes (the first row), trends (the second row), root-mean-square error (RMSE, the third row), and signal-to-noise ratio (SNR, the fourth row) of mass variations over South America (over the period January 2005 to December 2010) derived from the five solutions, respectively.

RMS Differences (cm) Ranges of RMS Differences (cm) GPD vs. CSR GPD vs. JPL GPD vs. GSFC
As shown in Figure 7, the spatial patterns of the annual amplitudes and trends for the five solutions are consistent with each other, but the differences between each other can still be observed if we carefully examine these results at basin scale. An apparent observation is that annual amplitudes and trends over the Amazon basin are significantly larger in the four mascon solutions as compared to those from the CSR RL06SH results. Because the four mascon solutions are provided on different sizes of grids (not the real spatial resolutions), for better comparisons, we resampled the GPD Mascon, JPL RL06M, and GSFC Mascon into 0.25 • × 0.25 • grids, consistent with the CSR RL06M grid's resolution. In terms of annual amplitudes, the correlation coefficients between GPD Mascon and the three other mascons (CSR RL06M, JPL RL06M, and GSFC Mascon) are~0.95, 0.93, and 0.97, respectively. In terms of trends, the corresponding correlation coefficients are~0.85, 0.79, and 0.88, respectively. The relatively lower correlation coefficients for the trends are attributed to the fact that temporal mass variations in South America are dominated by seasonal hydrological signals, and the long-term trend is usually small, and its uncertainty is relatively larger. Nevertheless, our results demonstrate that regional mass variations inverted from geopotential differences estimated by the RCR technique are comparable to those from the official GRACE mascon solutions, without using low degrees of SH solutions to complement the missing long-wavelength signal (unlike in the studies, e.g., Frédéric Frappart et al. [43] and Ramillien et al. [44]). In addition, the RMSE of mass variation residuals and SNR were used to evaluate the errors and accuracies of different solutions. As mentioned above, a simultaneous fit was made to the time series for each mascon block, and mass variation residuals were obtained by removing the fitted signals from the mass change time series. Then, the residual RMSE for each block was estimated, and the SNR was computed as the quotient between the annual amplitude and RMSE. As can be seen from Figure 7, the RMSE patterns from the five solutions are similar in most of the regions. For example, the mean RMSEs in the Amazon basin are ~4. 23, 4.24, 4.46, 4.11, and 3.34 cm for GPD Mascon, CSR RL06M, JPL RL06M, GSFC Mascon, and CSR RL06SH, respectively. Here, due to the spatial smoothing and filtering applied to the SH solutions, the mean RMSEs from CSR RL06SH in the Amazon basin are the smallest (compared with those of the mascon solutions). The SNR patterns from the five solutions are similar to those of the annual amplitudes shown in the top row of Figure  7. The mean SNR in the Amazon basin are ~4.65, 4.55, 4.56, 4.90, and 4.51 for GPD Mascon, CSR RL06M, JPL RL06M, GSFC Mascon, and CSR RL06SH, respectively. Again, these comparisons also suggest that the GPD Mascon solutions are comparable to the official GRACE mascons, and further demonstrate the effectiveness of solving regional mascon solutions using geopotential differences estimated by the RCR technique with a priori constraints. In addition, the RMSE of mass variation residuals and SNR were used to evaluate the errors and accuracies of different solutions. As mentioned above, a simultaneous fit was made to the time series for each mascon block, and mass variation residuals were obtained by removing the fitted signals from the mass change time series. Then, the residual RMSE for each block was estimated, and the SNR was computed as the quotient between the annual amplitude and RMSE. As can be seen from Figure 7, the RMSE patterns from the five solutions are similar in most of the regions. For example, the mean RMSEs in the Amazon basin are~4. 23, 4.24, 4.46, 4.11, and 3.34 cm for GPD Mascon, CSR RL06M, JPL RL06M, GSFC Mascon, and CSR RL06SH, respectively. Here, due to the spatial smoothing and filtering applied to the SH solutions, the mean RMSEs from CSR RL06SH in the Amazon basin are the smallest (compared with those of the mascon solutions). The SNR patterns from the five solutions are similar to those of the annual amplitudes shown in the top row of Figure 7. The mean SNR in the Amazon basin are~4.65, 4.55, 4.56, 4.90, and 4.51 for GPD Mascon, CSR RL06M, JPL RL06M, GSFC Mascon, and CSR RL06SH, respectively. Again, these comparisons also suggest that the GPD Mascon solutions are comparable to the official GRACE mascons, and further demonstrate the effectiveness of solving regional mascon solutions using geopotential differences estimated by the RCR technique with a priori constraints.

Basin-Scale Mass Changes over South America
In order to compare surface mass variations derived from different solutions over South America at the basin scale, mass variations over the eight largest basins (i.e., Amazon, Parana, Orinoco, Tocantins, San Francisco, Colorado, Rio Pamaiba, and Salado, see Figure 3a) were further analyzed. Figure 8 shows mass change time series for the eight basins over the period January 2005 to December 2010, derived from GPD Mascon, CSR RL06M, JPL RL06M, GSFC Mascon, and CSR RL06SH, respectively. Figure 9 presents the non-seasonal mass change time series for the eight basins, which were obtained by removing seasonal signals and long-term trends from the time series shown in Figure 8. Table 3 lists the corresponding annual amplitudes and trends of surface mass change time series for these eight basins, and the uncertainties correspond to a 95% confidence interval.

Basin-Scale Mass Changes over South America
In order to compare surface mass variations derived from different solutions over South America at the basin scale, mass variations over the eight largest basins (i.e., Amazon, Parana, Orinoco, Tocantins, San Francisco, Colorado, Rio Pamaiba, and Salado, see Figure 3a) were further analyzed. Figure 8 shows mass change time series for the eight basins over the period January 2005 to December 2010, derived from GPD Mascon, CSR RL06M, JPL RL06M, GSFC Mascon, and CSR RL06SH, respectively. Figure 9 presents the non − seasonal mass change time series for the eight basins, which were obtained by removing seasonal signals and long − term trends from the time series shown in Figure 8. Table 3 lists the corresponding annual amplitudes and trends of surface mass change time series for these eight basins, and the uncertainties correspond to a 95% confidence interval. As can be seen from Figure 8, basin mass change time series from the five solutions (i.e., GPD Mascon, CSR RL06M, JPL RL06M, GSFC Mascon, and CSR RL06SH) agree well with each other, especially for the large river basins with stronger seasonal signals (e.g., in Amazon, Parana, and Orinoco). However, discrepancies are also obvious, especially in some small basins with weaker mass change signals (e.g., in Colorado and Salado). This is because the uncertainty estimated in small basins with weaker mass change signals is relatively large compared to mass change signals themselves. In addition, due to the spatial filtering applied to the SH solutions, mass change signals from CSR RL06SH over the eight basins are slightly smaller than those of the mascon solutions (which can also be found in Table 3). three mascons (CSR RL06M, JPL RL06M, and GSFC Mascon) still reach up to 0.99, 0.98, and 0.97, respectively. In the Salado basin (the smallest of the eight basins), after removing the seasonal signals and linear trends, the correlation coefficients of mass change time series between GPD Mascon and CSR RL06M, JPL RL06M, and GSFC Mascon are ~0.91, 0.91, and 0.86, respectively. The excellent consistency among basin − scale mass changes from the four mascon solutions indicate that all these four mascon solutions perform fairly well at both seasonal and non − seasonal time scales. The linear trends of the mass changes over the eight basins are small, and there are obvious differences between the estimated trends from these five different solutions (see Table 3). Moreover, the uncertainties for the estimated trends are relatively larger, and in some cases the estimated uncertainty is even greater than the trend itself. The obvious differences in trends can also be confirmed from the lower spatial correlation of trends as analyzed in Section 4.2 (see Figure 7). As mentioned before, basin − scale mass variations over South America are dominated by seasonal TWS  Table 3. Statistics of the annual amplitudes and trends of surface mass change time series for the eight largest basins during the period from January 2005 to December 2010 over South America, and the uncertainties correspond to a 95% confidence interval. The mass change time series were derived from geopotential difference-based regional mascon solutions (GPD Mascon), CSR RL06 mascon solutions (CSR RL06M), JPL RL06 mascon solutions (JPL RL06M), GSFC mascon solutions (GSFC Mascon), and CSR RL06 SH solutions (CSR RL06SH), respectively. As can be seen from Figure 8, basin mass change time series from the five solutions (i.e., GPD Mascon, CSR RL06M, JPL RL06M, GSFC Mascon, and CSR RL06SH) agree well with each other, especially for the large river basins with stronger seasonal signals (e.g., in Amazon, Parana, and Orinoco). However, discrepancies are also obvious, especially in some small basins with weaker mass change signals (e.g., in Colorado and Salado). This is because the uncertainty estimated in small basins with weaker mass change signals is relatively large compared to mass change signals themselves. In addition, due to the spatial filtering applied to the SH solutions, mass change signals from CSR RL06SH over the eight basins are slightly smaller than those of the mascon solutions (which can also be found in Table 3).
As listed in Table 3, the annual amplitudes for the four mascon solutions agree well, especially in basins with stronger mass change signals (e.g., the Amazon, Parana, Orinoco, Tocantins, San Francisco, and Rio Pamaiba), and the differences of annual amplitudes are acceptable as compared to their uncertainties. Aside from the seasonal variability, basin-scale non-seasonal signals affected by extreme climate events, such as floods and droughts, are also of great interest. Figure 9 shows that the non-seasonal signals of the four mascon solutions for the eight basins are also consistent. The linear trends of the mass changes over the eight basins are small, and there are obvious differences between the estimated trends from these five different solutions (see Table 3). Moreover, the uncertainties for the estimated trends are relatively larger, and in some cases the estimated uncertainty is even greater than the trend itself. The obvious differences in trends can also be confirmed from the lower spatial correlation of trends as analyzed in Section 4.2 (see Figure 7). As mentioned before, basin-scale mass variations over South America are dominated by seasonal TWS changes, and trends are usually small and uncertainties are relatively larger, leading to poorer agreements in trends.
In order to further analyze the consistency among the four mascon solutions, a total of 30 river basins (with areas larger than 50,000 km 2 ) over South America were selected for comparisons. Mass change time series of the basin average for the 30 basins over the period January 2005 to December 2010 were calculated, and simultaneous fits of seasonal terms and trends (mentioned in Section 4.2) were made to these time series. Figure 10  As shown in Figure 10, if the two mascon solutions in one panel agree perfectly, the filled cycle would align on the y = x line (black solid lines). However, if the scattered points for any two mascon solutions are distributed around that line, the corresponding fitting curve is the linear regression equation y = kx + b (given by the red dashed line), and the relative positions of the black solid line and the red dashed line can describe the consistency of annual amplitudes between the corresponding two solutions. We can see that for every panel in Figure 10, the red dashed lines are very close to the black solid lines and the correlation coefficients are all greater than 0.96, which indicate that the consistencies of annual amplitudes for the four mascon solutions are very good, even though notable differences can also be observed. equation y = kx + b (given by the red dashed line), and the relative positions of the black solid line and the red dashed line can describe the consistency of annual amplitudes between the corresponding two solutions. We can see that for every panel in Figure 10, the red dashed lines are very close to the black solid lines and the correlation coefficients are all greater than 0.96, which indicate that the consistencies of annual amplitudes for the four mascon solutions are very good, even though notable differences can also be observed.  In general, when the slope of the red dashed line is less than 1.0, the y-axis values are smaller than those of the x-axis, e.g., GPD Mascon vs. JPL RL06M, GPD Mascon vs. GSFC Mascon, and CSR RL06M vs. JPL RL06M, as shown in Figure 10b-d, and vice versa from those in Figure 10a,e,f. In a particular case of GPD Mascon vs. CSR RL06M (Figure 10a), although the slope is greater than 1.0, the red dashed line is generally below the solid black line. In the sense of statistics, it means that the annual amplitudes of the 30 basins from GPD Mascon are slightly smaller than those from CSR RL06M. The relatively larger discrepancies in small basins among the four mascon solutions appear to be the main reason affecting the consistency.

Discussion
This study focused on applying the improved method to estimate regional mascon solutions using GRACE intersatellite geopotential differences, and the inversion results were examined through comprehensive comparisons with the official GRACE mascon solutions. A proper validation was not carried out through the use of independent in situ data. However, the comparative analysis shows that our geopotential difference-based regional mascon solutions (GPD Mascon) are comparable to the three official GRACE mascon solutions at the basin scale (i.e., the accuracies of these four mascon solutions are very similar), and have demonstrated the feasibility and effectiveness (and advantage) of the proposed method.
As can be seen from the comparative analysis, the spatial patterns and characteristics of mass variations over South America from GPD Mascon agree well with those from the three official GRACE mascon solutions (see Figures 6 and 7), and the annual amplitudes of mass change time series from the four mascon solutions at the basin scale over South America show good consistency with each other (see Table 3 and Figure 10), and the non-seasonal signals from these four mascon solutions in larger basins are also in good agreement with each other (see Figure 9). However, notable differences among the four mascon solutions still exist, in particular at the basin scale. The reasons for this are twofold. Firstly, the mass variations inversion methods are different, i.e., all the official GRACE mascon solutions are directly solved from instrument observations based on the dynamic method, while the GPD Mascon solutions are computed from intermediate observations (i.e., intersatellite geopotential differences), which are estimated based on the energy balance approach. Therefore, different approximate errors in inversion models will lead to differences in the final estimated mass variations. Secondly, different processing strategies (e.g., regularization constraints, background models, sizes of mascon blocks, and data pre-processing, etc.) can also induce discrepancies between these mascon inversion results. For instance, even though the three official GRACE mascon solutions are all solved by the dynamic method, differences among them are also obvious. The main reason accounting for this includes the different regularization constraints and different sizes of mascon blocks for the three mascon solutions as introduced in Section 2.4.3.
In our case study, surface mass changes over South America are dominated by seasonal hydrological signals, thus TWS changes from hydrological models (i.e., GLDAS) are chosen as a priori information for the estimation of mascon solutions. As mentioned in Section 4.1, the a priori information from GLDAS-derived TWS changes can contribute to the final estimated mascon solutions, but how the a priori constraints can affect the output needs to be assessed. It is well known that hydrological models suffer from large uncertainties in regions that have poor-quality data (e.g., in high-mountain Asia), or even no data (e.g., in Greenland and Antarctic) [61][62][63]. Therefore, it is necessary to analyze the influences of different a priori information on the estimation of mascon solutions, especially if only the rough a priori information can be obtained. To assess the effects of different a priori constraints on the inversion results, we employed CSR RL06 SH solutions post-processed with different Gaussian smoothing radius as a priori models to estimate the mascons, and compared them with the previous results, which use GLDAS-derived TWS changes as a priori models. Figure 11 shows the regional mass variations estimated from geopotential differences in September 2005 by using different a priori models, and the corresponding residuals between estimated mass variations and the a priori models are presented. Figure 12 shows the corresponding residuals between the geopotential difference-estimated mascon solutions by using different a priori models, and also the residuals between the a priori models. Remote Sens. 2020, 11, x FOR PEER REVIEW 23 of 29 In our case study, surface mass changes over South America are dominated by seasonal hydrological signals, thus TWS changes from hydrological models (i.e., GLDAS) are chosen as a priori information for the estimation of mascon solutions. As mentioned in Section 4.1, the a priori information from GLDAS − derived TWS changes can contribute to the final estimated mascon solutions, but how the a priori constraints can affect the output needs to be assessed. It is well known that hydrological models suffer from large uncertainties in regions that have poor − quality data (e.g., in high − mountain Asia), or even no data (e.g., in Greenland and Antarctic) [61][62][63]. Therefore, it is necessary to analyze the influences of different a priori information on the estimation of mascon solutions, especially if only the rough a priori information can be obtained. To assess the effects of different a priori constraints on the inversion results, we employed CSR RL06 SH solutions post − processed with different Gaussian smoothing radius as a priori models to estimate the mascons, and compared them with the previous results, which use GLDAS − derived TWS changes as a priori models. Figure 11 shows the regional mass variations estimated from geopotential differences in September 2005 by using different a priori models, and the corresponding residuals between estimated mass variations and the a priori models are presented. Figure 12 shows the corresponding residuals between the geopotential difference − estimated mascon solutions by using different a priori models, and also the residuals between the a priori models. (d-f) show the a priori models corresponding to (a-c); (g-i) show the corresponding residuals between estimated mass variations and the a priori models.
As shown in Figure 11a-c, there are only minor differences between the mascon solutions estimated by our improved method through using different a priori models. When taking the GLDAS-derived TWS changes, the CSR RL06 SH solution with 300-km Gaussian filtering, and the CSR RL06 SH solution with 750-km Gaussian filtering as a priori models, the RMS residuals between the estimated mass variations and the a priori models are~8.09, 4.43, and 8.45 cm, respectively. This is true that the a priori model of CSR RL06 SH solutions with 300-km Gaussian filtering is closer to the estimated mascons. However, as shown in Figure 12, although there are obvious differences between the three a priori models (e.g., RMS residuals is larger than 7.0 cm), the estimated mascon solutions agree well with each other (e.g., RMS residuals is~2.0 cm). This shows that our improved method can robustly estimate the mascon solutions, which are less affected by a priori information (even by using a rough a priori information). In other words, GRACE SH solution-derived TWS changes can be used as alternative a priori models when the quality of hydrological models is poor, or no a priori information can be obtained. Meanwhile, the resolution matrix [64] (i.e., (A T A + ασ 2 y C −1 x ) −1 ασ 2 y C −1 x derived from Equation (10)) was used to analyze the relative contribution of a priori constraints to the mascon solutions. Our results show that when using GLDAS-derived TWS changes, CSR RL06 SH solutions with 300-km Gaussian filtering, and CSR RL06 SH solutions with 750-km Gaussian filtering as a priori models, the average contributions of a priori constraints to the final mascon solutions are~11.8%, 12.9%, and 10.5%, respectively. It indicates that the contributions of a priori constraints to the mascon solutions are not dominant, and most contributions are from GRACE observations. Remote Sens. 2020, 11, x FOR PEER REVIEW 24 of 29 Figure 12. Residuals between the geopotential difference-estimated mascon solutions by: (a) using GLDAS-derived TWS and CSR RL06 SH solutions with 300-km Gaussian filtering as a priori models, (b) using GLDAS − derived TWS and CSR RL06 SH solutions with 750-km Gaussian filtering as a priori models; Residuals between the a priori models: (c) GLDAS − derived TWS and CSR RL06 SH solutions with 300-km Gaussian filtering, (d) GLDAS − derived TWS and CSR RL06 SH solutions with 750-km Gaussian filtering.
As shown in Figures 11a-c, there are only minor differences between the mascon solutions estimated by our improved method through using different a priori models. When taking the GLDAS-derived TWS changes, the CSR RL06 SH solution with 300-km Gaussian filtering, and the CSR RL06 SH solution with 750-km Gaussian filtering as a priori models, the RMS residuals between the estimated mass variations and the a priori models are ~8.09, 4.43, and 8.45 cm, respectively. This is true that the a priori model of CSR RL06 SH solutions with 300-km Gaussian filtering is closer to the estimated mascons. However, as shown in Figure 12, although there are obvious differences between the three a priori models (e.g., RMS residuals is larger than 7.0 cm), the estimated mascon solutions agree well with each other (e.g., RMS residuals is ~2.0 cm). This shows that our improved method can robustly estimate the mascon solutions, which are less affected by a priori information (even by using a rough a priori information). In other words, GRACE SH solution-derived TWS changes can be used as alternative a priori models when the quality of hydrological models is poor, or no a priori information can be obtained. Meanwhile, the resolution matrix [64] (i.e., ( + 2 −1 ) −1 2 −1 derived from Equation (10)) was used to analyze the relative contribution of a priori constraints to the mascon solutions. Our results show that when using GLDAS-derived TWS changes, CSR RL06 SH solutions with 300-km Gaussian filtering, and CSR RL06 SH solutions with 750 − km Gaussian filtering as a priori models, the average contributions of a priori constraints to the final mascon solutions are ~11.8%, 12.9%, and 10.5%, respectively. It indicates that the contributions of a priori constraints to the mascon solutions are not dominant, and most contributions are from GRACE observations.

Conclusions
In this study, we presented an improved method for the estimation of regional mass variations Figure 12. Residuals between the geopotential difference-estimated mascon solutions by: (a) using GLDAS-derived TWS and CSR RL06 SH solutions with 300-km Gaussian filtering as a priori models, (b) using GLDAS-derived TWS and CSR RL06 SH solutions with 750-km Gaussian filtering as a priori models; Residuals between the a priori models: (c) GLDAS-derived TWS and CSR RL06 SH solutions with 300-km Gaussian filtering, (d) GLDAS-derived TWS and CSR RL06 SH solutions with 750-km Gaussian filtering.

Conclusions
In this study, we presented an improved method for the estimation of regional mass variations from GRACE intersatellite geopotential differences using a priori constraints. Firstly, an alternative analytic formula was presented to incorporate the KBR range-rate into the improved energy balance equation, and precise intersatellite geopotential differences were estimated from GRACE Level-1B data based on the RCR technique, which can avoid the long-wavelength gravity signals to be absorbed by empirical parameters. Secondly, GLDAS-derived TWS changes were used as a priori constraints to reduce the ill condition for solving regional mass variations from geopotential difference data, and an adaptive algorithm was employed to determine the optimal regularization parameters based on iterative least-squares estimation. Thirdly, regional mascon solutions inverted from geopotential differences over South America on 2 • × 2 • grids at monthly intervals from January 2005 to December 2010 were compared to the official GRACE mascon solutions.
The results show that the regional mascon solutions inverted from geopotential differences estimated by the RCR technique using GLDAS hydrological models as a priori constraints can retain more signal energy and enhance regional mass variation inversion, and the final estimated mascon solutions are less affected by the a priori information. The spatial distribution characteristics of mass variations over South America from the geopotential difference-based regional mascon solutions (i.e., GPD Mascon) agree well with those from the three official GRACE mascon solutions (i.e., CSR RL06M, JPL RL06M, and GSFC Mascon). Additionally, the annual amplitudes of mass change time series derived from the four mascon solutions for the eight largest basins over South America show good consistencies with each other. Aside from the seasonal variability, the basin-scale non-seasonal signals from these four mascon solutions are also in good agreement with each other. However, obvious differences can be observed among the spatial patterns and trends derived from the four mascon solutions, especially in the small basins. The discrepancies are mainly due to different inversion methods with different error approximations and data processing strategies (regularization constraints, background models, sizes of mascon blocks, and data pre-processing, etc.) applied. In addition, differences in trends are also due to their small magnitudes, as mass variations over South America basins are dominated by seasonal hydrological signals. Despite the expected discrepancies, the GPD Mascon solutions are comparable to the three official mascon solutions at the basin scale, and without needing low degrees of SH solutions to complement the missing long-wavelength signal, which demonstrates the feasibility and advantage of the proposed method.
It should be noted that our GPD Mascon solutions did not account for the leakage effects across land-ocean boundaries (like the CRI filter used in JPL RL06M), which might affect the comparisons with the official GRACE mascon solutions. However, if the basin is large or far away from the ocean, the leakage effect is expected to be small. Our future work will include the implementation of leakage corrections to further improve the accuracy of geopotential difference-based regional mascon solutions. In addition, GRACE Follow-On is equipped with laser ranging interferometer (LRI), which can provide more accurate intersatellite ranging data, so we plan to use GRACE Follow-On data to further analyze and verify the geopotential difference-based method in the future.