WHU‐GRACE‐GPD01s: A Series of Constrained Monthly Gravity Field Solutions Derived From GRACE‐Based Geopotential Differences

To suppress the correlated noise of Gravity Recovery and Climate Experiment (GRACE) spherical harmonic (SH) solutions, we developed a series of constrained monthly gravity field solutions named WHU‐GRACE‐GPD01s from August 2002 to July 2016 using GRACE‐based geopotential differences. The constrained solutions were estimated using Kaula regularization, and the optimal regularization parameters were adaptively determined from the GRACE data itself through variance component estimation. The performance of the constrained WHU‐GRACE‐GPD01s solutions was validated against the official SH solutions (GFZ, JPL, and CSR RL06) and mass concentration (mascon) solutions (CSR RL06M) at global and regional scales. The results demonstrate that mass changes derived from the constrained solutions and official SH solutions with DDK4 filtering are in good agreement, but the constrained solutions present weaker longitudinal stripes and have a lower noise level at regional scales (e.g., in the Middle Pacific Ocean and Sahara Desert). Furthermore, regional mass changes (e.g., major river basins, Greenland, and Antarctic ice sheet) inferred from the constrained solutions agree with the ones derived from the official SH solutions with DDK4 filtering. The constrained solutions also have a higher signal intensity and smaller spatial leakages as compared to the CSR RL06 SH solutions with spatial filtering (Gaussian filtering plus de‐striping). For the problematic months, the constrained solutions are more reliable than the official SH solutions using spatial filtering or DDK4 filtering, and are closer to CSR RL06M mascon solutions. These validations demonstrate that our constrained solutions are comparable to the official SH solutions and can be used without postprocessing.

Another commonly used technique to improve the signal-to-noise ratio (SNR) in estimated mass variations is regularization. Generally, a regularization filter is implemented by reformulating the filtering problem into an estimation problem with stochastic constraints (Devaraju & Sneeuw, 2017). A number of regularization filters have been proposed, such as the Wiener filter (Klees et al., 2008;Sasgen et al., 2006), DDK filter (Kusche, 2007;Kusche et al., 2009), and the standard regularization as a postprocessing step to the unconstrained solutions (Swenson & Wahr, 2011). For instance, Kusche (2007) constructed the DDK filter by using a priori signal covariance and approximate error covariance of SH coefficients; Klees et al. (2008) designed an optimal anisotropic and nonsymmetric (ANS) filter based on the expected signal in terms of SH coefficients and their error estimates; Swenson and Wahr (2011) obtained the regularized solutions by using the information available from the covariance matrix of the unconstrained solutions. The reliability of the regularization filter is based on the assumption that the external information (e.g., signal and error covariance matrices of SH coefficients) used is sufficiently accurate and reliable. Nevertheless, the signal and error covariance are not publicly available products for most of the GRACE Level-2 monthly solutions. Therefore, an approximate or empirical signal and error models are typically used, which degrade the performance of regularization filters.
In general, both spatial filtering and regularization filtering can generate comparable filtered solutions (Swenson & Wahr, 2011). However, the regularization filter can also be integrated into the processing of GRACE Level-1B data directly, which is equivalent to using regularization as a postprocessing step for unconstrained solutions (Klees et al., 2008). A major benefit of this technique is that with spatio-temporal constraints applied in the regularized least-squares estimation, we can improve signal recovery and reduce spatial leakage as well. In particular, the inversion of mass concentration (mascon) solutions makes use of spatio-temporal constraints derived from ancillary geophysical information (Loomis et al., 2019;Watkins et al., 2015) or the GRACE data itself (Save et al., 2016), thus reducing spatial leakage and enhancing spatial resolution of GRACE estimates. Essentially, the unconstrained mascon solutions are equivalent to the unconstrained SH solutions (Rowlands et al., 2010;Watkins et al., 2015), and the improvements in mascon solutions are primarily attributed to the regularization constraints applied in mascon modeling (Q. Chen et al., 2021;Ran et al., 2018Ran et al., , 2021. Therefore, the SH solutions estimated with regularization constraints are also equivalent to the mascon solutions, and the correlated noise and spatial leakage can be reduced through regularized least-squares estimation. For example, Bruinsma et al. (2010) established regularization matrices through the Kaula rule and estimated a series of 10-day gravity field solutions based on Tikhonov regularization; Save et al. (2012) used the error characteristics of unconstrained CSR RL04 monthly solutions to design constraint matrices and solved the regularized solutions using the Tikhonov regularization with L-curve method; Q. Chen et al. (2021) transformed the spatial constraints from GRACE-based ZHONG ET AL.

10.1029/2022EA002699
3 of 26 filtered mass changes into spectral domain (i.e., expressed as SH coefficients instead of mascons) and solved the regularized high-resolution (degree 180) monthly gravity field solutions based on the mean square error (MSE) method. These studies reveal that the constrained regularized solutions show markedly reduced stripes compared with the unconstrained solutions; in particular, the signal leakage is evidently suppressed for high-resolution regularized solutions. Moreover, such regularized solutions can be consistently used by end-users to estimate mass changes without further postprocessing.
In this study, we recovered constrained monthly gravity field solutions using GRACE-based geopotential difference (GPD) data determined by the energy integral method (Han, 2004;Shang et al., 2015), which is different from the methods adopted by official and other processing centers. The GRACE-based GPD is a type of pseudo-observation estimated from GRACE Level-1B RL03 data with updated background models, in which the improved energy balance equation and remove-compute-restore (RCR) technique are employed to improve the estimation accuracy (Zhong et al., 2020). In addition, Kaula regularization was used to constrain the monthly gravity field solutions to avoid introducing any bias or error in the external geophysical information; this is because Kaula's rule of thumb (Kaula, 1966) is a powerful tool for the quantitative analysis of SH spectral estimates of the gravity field without requiring any external information. Furthermore, to avoid over-smoothing caused by the L-curve and MSE methods as noted in previous studies (Q. Chen et al., 2021;Kusche & Klees, 2002;Xu, 1998), the optimal regularization parameters for the monthly solutions were adaptively determined from the GRACE data itself using variance component estimation (Koch & Kusche, 2002). Finally, a series of constrained monthly gravity field solutions named WHU-GRACE-GPD01s (complete to degree and order 96) from August 2002 to July 2016 were recovered and compared with the official SH solutions (e.g., GFZ, JPL, and CSR RL06) and mascon solutions (e.g., CSR RL06M) at global and regional scales.
The remainder of this paper is organized as follows. In Section 2, the inversion methods and data processing strategies for constrained monthly gravity field solutions based on GRACE-based GPD observations are described. In Section 3, the comparisons of constrained monthly gravity field solutions at global scale are presented. In Section 4, the regional evaluations of the constrained monthly gravity field solutions are presented. In Section 5, the advantages of regularization filtering and the reliability of the constrained solutions are discussed. Finally, the conclusions are presented in Section 6.

Estimation of GRACE-Based GPDs
Precise estimation of GRACE-based GPD is a prerequisite for the inversion of monthly gravity field solutions. In this study, we determined the residual GPD between the two GRACE satellites using an improved energy balance equation, which can be formulated in the inertial reference frame as follows (Shang et al., 2015;Zhong et al., 2020): (1) in which, 12 is the residual GPD between the two GRACE satellites; 1 , 2 , and ̇1 , ̇2 are the corresponding position and velocity vectors of the two satellites; 1 and 2 are the sums of the residual geopotential acceleration and nongravitational acceleration vectors for the two satellites; is the angular velocity vector of the Earth-fixed reference frame relative to the inertial reference frame; 12 is the kinetic energy difference between the two satellites; 12 is the reference geopotential difference; integration in Equation 1 ranges from 0 to , and 0 12 is the difference in the integral constants.
To accurately calculate 12 in Equation 1, a rigorous formula expressed by the GRACE K-band range-rate (KBRR) was adopted as follows (Zhong et al., 2020(Zhong et al., , 2023: where is the KBRR measurement, ̇1 2 is the relative velocity vector of the two GRACE satellites, 12 = 12 | 12 | is the line-of-sight (LOS) unit vector, and 12 is the relative position vector of the two GRACE satellites. ZHONG ET AL.
10.1029/2022EA002699 4 of 26 Generally, the residual GPD 12 mainly includes the time-variable components caused by mass transport and redistribution in the Earth system. The monthly gravity field solutions can be directly estimated from GRACE-based GPD observations using the least-squares method (Han, 2004) where is the geocentric gravitational constant of the Earth (3.986004415 × 10 14 m 3 s −2 ), and is the Earth's mean radius (6378136.3 m); ( 1, 1, 1) and ( 2, 2, 2) are the corresponding spherical coordinates of the two satellites in the Earth-fixed reference frame; is the fully normalized associated Legendre function of degree and order ; and are the unknown SH coefficients of degree and order .

Monthly Solutions Estimated With Kaula Regularization
As it can be seen from Equation 3, the residual GPD 12 depends linearly on the SH coefficients. To estimate unknown SH coefficients from GPD observations based on the least-squares method, the observation equation can be expressed as follows: where is the vector of the SH coefficient parameters to be estimated, is the vector of residual GPD observations, is the design matrix with respect to , is the correction vector to the GPD observations, denotes the dispersion operator, 2 0 and are the unknown variance factor and positive definite weight matrix of the observations .
For the Tikhonov regularization, the expected value and the covariance matrix of the SH coefficient parameters vector are introduced as the a priori information, which gives the observation equation as follows: where is the correction vector to the SH coefficient parameters, 0 is the a priori expected value (or mean) of the SH coefficient parameters vector, 2 and are the unknown variance factor and weight matrix of the SH coefficient parameters.
With 0 = 0 , the solution vector for the Tikhonov regularization is obtained by minimizing the functional where is also called the regularization matrix, is the regularization parameter, and = 2 0 ∕ 2 . On the left side of Equation 6, the first term forces the model to fit the data, the second term represents the regularization, and the regularization parameter balances the weights between these two terms. Subsequently, the solution of the minimization problem in Equation 6 is given bŷ For the Kaula regularization, is a diagonal matrix with the reciprocal of the degree variance (i.e., 1∕ 2 ) of the corresponding SH coefficients on its diagonal. The degree variance 2 determined using Kaula's rule of thumb is defined as (Reigber, 1989) Usually, the constant factor of 5.0 × 10 −11 in Equation 8 is omitted; its reciprocal is incorporated in the regularization parameter , and then the Kaula regularization matrix is defined as (Ditmar et al., 2003) ZHONG ET AL.
10.1029/2022EA002699 5 of 26 where and denote the row and column numbers of the Kaula regularization matrix , respectively; ( ) is the degree related to the row number ; and is the Kronecker symbol.
To obtain an optimal constrained monthly gravity field solution without any over-smoothing, we adaptively determined the regularization parameter from the GRACE data through an iterative process based on variance component estimation. The regularized solution of Equation 7 can be explicitly expressed with variance factors 2 0 and 2 as follows: Then, the unknown variance factors 2 0 and 2 are estimated iteratively, as follows: where ̂ and ̂ are estimated by 0 and are the redundancy numbers, which are computed as follows: where is the total number of GPD observations, is the number of the unknown SH coefficient parameters.
Therefore, starting with the initial variance factors, the initially regularized solution can be determined through Equation 10. The updated variance factors (Equation 11) can then be derived from the initially regularized solution with the corresponding correction vectors (using Equation 12) and redundancy numbers (using Equation 13). As a result, the variance factors and the regularized solution are iteratively determined, and the optimally regularized solution ̂ can be obtained until the variance factors 2 0 and 2 converge. For more details on the regularization of the gravity field solution using variance component estimation, please refer to the study by Koch and Kusche (2002).

Data Processing Strategies
Based on the energy integral method, GPD observations were estimated from GRACE Level-1B RL03 data, including the KBRR measurements (KBR1B), satellite orbits (GNV1B), accelerometer measurements (ACC1B), and satellite attitude measurements (SCA1B). The background models used to estimate the GPD observations are shown in Table 1, which are mainly consistent with the GRACE Level-2 Processing Standards Document for Level-2 Product RL06 (Bettadpur, 2018;Dahle et al., 2018;Yuan, 2018). For example, the static earth gravity filed was modeled by the GGM05C model (Ries et al., 2016), and the third-body perturbations were computed 6 of 26 based on JPL planetary and lunar ephemeris DE430 (Folkner et al., 2014); the solid earth tide, solid earth pole tide, ocean pole tide (Desai, 2002), as well as the contribution of general relativity were computed according to International Earth Rotation Service (IERS) Conventions 2010 (Petit & Luzum, 2010). The ocean tide was modeled by the EOT11 model (Savcenko & Bosch, 2012), and the S1 and S2 air tidal contributions were calculated by the Ray/Ponte model (Ray & Ponte, 2003). Additionally, the Atmosphere and Ocean De-aliasing Level-1B (AOD1B) RL06 products were used to reduce the short period nontidal variability in the atmosphere and oceans.
To reduce the impact of satellite orbit and background model errors on GPD estimation, we estimated purely dynamic orbits for each day by using the GNV1B orbits as pseudo-observations and calibrated accelerometer data simultaneously based on the dynamic integral approach using the known background models. Subsequently, based on the improved energy balance equation (Equation 1), the residual GPD observations can be calculated by using the KBR1B data, re-estimated purely dynamic orbits, calibrated accelerometer data, SCA1B data, and the known background models. This implies that no new time-variable gravity information can be obtained from Equation 1 if the KBRR measurements were absent. The advantage of the above calculation is that the estimated residual GPD only includes the time-variable gravity signals from the KBRR data, which effectively reduces the error influences of the satellite orbit and background model (Luthcke et al., 2006;Shang et al., 2015;Zhong et al., 2020).
Note that the residual GPD observations estimated using Equation 1 contain many systematic errors (mainly caused by the KBRR low-frequency noises), which are conventionally reduced by applying empirical parameters (Han et al., 2006;Zhao et al., 2011;Zhou et al., 2018). However, the empirical parameters not only reduce systematic errors but also absorb long-wavelength parts of temporal gravity signals (Shang et al., 2015;Zhong et al., 2020). To retain the time-variable gravity information to the full extent, we employed RCR technique to estimate residual GPD observations. Our previous study (Zhong et al., 2020(Zhong et al., , 2023 demonstrated that the RCR technique can preserve both low-frequency and high-frequency temporal gravity signals, and the accuracies of GPD-inferred regional mascon solutions in South America and the Yangtze River Basin of China are comparable to the official mascon solutions (e.g., CSR and JPL RL06 mascons). Therefore, based on the accurately estimated GPD observations, the inversion of time-variable gravity field models at a full spectrum (both low-frequency and high-frequency) can be obtained. It should be pointed out that the RCR technique employs the low-degree SH solutions (i.e., CSR RL06) as the reference field to restore the low-frequency temporal gravity signals. Therefore, the monthly gravity field solutions recovered by the GPD observations in this study are not strictly independent products, which are different from the SH solutions developed based on dynamic approach. For more details on the precise estimation of the GRACE-based GPD observations, please refer to the study by Zhong et al. (2020). In addition, the data set of GRACE-based residual GPD observations from April 2002 to July 2016 is available from Zhong et al. (2022). Figure 1 illustrates the residual GPD observations in March, June, September, and December 2005 estimated from GRACE Level-1B RL03 data based on the improved energy balance equation and RCR technique, in which the GGM05C model up to degree and order 360 was used as the reference gravity field. It can be seen that the GRACE-based GPD, as a quantity with explicit geophysical interpretation (compared with the geometric KBRR observation), is a direct reflection of the time-variable gravity field and has evident and characteristic seasonal and regional variations (see Figures 1a-1d). As for all estimates from April 2002 to July 2016, we found that the root mean square (RMS) of the residual GPD for each month was ∼2.0-3.0 ×10 −3 m 2 /s 2 . Therefore, to reduce the influence of outliers on the estimation of monthly gravity field models, we screened out the GPD observations with the values greater than 0.01 m 2 /s 2 (slightly greater than 9.0 ×10 −3 m 2 /s 2 ) based on the three times standard deviation criterion during the estimation process. Additionally, to accelerate the convergence speed of the regularized solution of Equation 7 based on variance component estimation, we set the initial variance factors 2 0 = 1.0 × 10 −6 and 2 = 1.0 in this study. Our tests indicate that if the convergence condition is defined as the absolute value of the difference between the variance factors calculated by two adjacent iterations being less than 1.0× 10 −12 , then the regularized solutions can converge within 10 iterations for most of the months.

Global Comparisons of Constrained Monthly Solutions
Based on the inversion method and data processing strategies presented in Section 2, a series of constrained monthly gravity field solutions (named WHU-GRACE-GPD01s, published at: http://icgem.gfz-potsdam.de/ series/03_other/WHU/WHU-GRACE-GPD01s) complete to degree and order 96 from August 2002 to July 2016 were recovered using GRACE-based GPD observations. In this section, we verify the performance of the WHU-GRACE-GPD01s solutions by comparing them with the official SH solutions (including GFZ, JPL, and CSR RL06), especially for inferring surface mass changes at the global scale. It should be noted that we only considered the period from January 2005 to December 2010 to avoid problematic months (e.g., reduced observation number, deep orbit repeat cycles, and poor observation accuracy) distorting the mass change estimates. Therefore, the homogeneous data quality during this period allowed us to conduct a fair comparison between the WHU-GRACE-GPD01s solutions and official SH solutions. In addition, to infer reliable mass changes from the official SH solutions, we used the published SH solutions with DDK4 filtering for comparisons. This is because the DDK filter is a reliable regularization filter, which can well reduce the signal leakages and improve the SNR of SH solutions compared to the commonly used spatial filtering (e.g., Gaussian filtering).

Spatial Patterns of Global Mass Changes
Investigating the spatial patterns of the global mass changes allows us for better understanding of the performance of WHU-GRACE-GPD01s solutions. We calculated global mass changes in terms of equivalent water heights (EWH) at 1° × 1° grid based on the constrained SH solutions with maximum degree and order 96. Meanwhile, the global mass changes derived from the official unconstrained SH solutions with DDK4 filtering were also presented for a comparison. It should be noted that while calculating mass changes in this study (unless otherwise mentioned), the C 20 coefficients of the official SH solutions and WHU-GRACE-GPD01s SH solutions were replaced by the SLR estimates (Loomis et al., 2020), and the geocenter motion and Glacial Isostatic Adjustment (GIA) effects were corrected using the degree-1 coefficients from Landerer (2019) and the ICE-6G_D model (Peltier et al., 2018), respectively. Figure 2 presents the spatial patterns of global mass changes derived from the three official unconstrained SH solutions (GFZ, JPL, and CSR RL06 with DDK4 filtering) and the WHU-GRACE-GPD01s solutions for April and August 2005. As shown in Figure 2, the temporal signal features over South America, South Africa, and Greenland can be well recovered by these four SH solutions and the signal amplitudes are comparable to each other. However, although the correlated noises in the official SH solutions have been reduced through DDK4 filtering, there are still residual longitudinal stripe errors in the mid-low latitudes. In particular, the longitudinal stripes of the GFZ RL06 and JPL RL06 solutions were more significant than those of the CSR RL06 solutions. On the contrary, WHU-GRACE-GPD01s solutions present almost no or very small longitudinal stripes. This demonstrates that our constrained SH solutions with Kaula regularization can effectively suppress the longitudinal stripe errors compared to the official unconstrained SH solutions with DDK4 filtering as a postprocessing strategy.
To assess the performance of different SH solutions for inferring global mass changes on the long time scales, the annual amplitudes, annual phases, and long-term trends of mass changes at the global scale from January 2005 to December 2010 were calculated and presented in Figure 3. The spatial patterns of annual amplitudes from WHU-GRACE-GPD01s solutions are similar to those of official unconstrained SH solutions with DDK4 filtering, and the evident seasonal characteristics of mass changes in South America, Southeast Asia, South Africa, as well as other typical regions are well revealed. For the annual phases, the WHU-GRACE-GPD01s solutions and official SH solutions still present similar characteristics and the annual phases of the ocean areas lagged significantly behind the land areas. For the long-term trends, both WHU-GRACE-GPD01s solutions and official SH solutions detect significant trend changes, especially in the regions of South Africa, Greenland, and Antarctic ice sheet. Nevertheless, there are still obvious differences in the regional areas (e.g., in the Greenland and Antarctic ice sheet), more detailed comparisons will be presented in the regional evaluations.
Additionally, the annual amplitudes, annual phases, and long-term trends of global mass changes from the official unconstrained SH solutions with DDK4 filtering are still contaminated by the residual correlated noises (see Figure 3). Among them, the annual phases of the official SH solutions are more contaminated by longitudinal stripes than the annual amplitudes and long-term trends. In contrast, the constrained WHU-GRACE-GPD01s solutions can better suppress the longitudinal stripe errors than the official SH solutions with DDK4 filtering (see the RMS residuals statistics in Section 3.2), which preliminarily verifies the effectiveness of our regularized monthly gravity field solutions.

Quantification of Noise Level at Global Scale
In general, surface mass variations (mainly including terrestrial water, ice sheets, and ocean mass changes) in land areas are more remarkable than those in ocean areas (see the annual amplitudes of global mass changes presented in Figure 3). A commonly used method to evaluate the uncertainty of GRACE-inferred surface mass changes is to use RMS residuals over the oceans, where true long-term mass variations are expected to be very small (J. L. Kvas et al., 2019). To this end, we simultaneously fitted the bias, trend, annual and semiannual components to the mass change time series derived from the official unconstrained SH solutions (GFZ, JPL, and CSR RL06 with DDK4 filtering) and WHU-GRACE-GPD01s solutions, and then obtained the RMS residuals by removing the fitted signals from the mass change time series. Figure 4 presents the spatial patterns of the RMS residuals of global mass changes derived from different monthly gravity field solutions, i.e., the three official unconstrained SH solutions with DDK4 filtering and the constrained WHU-GRACE-GPD01s solutions. As shown in Figure 4, the GFZ RL06 and JPL RL06 solutions demonstrated comparable levels of RMS residuals over ocean areas. Nevertheless, the magnitudes of RMS residuals for constrained WHU-GRACE-GPD01s solutions are comparable to those of the CSR RL06 solutions, which are smaller than those of the GFZ RL06 and JPL RL06 estimates. The results demonstrate that our constrained monthly gravity field solutions (WHU-GRACE-GPD01s) are more reliable than the GFZ and ZHONG ET AL.
10.1029/2022EA002699 9 of 26 JPL RL06 solutions using DDK4 filtering and present comparable or even better performance than the CSR RL06 SH solutions with DDK4 filtering. It should be noted that the magnitude of the RMS residuals over the Amazon basin is larger than other regions. This is because the main components of the surface mass changes are terrestrial water and ice sheets distributed in the land area, the Amazon basin possesses the most significant

Regional Evaluations of Constrained Monthly Solutions
Although the global comparison results show that the WHU-GRACE-GPD01s solutions are consistent with or even better than the official unconstrained SH solutions with DDK4 filtering in the spatio-temporal domains, the performance of our constrained solutions in recovering mass changes at regional scales needs to be further evaluated. For this purpose, we calculated the time series and spatial patterns of mass changes in typical regions (including major river basins, deserts, ocean areas, Greenland, and Antarctic ice sheet) from the WHU-GRACE-GPD01s solutions, and compared them with the corresponding estimates from the official SH solutions. Figure 6 presents the time series of mass changes for the four representative river basins (i.e., Amazon, Congo, Mississippi, and Yangtze) derived from the official unconstrained SH solutions with DDK4 filtering and the constrained WHU-GRACE-GPD01s solutions. As shown in Figure 6, the mass changes derived from the WHU-GRACE-GPD01s and official SH solutions are in good agreement, and both well reveal evident seasonal variations among these four river basins. It is noteworthy that there are still some differences between these four solutions in individual months over the four river basins due to the different inversion method and data processing strategies. Table 2 lists the RMS differences and correlation coefficients (CC) of mass change time series between the WHU-GRACE-GPD01s and official SH solutions in above four river basins. It demonstrates that the RMS differences and CC between the WHU-GRACE-GPD01s and official SH solutions are 0.19-0.91 cm and 0.983-0.999, respectively. The RMS differences and CC between the WHU-GRACE-GPD01s and CSR RL06 solutions are better than those between the WHU-GRACE-GPD01s and GFZ RL06, WHU-GRACE-GPD01s and JPL RL06 estimates. This means that our constrained WHU-GRACE-GPD01s solutions are more consistent with the CSR RL06 SH solutions with DDK4 filtering. One major reason for this is that, as mentioned earlier, the GPD observations estimated by RCR technique use the low-degree CSR RL06 SH solutions as the reference field, so the recovered WHU-GRACE-GPD01s solutions are close to CSR RL06 SH solutions in the long-wavelength components.

Mass Changes at River Basin Scales
In addition, Table 3 lists the long-term trends, annual amplitudes, and corresponding SNRs (the computation method can be found in Zhou et al. (2019)) of the mass changes derived from the official SH solutions with DDK4 filtering and WHU-GRACE-GPD01s solutions over 20 major river basins. As listed in Table 3, the longterm trends of these river basins estimated by the constrained WHU-GRACE-GPD01s solutions are agree with the official SH solutions. Additionally, compared with the official SH solutions, the WHU-GRACE-GPD01s has comparable or stronger (∼35% of the river basins) annual amplitudes over these river basins. Meanwhile, the SNR values of the WHU-GRACE-GPD01s solutions are comparable to or even better (∼65% of the river basins) than those of official SH solutions. Hence, these results demonstrate the stability and reliability of our constrained monthly gravity field solutions for inferring mass changes at the river basin scales.

Mass Changes Over Greenland and Antarctic Ice Sheet
The Greenland and Antarctic ice sheet suffer significant mass loss from land-based glaciers due to human-induced global warming. The Greenland is the largest non-continental island and the second largest ice sheet in the world, with a total area of 2,166,086 km 2 . The Antarctic ice sheet is the largest single mass of ice on Earth, covering an area of ∼14 million km 2 and containing 30 million km 3 of ice (Rémy & Parouty, 2009). In the context of global climate change, the melting of Greenland and Antarctic ice sheet has become an important contributor to sea-level rise (J. L. Chen et al., 2006;Sasgen et al., 2012). Accurate quantification of the spatio-temporal distributions of ice sheets and their increasing and decreasing rates is of great significance for dealing with sea-level and global climate changes.    Figure 7, for the Greenland, the long-term trends of mass changes derived from these four solutions show consistent spatial patterns, a significant decreasing trend of mass changes in the west and south but a certain increasing trend of mass changes in the northeast of Greenland can be observed. In addition, the strength of the trend signals derived from the WHU-GRACE-GPD01s solutions are slightly smaller than those of official SH solutions. This is mainly owing to the smoothing effect of Kaula regularization, which will be discussed later. For the Antarctic ice sheet, the spatial patterns of mass change trends derived from the four solutions are in good agreement, and they present an evident decreasing trend in mass changes in the Antarctic Peninsula and West Antarctica, but with an increasing trend of mass changes over most parts of East Antarctica. However, it is worth noticing that both the official SH solutions (with DDK4 filtering) and constrained WHU-GRACE-GPD01s solutions still present signal leakages at the edge of Greenland and Antarctic ice sheet.    Note. Note that the three numbers separated by the slashes correspond the long-term trend, annual amplitude, and SNR, respectively; the bold numbers indicate that WHU-GRACE-GPD01s has larger long-term trends, stronger annual amplitudes, or higher SNRs than all the three official SH solutions.
ZHONG ET AL.

10.1029/2022EA002699
14 of 26 and the GFZ RL06 solutions (−0.38 ± 0.11 cm/yr) are the smallest. Again, these comparisons also suggest that the constrained WHU-GRACE-GPD01s solutions are comparable to the official unconstrained SH solutions with DDK4 filtering, and further demonstrate the effectiveness of Kaula regularization filtering.

Quantification of Noise Level at Regional Scale
Another approach to quantify the performance of the monthly gravity field solutions is to compare the surface mass changes over the regions with minimal signal amplitude (Zhou et al., 2018). Figure 9 shows the mass change time series in the Middle Pacific Ocean and Sahara Desert, where the temporal gravity signals are expected to be very small.
As shown in Figure 9, the mass changes in the Middle Pacific Ocean and Sahara Desert derived from the official SH solutions with DDK4 filtering and WHU-GRACE-GPD01s solutions are in good agreement, but there are still evident differences when compared month-by-month. Specifically, the RMS residuals in the Middle Pacific Ocean are 1.68 cm (GFZ RL06), 1.32 cm (JPL RL06), 1.44 cm (CSR RL06), and 1.36 cm (WHU-GRACE-GPD01s), and the corresponding values in the Sahara Desert are 1.76 cm (GFZ RL06), 1.94 cm (JPL RL06), 1.80 cm (CSR RL06), and 1.76 cm (WHU-GRACE-GPD01s), respectively. This demonstrates that our constrained WHU-GRACE-GPD01s solutions generally suffered from lower noise levels than the official unconstrained SH solutions with DDK4 filtering during the study period (i.e., from January 2005 to December 2010).

Comparison of Regularization and Spatial Filtering
From the global comparisons and regional evaluations presented in Sections 3 and 4, we found that mass changes derived from the constrained WHU-GRACE-GPD01s solutions presented good consistency with the official unconstrained SH solutions using DDK4 filtering in the spatiotemporal  domains. Note that the DDK filter is also a regularization filter, which uses the regularization as a postprocessing step for the unconstrained SH solutions. Unlike the DDK filtering, our regularization filtering is integrated into the processing of GRACE Level-1B data in regularized least-squares estimation. Theoretically, these two regularization filters are equivalent. The differences between their filtered SH solutions are mainly related to the selection of regularization constraint matrix and the determination of regularization parameters.
To reduce the correlated noise, another classical method is the spatial filtering. For instance, the commonly used spatial filtering is Gaussian filtering plus de-striping (combined spatial filtering) (Swenson & Wahr, 2006;Wahr et al., 1998), which has been widely used as the postprocessing of the official unconstrained SH solutions. Compared with spatial filtering, regularization filtering can not only effectively suppress correlated noises, but also help to improve signal recovery and reduce spatial leakage. To illustrate this point, we compared the constrained WHU-GRACE-GPD01s solutions with the unconstrained GPD-inverted monthly gravity field solutions (called unconstrained WHU-GRACE-GPD01s with the maximum degree and order 96) using the combined spatial filtering. Figure 10 shows the spatial patterns of annual amplitudes, annual phases, and long-term trends of global mass changes from the constrained (top row) and unconstrained (bottom row) GPD-inverted monthly SH solutions over the period of January 2003 to December 2013. Note that the unconstrained WHU-GRACE-GPD01s solutions were postprocessed by 300-km Gaussian filtering and P4M6 de-striping (J. L. Chen et al., 2009).
As shown in Figure 10, the annual amplitudes, annual phases, and long-term trends of surface mass changes derived from the constrained and unconstrained WHU-GRACE-GPD01s solutions are generally in good agreement. However, although both constrained and unconstrained solutions are estimated from the same GPD observations, there are still significant differences between them due to different filtering processing strategies (i.e., regularization filtering versus spatial filtering). For instance, the signal amplitudes of unconstrained solutions using the combined spatial filtering are obviously weaker than those of constrained solutions (e.g., in the Amazon basin, South Asia, and central Africa), because the spatial filtering can reduce both the correlated noise and temporal signals. For the annual phases, the constrained and unconstrained solutions have obvious differences in the ocean areas (e.g., the western part of the Atlantic, and Southwestern part of the Indian ocean) due to different filtering strategies. Furthermore, the long-term trends derived from the unconstrained solutions using the combined spatial filtering suffer from larger spatial leakages than those of constrained solutions.
To compare the regularization filtering and spatial filtering on the regional scales, we further investigated the performance of the constrained WHU-GRACE-GPD01s solutions (corresponding to regularization filtering) and the unconstrained CSR RL06 SH solutions (with combined spatial filtering) in recovering mass changes over Antarctic ice sheet and Amazon basin. Here, to effectively reduce the correlated noise, we also simultaneously performed 300-km Gaussian filtering and P4M6 de-striping on the unconstrained CSR RL06 SH solutions as the combined spatial filtering (unless otherwise mentioned). Figure 11 shows the patterns of mass change trends over Antarctic ice sheet derived from the CSR RL06 SH solutions with combined spatial filtering and constrained WHU-GRACE-GPD01s  As shown in Figure 11, the strength of the trend signals over the Antarctic ice sheet derived from CSR RL06 with combined spatial filtering (Figure 11a) are significantly smaller than those of other three solutions (Figures 11b-11d), and obvious signal leakages can be observed. The spatial patterns of mass change trends derived from the CSR RL06 SH solutions with DDK4 filtering and constrained WHU-GRACE-GPD01s solutions are highly consistent with each other, and they are closer to the CSR RL06M mascon solutions than the CSR RL06 SH solutions with combined spatial filtering. This is because the spatial filtering reduces both the temporal signals and the high-frequency noises, while the regularization filtering can improve the recovery of temporal signals and suffer from lower leakage errors. However, the spatial resolution of the CSR RL06M estimates was significantly higher than those of the constrained WHU-GRACE-GPD01s solutions and CSR RL06 SH solutions with DDK4 filtering, with almost no spatial leakage. This is because the CSR RL06M mascon solutions were also estimated through regularization filtering, but with higher spatial resolution (an equal-area geodesic grid of size ∼1° × 1° at the equator (see Save et al., 2016). In addition, there is a narrow positive trend signal in western Antarctica derived from CSR RL06M mascon solutions (see Figure 11d), but this trend signal is very significant in Figures 11a-11c, which may be the signal leakage caused by spatial or regularization filtering.
The Amazon basin is the largest tropical rainforest area on the earth that contains the Amazon river (the second longest river in the world) and its tributaries. The seasonal mass changes over the Amazon basin were more remarkable than those in other regions (see the left column of Figure 3), therefore investigating the spatial patterns and signal amplitudes of the mass changes over the Amazon basin can well verify the performance of temporal gravity field solutions. Figure 12 presents the spatial patterns of annual amplitudes of mass changes for the Amazon basin derived from the CSR RL06 SH solutions with combined spatial filtering, CSR RL06 SH solu tions with DDK4 filtering, constrained WHU-GRACE-GPD01s solutions, and CSR RL06M mascon solutions.
As shown in Figure 12, the patterns of annual amplitudes derived from the four solutions are consistent with each other, but there are still some differences in the signal strengths. We further calculated the mean annual amplitudes of Amazon basin from these four solutions, the mean annual amplitude of the CSR RL06 SH solutions with combined spatial filtering is the smallest (21.89 ± 2.03 cm); although the signal strength of WHU-GRACE-GPD01s (23.74 ± 2.18 cm) is smaller than that of CSR RL06 SH solutions with DDK4 filtering (24.59± 2.35 cm), the WHU-GRACE-GPD01s estimates are closer to the CSR RL06M mascon solutions (24.12 ± 2.02 cm), and the differences of mean annual amplitudes between CSR RL06M and WHU-GRACE-GPD01s, CSR RL06M and CSR RL06 SH solutions with DDK4 filtering are 0.38 and 0.47 cm, respectively. This further demonstrates that the WHU-GRACE-GPD01s solutions solved with regularization filtering improve the signal recovery better than the unconstrained CSR RL06 SH solutions with combined spatial filtering, and are closer to the CSR RL06M mascon solutions than the CSR RL06 SH solutions with DDK4 filtering.

Effectiveness of Kaula Regularization
In this study, Kaula regularization was employed to constrain the WHU-GRACE-GPD01s solutions and variance component estimation was used to adaptively determine the optimal regularization parameters from the GRACE data itself. The results demonstrate that, due to the Kaula regularization filtering, mass changes derived from the constrained WHU-GRACE-GPD01s solutions generally suffered from lower noise levels when compared with the official unconstrained SH solutions using DDK4 filtering (see Figures 4, 5, and 9). Compared with the DDK filtering, Kaula regularization filtering effectively reduces the high-frequency noises, but it also suppresses the 18 of 26 signal strength of mass changes, especially in the regions with strong temporal gravity signals (e.g., the Greenland and individual river basins (such as the Amazon and Zambezi rivers), see Figure 7, and Tables 3 and 4). This is because the Kaula regularization leads to a globally average damping of the estimated gravity field, especially the regions with large mass change signals would be suppressed to a certain extent. However, previous analysis demonstrates that the temporal signal attenuations caused by Kaula regularization are not evident (see Tables 3  and 4).
To illustrate this point, we conducted a numerical simulation experiment to quantify the impact of Kaula regularization filtering on the temporal gravity signal. We used the Global Land Data Assimilation System (GLDAS; Rodell et al., 2004) hydrological model outputs as the original temporal signal, and simulated the residual GPD observations over the period from January to December 2005 based on the real GRACE orbits. We estimated two sets of monthly gravity field solutions (with maximum degree and order 96) using the simulated GPD data without any observation errors. One set of unconstrained SH solutions were estimated directly, and the other set of constrained SH solutions were estimated through Kaula regularization. Note that the differences between the estimated solution without observation errors and the original GLDAS signal can be attributed to the errors of inversion model and numerical calculation, which is helpful for quantifying the impact of regularization filtering and spatial filtering on the original signal. Figure 13 shows the mass changes derived from the original GLDAS signal, unconstrained and Kaula constrained solutions without observation errors in January and September 2005, respectively. It is noteworthy that the

10.1029/2022EA002699
19 of 26 unconstrained solutions are postprocessed by the combined spatial filtering (as described previously), and the RMS residuals for each solution are also presented for comparison. We can see that both the unconstrained and Kaula constrained solutions estimated from the GPD data without observation errors can well recover the original GLDAS signals, but the Kaula constrained solutions present stronger temporal signals and smaller spatial leakages compared with the unconstrained solutions using combined spatial filtering, and the RMS residuals of the unconstrained and Kaula constrained solutions are 1.48 and 0.63 cm in January 2005 and 1.81 and 0.97 cm in September 2005, respectively. The results demonstrate that the Kaula constrained solutions are closer to the original GLDAS signals than the unconstrained solutions with combined spatial filtering, especially the Kaula constrained solutions estimated from the GPD data without observation errors show almost no signal attenuations and spatial leakages, which verifies the effectiveness of Kaula regularization. Figure 14 presents the spatial patterns of annual amplitudes, annual phases, and long-term trends of mass changes from the original GLDAS signals, unconstrained and Kaula constrained solutions estimated from GPD data without observation errors for the period of January-December 2005, in which the unconstrained solutions were postprocessed by the combined spatial filtering. As shown in Figure 14, the annual amplitudes, annual phases, and long-term trends of mass changes inferred from the Kaula constrained solutions are more consistent with those of the original GLDAS signals. Again, the Kaula constrained solutions present smaller To further quantify the attenuation or distortion effect of Kaula regularization filtering on the temporal gravity signal, we compared the annual amplitudes of mass changes inferred from different solutions in the river basins with large mass change signals. Table 5 lists the annual amplitudes of mass changes of 12 river basins (with annual amplitudes greater than 10 cm) inferred from the original GLDAS signals, unconstrained and Kaula constrained solutions (estimated from the GPD data without observation errors) from January to December 2005. As listed in Table 5, the annual amplitudes of mass changes of these river basins estimated by the unconstrained and Kaula constrained SH solutions are generally weakened, except for individual basins (e.g., the Chari river basin). However, the annual amplitudes of Kaula constrained solutions are closer to the original GLDAS signals than those of unconstrained solutions using combined spatial filtering. For instance, the deviation rates of annual amplitudes for Kaula constrained solutions in the 12 river basins range from 0.86% (Mekong) to 9.22% (Volta), while the corresponding values for the unconstrained solutions range from 3.41% (Volta) to 31.61% (Xi). This demonstrates that the signal attenuations or distortions caused by Kaula regularization filtering are much smaller than those caused by spatial filtering. Therefore, we believe that the Kaula regularization is applicable and effective to the recovery of temporal gravity field.

Reliability of Constrained Monthly Solutions
The precision of GRACE monthly gravity field solutions is not only related to the inversion methods and data processing strategies, but is also affected by the accuracy of observation data, orbit repeat cycles, and observation number (Q. Chen et al., 2021;Kvas et al., 2019;Shang et al., 2015). When the data have poor observation accuracy, reduced observation number, and poor data coverage, the stronger constraints (i.e., larger regularization parameters) are required to stabilize the inversion results. To verify the reliability of the constrained WHU-GRACE-GPD01s solutions estimated by poor GPD observations, we analyzed the variance factors of GPD observations and priori information (Kaula constraints), as well as the data availability of GPD observations (with outliers removed) and corresponding regularization parameters for each month. Furthermore, we also compared the spatial patterns of mass changes derived from the constrained WHU-GRACE-GPD01s solutions, CSR RL06 SH solutions with combined spatial filtering, CSR RL06 SH solutions with DDK4 filtering, and CSR RL06M mascon solutions, especially in the representative months that were affected by problematic data (e.g., reduced observation number, deep orbit repeat cycles, and poor observation accuracy). Figure 15 presents (a) the prior variance factor time series of the raw GPD residuals and the corresponding time series with outliers are screened out, (b) the data availability for each month after the outliers are removed, (c) the posterior variance factors 2 0 (for GPD observations) and 2 (for Kaula constraints) estimated along with the constrained WHU-GRACE-GPD01s solutions, and (d) the corresponding regularization parameters ( = 2 0 ∕ 2 ), over the period of January 2003 to December 2013. As shown in Figure 15a, we can see that the accuracy of raw GPD observations (black curve) between 2005 and 2010 is significantly better than those of other time periods, and the variance is at the level of 10 −4 . However, for the variance of the GPD residuals (red curve) with the outliers is removed (see Section 2.3), the statistical accuracy of GPD before 2005 and after 2010 is slightly better than that between 2005 and 2010 (the accuracy difference between different periods is actually very small), and the variance is improved to the level of 10 −6 . This seems to be contrary to the statement that the quality of GRACE data was better between 2005 and 2010. This is because the outliers with values greater than 0.01 m 2 /s 2 have been screened out for statistical GPD residuals. However, for the data quality of each month, in addition to the observation accuracy, it also depends on the data availability and data coverage. Although the statistical data accuracy between 2005 and 2010 is slightly lower than those periods before 2005 and after 2010 after the outliers are removed, the data availability (see Figure 15b) between 2005 and 2010 is significantly higher than those of other periods (especially before 2004 and after 2010). Benefitting from the larger amount of available GPD data between 2005 and 2010 than in other time periods, the estimation accuracy of the final SH solutions can be improved. Therefore, we believe that the overall data quality for the recovery of monthly gravity field models between 2005 and 2010 is still better than those of other periods.
As shown in Figure 15c, the posterior variance factors 2 0 and 2 for the GPD observations and Kaula constraints generally present opposite patterns, which indicates the relative contributions of the GPD data and priori information to the final SH solutions. The mean value of the posterior variance factor 2 0 over the period of January 2003 to December 2013 is ∼2.2 × 10 −6 , and the corresponding posterior standard deviation of the GPD data is ∼1.5 × 10 −3 m 2 /s 2 , which is less than the a priori RMS of residual GPD observations 2.0-3.0 ×10 −3 m 2 /s 2 . It should be noted that the shape of the monthly prior variance factors (the red curve in Figure 15a) and posterior variance factors (the red curve in Figure 15c) of GPD residuals is highly similar to each other, this is because the balance between the observation data and priori information is mainly determined by adjusting the variance factor of the Kaula constraints (the black curve in Figure 15c) when the prior variance of the observation data is relatively accurate. As shown in Figure 15d, the time series of the regularization parameters shows features similar to the posterior variance factors of the GPD data (see the red curve in Figure 15c), and they are highly correlated with a CC of 0.77. Although many factors (e.g., observation accuracy, number of observations, data coverage, etc.) can affect the determination of regularization parameter, here the results demonstrate that the overall observation accuracy of GPD data (with outliers removed) is the main factor affecting the determination of regularization parameter. Also, the regularization parameter = 2 0 ∕ 2 shows that is proportional to 2 0 , which means that higher observation accuracy usually corresponds to a smaller regularization parameter. It is important to note that the regularization parameters were determined adaptively from the data itself based on the variance component estimation in this study, no external interference has been introduced into the calculation process. Therefore, we believe that these estimated regularization parameters are the optimal values determined by the adaptive estimation under the consideration of the quality of GPD data (including the observation accuracy, ZHONG ET AL.  16a, 16d, 16g, and 16j), all the four solutions effectively recover the mass changes at large and medium scales. However, the CSR RL06 SH solutions with combined spatial filtering present more evident longitudinal stripes at low latitudes and the temporal signals are also attenuated more seriously due to the postprocessing smooth filtering. In addition, the residual longitudinal stripes of constrained WHU-GRACE-GPD01s solutions are smaller than those of the CSR RL06 SH solutions with DDK4 filtering, and the WHU-GRACE-GPD01s solutions are closer to the CSR RL06M mascon solutions.
Additionally, for the months with deep orbit repeat cycles (see Figures 16b, 16e, 16h, and 16k), although the CSR RL06 SH solutions have been processed by combined spatial filtering or DDK4 filtering, there are still evident longitudinal stripe errors between mid-latitude and low latitude; on the contrary, the constrained WHU-GRACE-GPD01s solutions are less contaminated by longitudinal stripes, and are closer to the CSR RL06M estimates than those CSR RL06 SH solutions. For months with poor observation accuracy (see Figures 16c, 16f, 16i, and 16l), the signal amplitudes of mass changes derived from CSR RL06 SH solutions with combined spatial filtering are significantly smaller than the other three solutions. The signal amplitudes of the WHU-GRACE-GPD01s solutions are slightly smaller than the CSR RL06 SH solutions with DDK4 filtering, but the longitudinal stripe errors of constrained WHU-GRACE-GPD01s solutions are smaller. In addition, the CSR RL06M mascon solutions present the largest signal amplitudes and less error contamination than the three SH solutions. The above results demonstrate that even for months with problematic data, mass changes derived from the constrained WHU-GRACE-GPD01s solutions are generally more reliable than the official unconstrained SH solutions (e.g., CSR RL06) with combined filtering or DDK4 filtering, especially for the months with deep orbit repeat cycles.

Conclusions
In this study, we developed a series of constrained monthly gravity field solutions named WHU-GRACE-GPD01s (complete to degree and order 96) from August 2002 to July 2016 by using GRACE-based GPD observations, which were accurately estimated from the latest GRACE Level-1B (RL03) data using updated background models based on the improved energy balance equation and RCR technique. To suppress the correlated noises of SH coefficients, the WHU-GRACE-GPD01s solutions were estimated using Kaula regularization constraints, during which the optimal regularization parameters were adaptively determined from GRACE data itself using variance component estimation. The effectiveness of Kaula regularization for the recovery of temporal gravity field was verified through a numerical simulation experiment, and the performance of the constrained WHU-GRACE-GPD01s solutions was verified by comparisons with the official unconstrained SH solutions (GFZ, JPL, and CSR RL06) and mascon solutions (CSR RL06M) at global and regional scales. At the global scale, the spatial patterns of annual amplitudes, annual phases, and long-term trends of global mass changes derived from the WHU-GRACE-GPD01s solutions and official SH solutions with DDK4 filtering are highly consistent with each other. However, the WHU-GRACE-GPD01s solutions are less contaminated by longitudinal stripes than the official SH solutions. Furthermore, the RMS residuals of the WHU-GRACE-GPD01s solutions are significantly smaller than those of official SH solutions with DDK4 filtering, and the mean values of the RMS residual time series derived from the GFZ RL06, JPL RL06, CSR RL06, and WHU-GRACE-GPD01s solutions over the ocean areas are 2.31, 2.37, 2.15, and 1.73 cm, respectively. At the regional scale, the constrained WHU-GRACE-GPD01s solutions revealed evident seasonal variation of mass changes over 20 major river basins and demonstrated good consistency with the official SH solutions with DDK4 filtering. However, the WHU-GRACE-GPD01s solutions generally have comparable annual amplitudes and better SNRs (∼65% of the river basins) compared with the three official SH solutions. In addition, the spatial patterns and time series of mass changes derived from WHU-GRACE-GPD01s solutions over typical regions (e.g., Greenland and Antarctic ice sheet, Amazon basin) are consistent with the official SH solutions using DDK4 filtering, and are closer to the CSR RL06M mascon solutions than the CSR RL06 SH solutions with combined spatial filtering (i.e., Gaussian filtering plus de-striping). Furthermore, the constrained WHU-GRACE-GPD01s solutions generally suffer from lower noise levels than the official SH solutions with DDK4 filtering at the regional scales. For example, the RMS residuals in the Middle Pacific Ocean are 1.68 cm (GFZ RL06), 1.32 cm (JPL RL06), 1.44 cm (CSR RL06), and 1.36 cm (WHU-GRACE-GPD01s), and the corresponding values in the Sahara Desert are 1.76 cm (GFZ RL06), 1.94 cm (JPL RL06), 1.80 cm (CSR RL06), and 1.76 cm (WHU-GRACE-GPD01s). Moreover, for the problematic months, mass changes derived from the constrained WHU-GRACE-GPD01s solutions (corresponding to regularization filtering) are more reliable than the official SH solutions (e.g., CSR RL06) with spatial filtering or DDK4 filtering.
Aliasing and spatial leakage are still an important issue when using GRACE data because necessitate selection of appropriate filters to maintain a good spatial resolution, while removing as much of the correlated noises as possible in the SH solutions. The design criteria for most postprocessing filtering are empirical or time-dependent and model-dependent, resulting in different filtered solutions even when using SH solutions with the same filtering technique. In contrast, regularization filtering is integrated into the least-squares estimation process, which can improve temporal signal recovery and reduce leakage errors. Furthermore, regularized SH solutions can be directly used to infer mass changes without postprocessing, which ensures the consistency of the inversion results obtained by different end-users. Next, to maintain the continuity of our monthly gravity field solutions, we will recover constrained SH solutions from GFO-based GPD observations in the near future.