Evaluation of hydrological and cryospheric angular momentum estimates based on GRACE, GRACE-FO and SLR data for their contributions to polar motion excitation

In geodesy, a key application of data from the Gravity Recovery and Climate Experiment (GRACE), GRACE Follow-On (GRACE-FO), and Satellite Laser Ranging (SLR) is an interpretation of changes in polar motion excitation due to variations in the Earth’s surficial fluids, especially in the continental water, snow, and ice. Such impacts are usually examined by computing hydrological and cryospheric polar motion excitation (hydrological and cryospheric angular momentum, HAM/CAM). Three types of GRACE and GRACE-FO data can be used to determine HAM/CAM, namely degree-2 order-1 spherical harmonic coefficients of geopotential, gridded terrestrial water storage anomalies computed from spherical harmonic coefficients, and terrestrial water storage anomalies obtained from mascon solutions. This study compares HAM/CAM computed from these three kinds of gravimetric data. A comparison of GRACE-based excitation series with HAM/CAM obtained from SLR is also provided. A validation of different HAM/CAM estimates is conducted here using the so-called geodetic residual time series (GAO), which describes the hydrological and cryospheric signal in the observed polar motion excitation. Our analysis of GRACE mission data indicates that the use of mascon solutions provides higher consistency between HAM/CAM and GAO than the use of other datasets, especially in the seasonal spectral band. These conclusions are confirmed by the results obtained for data from first 2 years of GRACE-FO. Overall, after 2 years from the start of GRACE-FO, the high consistency between HAM/CAM and GAO that was achieved during the best GRACE period has not yet been repeated. However, it should be remembered that with the systematic appearance of subsequent GRACE-FO observations, this quality can be expected to increase. SLR data can be used for determination of HAM/CAM to fill the one-year-long data gap between the end of GRACE and the start of the GRACE-FO mission. In addition, SLR series could be particularly useful in determination of HAM/CAM in the non-seasonal spectral band. Despite its low seasonal amplitudes, SLR-based HAM/CAM provides high phase consistency with GAO for annual and semiannual oscillation.


Introduction
One of the important tasks in contemporary geodesy is to determine and interpret changes in polar motion (PM), which is one of the parameters describing the rotation of our planet. PM is influenced by wide range of external factors with different frequency and magnitude. In previous works, particular attention was paid to assessing the role of global mass redistribution of liquid and gaseous surface components of the Earth, namely the atmosphere, ocean, continental hydrosphere and cryosphere (Brzeziński et al. 2009;Seoane et al. 2009). Such impacts were examined using information on temporal variations Open Access *Correspondence: jsliwinska@cbk.waw.pl 1 Space Research Centre, Polish Academy of Sciences, Bartycka 18A, 00-716 Warsaw, Poland Full list of author information is available at the end of the article of the gravity field provided by the Gravity Recovery and Climate Experiment (GRACE) mission (Tapley et al. 2004(Tapley et al. , 2019Wahr et al. 1998;Wouters et al. 2014), which was operational between 2002 and 2017. After the end of the GRACE activity, a successor of the GRACE mission, GRACE Follow-On (GRACE-FO), was launched in 2018 and continues its measurements to the present (Kornfeld et al. 2019;Landerer et al. 2020).
In the studies on PM variations due to global mass redistribution, particular emphasis was placed on determining the contribution of the continental hydrosphere and cryosphere (e.g., Göttl et al. 2018;Jin et al. 2010Jin et al. , 2012Meyrath and van Dam 2016;Nastula and Śliwińska 2020;Nastula et al. 2007Nastula et al. , 2011Nastula et al. , 2019Seoane et al. 2011Seoane et al. , 2012, 2020aWińska et al. 2016Wińska et al. , 2017. The sum of land hydrosphere and cryosphere components, which include soil water, surface water, groundwater, water in snow cover and ice, and water in vegetation, is commonly known as terrestrial water storage (TWS). It has been shown that variation in TWS mainly affects the annual changes in PM (Chen and Wilson 2005), which may be related to the seasonal change in rainfall and the associated variations in soil moisture. However, groundwater storage, which is a component of TWS with more long-term variability than soil moisture, may affect PM changes at interannual and longer time scales (Youm et al. 2017).
The contribution of continental hydrosphere and cryosphere to changes in PM is typically examined by analyzing two equatorial components (χ 1 and χ 2 ) of the so-called hydrological and cryospheric excitation of PM or hydrological and cryospheric angular momentum (HAM/CAM). For convenience, a complex form (χ 1 + iχ 2 ) can be used. Of the two levels of GRACE data processing, Level-2 (Bettadpur 2018a, b;Dahle et al. 2018;Yuan 2018) and Level-3 (Cooley and Landerer 2020), the former has most often been used to determine HAM/CAM (e.g., Fernández et al. 2007;Meyrath and van Dam 2016;Nastula and Śliwińska 2020;Nastula et al. 2007;Seoane et al. 2009;Śliwińska et al. 2020a, b). This method involves using ΔC 21 and ΔS 21 (degree-2 order-1) spherical harmonic (SH) coefficients of geopotential. Another method of HAM/CAM determination is the use of Level-3 GRACE data, which are represented by TWS anomalies, as presented in Nastula et al. (2019), Seoane et al. (2011), and Śliwińska et al. (2019, 2020b). TWS variations in Level-3 data can be computed either from SH coefficients of geopotential or directly from range-rate measurements by using mascon (MAS) functions.
The great importance and usefulness of the GRACE and GRACE-FO data has been emphasized by the scientific community. However, there is a gap in GRACE/ GRACE-FO estimates, spanning from the end of the GRACE mission (July 2017) to the start of GRACE-FO (May 2018). Therefore, researchers began searching for alternative data sources that could be used to fill this gap. For example, for the purpose of HAM/CAM estimation, temporal gravity field models obtained from kinematic orbits of low-Earth-orbit satellites were used by .
Another source of data for determination of HAM/ CAM is gravity field models obtained from Satellite Laser Ranging (SLR) measurements. It is commonly known that the SLR technique is critical in the precise determination of second zonal parameter (ΔC 20 ) of geopotential, which is related to the Earth's flattening (Chen and Wilson 2008;Cheng et al. 2011). Moreover, ΔC 21 and ΔS 21 coefficients, which are proportional to changes in the mass term of PM excitation, are also quite well determined by the SLR method. Another advantage of the SLR approach is that the data available cover the entire period of GRACE and GRACE-FO activity, including a one-year gap between the two missions, so there is a good potential to validate the SLR datasets through comparison with GRACE/GRACE-FO estimates. A disadvantage of the SLR approach is that the currently available SLR solutions have a relatively low resolution (maximum degree and order of 10), and thus it is not possible to obtain highresolution maps of TWS distribution. Therefore, HAM/ CAM estimates based on TWS from SLR may not be sufficiently accurate, and only ΔC 21 and ΔS 21 coefficients are commonly used for computation of precise HAM/CAM series from SLR.
The objective of this study is to estimate the contribution of the land hydrosphere and cryosphere to PM excitation as observed by GRACE. To do this, we calculate the HAM/CAM series based on different types of GRACE data, namely ΔC 21 and ΔS 21 coefficients, TWS anomalies computed from SH coefficients of geopotential, and TWS anomalies obtained from MAS data. The analysis of the near 14-year-length GRACE-derived HAM/CAM time series allows us to select the most appropriate method for determining HAM/CAM from GRACE-like mission data in various spectral bands. We also assess the quality of the determined HAM/CAM by comparing them with geodetic residual time series (GAO) data obtained from observations and geophysical models. Comparison and evaluation of HAM/CAM series is conducted separately for trends and for seasonal and non-seasonal variations. As previous studies have demonstrated the highest usefulness of MAS data in the studies of global and regional mass variations of land water and ice (e.g., Jing et al. 2020;Long et al. 2017;Luthcke et al. 2013;Rateb et al. 2017;Schrama et al. 2014), we check whether this also applies to PM excitation. In particular, we assess the improvement of HAM/CAM determined when MAS data are used instead of other GRACE data types. This study also provides preliminary HAM/CAM estimates derived from data from the first two years of the new GRACE-FO mission. As there is a gap in data between the end of the GRACE and start of the GRACE-FO mission that is nearly a year in duration, we also test the usefulness of the SLR solution in the determination of HAM/CAM. The results of this study could be used to recommend the GRACE and GRACE-FO data that are most suitable for PM excitation determination for certain oscillations.

Hydrological and cryospheric signal in observed PM excitation
Contemporary space geodesy techniques, such as VLBI (Very Long Baseline Interferometry), SLR, and GNSS (Global Navigation Satellite Systems), provide precise measurements of variations in pole coordinates (x, y). The method described in Brzezinski (1992) and Eubanks (1993) allows the determination of equatorial components (χ 1 , χ 2 ) of observed (geodetic) excitation of PM (geodetic angular momentum, GAM) from x and y. GAM series describe total variation in PM resulting from external effects, including hydrological, cryospheric, oceanic, and atmospheric signals. In order to separate the HAM/ CAM signal from GAM, the effects of atmospheric pressure and winds (described by mass and motion terms of atmospheric angular momentum, AAM), and ocean bottom pressure and currents (described by mass and motion terms of oceanic angular momentum, OAM), should be removed from GAM. The resulting series are usually denoted as geodetic residual time series and abbreviated as GAM-AAM-OAM or simply GAO. AAM and OAM, needed for computation of GAO, are calculated using data from geophysical models of atmosphere and oceans. Detailed information on determining GAO can be found in Śliwińska et al. (2020a) or Wińska and Śliwińska (2019).
GAO series represent mainly the impact of land hydrosphere and cryosphere on PM excitation; however, they also contain signals from barystatic sea-level changes, large earthquakes, post-glacial rebound, and other unmodeled signals, i.e., errors in AAM and OAM. The post-glacial rebound process mainly affects the trend in PM excitation and can be removed by applying a glacial isostatic adjustment (GIA) model Peltier et al. 2015). It should be also noted that in GAO, both mass and motion parts of HAM/CAM are present, while GRACE and GRACE-FO measurements provide only mass term of HAM/CAM. It has been stated in the work of Dill (2008) that the motion term of HAM is caused mainly by high-velocity river flows in the areas with high topographic slopes, and such regions do not necessarily correspond to the areas with high sensitivity to the mass part of HAM. The motion term of HAM shows seasonality and the maximum river flows occurs in late spring, which is typical of melting snow and ice. However, Dobslaw et al. (2010) computed the annual amplitudes of motion part of HAM, estimated from Land Surface Discharge Model (LSDM), and showed that they are approximately zero. Therefore, it can be assumed that motion term of HAM can be omitted and its absence should not affect the final results. Nevertheless, it should be kept in mind that the HAM motion term obtained from LSDM considers only the fresh water flow. The motion part of CAM is not considered in LSDM as no sophisticated model of ice mass transport and precipitation over glaciated regions is applied (Dill 2008). LSDM includes only the annual snow accumulation and melting, while the long-term ice mass is kept constant. As the integration of complex thermo-mechanical model of ice sheets is still an unsolved task, the absence of CAM motion term has not been validated yet.
In this study, for computation of GAO, we used the following datasets: There are a number of atmosphere and ocean models that can be used for AAM and OAM determination. We selected the ECMWF and MPIOM models because their use provides results that are consistent with GRACE/ GRACE-FO estimations. This is because GRACE/ GRACE-FO Atmosphere and Ocean De-Aliasing Level-1B dataset (AOD or AOD1B) (Dobslaw et al. 2017), which was used by the computing centers to remove non-tidal atmospheric and oceanic contributions from GRACE/GRACE-FO gravity fields, is based on the same two models. Therefore, the highest possible consistency between HAM/CAM and GAO in terms of AAM and OAM errors is maintained.

HAM/CAM determination from the GRACE/GRACE-FO Level-2 data
In general, two types of GRACE/GRACE-FO data can be used to calculate mass terms of HAM/CAM, namely Level-2 data and Level-3 data. Level-2 GRACE/GRACE-FO datasets have a form of SH coefficients of geopotential with a specified maximum degree and order. However, for the study of polar motion excitation, only degree-2 order-1 (ΔC 21 , ΔS 21 ) coefficients are needed. These data are also known as GSM (GRACE Satellite-only Model) coefficients. The changes in χ 1 and χ 2 equatorial components of HAM/CAM are proportional to the variations in ΔC 21 and ΔS 21 , which allows the calculation of HAM/ CAM independently from other coefficients (Chen and Wilson 2008). The relationship between χ 1 and χ 2 components of HAM/CAM and ΔC 21 and ΔS 21 coefficients of geopotential can be described as follows (Gross 2015;Eubanks 1993): where R e is the Earth's mean radius; M is the Earth's mass; A, B, and C are the principal moments of inertia for Earth (A = 8.0101 × 10 37 kg·m 2 , B = 8.0103 × 10 37 kg·m 2 , and C = 8.0365 × 10 37 kg·m 2 , from Gross 2015); A′ = (A + B)/2 is an average of the equatorial principal moments of inertia; and ΔC 21 and ΔS 21 are the normalized SH coefficients of geopotential.
Since the atmospheric and oceanic impacts are removed in all GRACE/GRACE-FO datasets using the AOD product, HAM/CAM computed in this way contains mainly the hydrological and cryospheric signals (Bettadpur 2018a;Dobslaw et al. 2017). However, similar to GAO, HAM/CAM computed from ΔC 21 and ΔS 21 contains also other signals such as the effects of post-glacial rebound, barystatic sea-level changes, or tectonic events related to large earthquakes (Śliwińska et al. 2020a). In addition, some residual atmospheric and ocean signals that remain after applying AOD data to GRACE/ GRACE-FO solutions by processing centers may not be completely removed. Such residual effects are mainly related to non-removal of the entire atmosphere and (1) ocean mass, errors in the atmospheric and ocean models used, or leakage between continents and oceans (Nastula 2017). However, the errors of the atmospheric and ocean models are also associated with GAO.

HAM/CAM determination from the GRACE Level-3 data
In contrast to the previous data type, Level-3 GRACE data have a form of terrestrial water storage (TWS) anomalies related to the average value. In the Level-3 GRACE data, a number of corrections are applied for better representation of the TWS anomalies. These corrections are: degree-2 zonal coefficient (ΔC 20 ) or Earth's flattening correction (Cheng et al. 2011), degree-1 coefficients (ΔC 10 , ΔC 11 , ΔS 11 ) or geocenter correction (Swenson et al. 2008), and glacial isostatic adjustment (GIA) correction (Peltier et al. 2018). TWS variation can be determined from gravity field observations in several ways. The most common methods involve the use of either SH coefficients of Level-2 data (see, for example, Eqs. 1 and 2 in Śliwińska et al. 2019), or MAS solutions (Rowlands et al. 2010;Watkins et al. 2015). The calculation of TWS based on SH coefficients requires summing up the effects represented by all available SH coefficients up to their maximum degree. The resulting TWS variations suffer from limitations associated with the constellation of GRACE satellites, which causes vertical stripes on the maps and requires data filtering (Swenson and Wahr 2006). However, the filtration of TWS data removes part of the actual geophysical signal, which then needs to be restored using gain factors (Landerer and Swenson 2012). Nevertheless, Göttl et al. (2019) showed that filtration has a minor effect on PM excitation obtained from TWS data. In the newest GRACE Level-3 data, TWS anomalies computed from the SH method are filtered to remove north-south stripes, and have gain factors applied to restore geophysical signals removed during filtration (Landerer and Swenson 2012).
In the MAS method, TWS anomalies are derived directly from range-rate and acceleration measurements using mascon functions (Rowlands et al. 2010;Watkins et al. 2015). This approach allows the production of TWS maps without north-south stripes and therefore does not require filtering. In this method, TWS anomalies are determined in equal-area mass concentration blocks, called mascons. In contrast to SH representation, all mascon blocks have a known geographical location.
Equatorial components of HAM/CAM can be computed from TWS anomalies (obtained from both SH and MAS solutions) using the following formulas (Barnes et al. 1983;Eubanks 1993): where C and A are Earth's principal moments of inertia; R e is Earth's mean radius; ∆q is a change in TWS; (ϕ, , t) are latitude, longitude, and time, respectively; and dS is the surface area. The factor 1.0966 accounts for the yielding of the solid Earth to rotational deformation, surface load, and core-mantle decoupling (Eubanks 1993). Here, for the determination of this relationship, we used Eqs 3-13, 3-33, 3-36 from Eubanks (1993) and the geodetic parameters of the Earth taken from Gross (2015).
The determination of TWS anomalies from raw GRACE/GRACE-FO observations requires many processing steps, including calibration and weighting of observations, removal of unknown signals using background models, or inclusion of appropriate corrections. The data processing algorithms and background models used in the newest release of GRACE solutions and first GRACE-FO datasets are similar. Moreover, for both missions, the three data types (GSM, TWS, and MAS) are available for scientific use. In this study, we consider HAM/CAM series computed using all three GRACE/ GRACE-FO data types. To do so, we used the newest solutions provided by the Jet Propulsion Laboratory (JPL): • for SH coefficients of geopotential-JPL RL06 solution; • for TWS anomalies from SH coefficients of geopotential-JPL RL06 solution; • for TWS anomalies from MAS-JPL RL06M v02 solution.
In the remainder of this paper, we use "GSM-based HAM/CAM" for HAM/CAM computed from ΔC 21 and ΔS 21 SH coefficients of geopotential, "TWS-based HAM/ CAM" for HAM/CAM computed from TWS anomalies based on SH coefficients of geopotential, and "MASbased HAM/CAM" for HAM/CAM computed from TWS anomalies based on MAS data.
At this point, it should be recalled that although nontidal atmospheric and ocean signals are mostly removed from all GRACE/GRACE-FO data products using AOD, some residual atmospheric and ocean signals are retained in the gravity fields because of the errors in models used. The uncertainties of ocean models affect the ocean areas, (3) while the errors of atmosphere models can affect both the oceans and continents. The signals that are retained over the ocean areas after applying AOD dataset are referred to in this study as residual ocean signals. The use of Eqs. 3 and 4 makes it possible to place a mask on the ocean area, so the resulting HAM/CAM series should be free of signals over oceans. In the case of TWS data, which are determined from the SH expansion, the mask was placed by computing centers (the grids for land and ocean are provided as separate datasets). For MAS-based HAM/ CAM, we may consider the so-called global series (global MAS-based HAM/CAM), in which no mask was placed on ocean areas, and the land series (land MAS-based HAM/CAM), in which the mask is applied to ocean regions in order to remove residual ocean signals. In both global and land variants, corrections listed above in this section (GIA, flattening, and geocenter) are included and the only difference between global MAS-based HAM/ CAM and land MAS-based HAM/CAM is the removal of residual ocean signals in the latter by placing the mask. By comparing these two variants of MAS-based HAM/ CAM, we can assess the effect of residual ocean signals on HAM/CAM variability and on the compatibility between HAM/CAM and GAO. In contrast, in GSM-based HAM/CAM, whose calculation is based on the relationship between ΔC 21 and ΔS 21 and χ 1 and χ 2 , the residual ocean signatures are not removed as it is not possible to use a mask for SH data. In addition, GSM-based HAM/CAM has a signal resulting from the lack of removal of GIA. However, for GAO, which is used here for evaluation of different HAM/ CAM, residual ocean signatures (due to the errors of the subtracted OAM model) and the GIA trend remain.

HAM/CAM determination from SLR
The SLR data have been extensively used in determination of Earth's oblateness, the changes of which are proportional to the variations of ΔC 20 coefficient, precisely determined with this technique (e.g., Cheng et al. 2011;Chen and Wilson 2008). Moreover, this method delivers also other low-degree coefficients, including ΔC 21 and ΔS 21 , that are necessary for PM excitation computation.
In this study, for HAM/CAM computation, in addition to GRACE/GRACE-FO datasets, we also use the latest monthly SLR solution provided by Center for Space Research (CSR) (available at http://downl oad.csr. utexa s.edu/pub/slr/degre e_5/CSR_Month ly_5x5_Gravi ty_Harmo nics.txt). The set contains coefficients up to degree and order 5 (plus ΔC 61 and ΔS 61 ). The inclusion of ΔC 61 and ΔS 61 coefficients, not present in previous releases, should improve the consistency between ΔC 21 and ΔS 21 coefficients from SLR with corresponding values determined from GRACE/GRACE-FO. The SLR series provided by CSR are estimated from analysis of five geodetic satellites, namely LAGEOS-1, LAGEOS-2, Ajisai, Starlette, and Stella (Cheng et al. 2011). From March 2012, the data from an additional satellite (LARES) were included, which contributed to the subsequent improved quality of the SLR solution. It should be also noted that background models used in SLR solution are consistent with those applied in the GRACE/GRACE-FO Release 06 datasets. The series of SH coefficients are provided both with AOD removed and retained. In this study, to maintain consistency with GRACE/GRACE-FO and GAO, we used SLR coefficients in which these effects are eliminated. The data available span from January 2002 to the present.
For computation of HAM/CAM from SLR data, we used ΔC 21 and ΔS 21 coefficients that had been recomputed into χ 1 and χ 2 according to Eqs. 1 and 2. For the resulting series, we use "SLR-based HAM/CAM" in the remainder of this paper.

Time series processing and methods of evaluation
Taking into consideration different time resolution of various datasets (GAM are daily data, AAM and OAM are 3-h data, GRACE and GRACE-FO are monthly), we downsampled the series to remove oscillations with periods shorter than one month. To do this, we used a Gaussian filter with full width at a half-maximum (FWHM) of 60 days. We conducted separate analyses for the GRACE period (between January 2003 and July 2016; we excluded periods of lower GRACE accuracy at the beginning and at the end of the mission) and for the first 2 years of GRACE-FO activity (between June 2018 and July 2020).
For the GRACE period, we considered the overall time series (without considering specific oscillations), trends, and seasonal and non-seasonal oscillations. To compute seasonal signals in HAM/CAM and GAO, we fit to the series a model consisting of a first-degree polynomial (trend) and a sum of sinusoids with periods of 365.25 days (annual oscillation), 182.625 days (semiannual oscillation), and 121.75 days (terannual oscillation). Then, to focus only on seasonal variations, we separated the trends and considered each one alone. Seasonal oscillations are presented as time series plots of χ 1 and χ 2 components of seasonal oscillation and phasor diagrams. Phasor diagrams show the amplitudes and phases of seasonal oscillation (separately for annual, semiannual, and terannual) by plotting vectors whose length corresponds to the magnitude of the amplitude, while their direction reflects the phase of the oscillation. Bizouard (2016), Lambeck (1980), and Munk and MacDonald (1960) showed that PM has an elliptic character with two circular terms: retrograde (clockwise) and prograde (counter-clockwise). Therefore, phasor diagrams are usually shown separately for prograde and retrograde oscillations (Brzeziński et al. 2009;Dobslaw et al. 2010;Nastula et al. 2011;Seoane et al. 2009Seoane et al. , 2011, 2020aWińska et al. 2016). Non-seasonal oscillations are obtained after removing trends and seasonal oscillations from the time series and they are presented as time series plots of χ 1 and χ 2 components of non-seasonal oscillation.
For the GRACE-FO period, because of the short length of data (2 years), we do not study trends and seasonal and non-seasonal oscillations, but only the overall time series.
A validation of different HAM/CAM was performed by comparison with GAO. The quality of HAM/CAM was assessed by analyzing the following parameters: standard deviation (STD) of time series, error of STD (STDerr), correlation coefficient (Corr) between HAM/CAM and GAO, percent of variance of GAO explained by HAM/ CAM (relative explained variance, Var exp ), and root mean square deviation (RMSD).
STDerr (%) values show an agreement between STD of HAM/CAM and the STD of GAO and were computed as follows: The positive values indicate higher STD for HAM/ CAM series than for GAO, whereas negative values indicate higher STD for GAO than for HAM/CAM. The most favorable value of STDerr is 0%.
Var exp shows the variability agreement between two time series and was computed using the following formula: The values of this parameter are given in % and the most favorable value is 100%.
The centered RMSD was computed as follows (Taylor 2001): where f n are reference (GAO) time series, − f is the mean value of f n , r n are evaluated time series (HAM/CAM), and − r is the mean value of r n . The most satisfactory value of RMSD is 0.
A Taylor diagram was chosen in this study to display STD of the time series, Corr, and RMSD. This method of presentation provides a statistical summary of how well time series match each other and is typically used to compare simulated data with observations (Taylor 2001). The values of STDerr and Var exp and the exact values of (5) , the parameters shown in Taylor diagrams are presented  in Tables 3, 4, 5, 6 in Appendix. Later in the paper, we analyze in detail the extent to which the MAS datasets are better than the SH data in the study of HAM/CAM. To achieve this, we analyzed the relative change of STDerr, Corr, Var exp , and RMSD. This change can be calculated as a difference between the parameter value (STDerr, Corr, Var exp , or RMSD) obtained for the MAS solution and the same parameter value obtained for the SH solution. Subsequently, for each of the four parameters, the ratio of the differences to the parameter value obtained for SH data was calculated as follows: where abs indicates the absolute value. The values of STDerr change, Corr change, Var exp change, and RMSD change are given in %. These quantities indicate by how many percent STDerr, Corr, Var exp , and RMSD increase (or decrease) when MAS data are used instead of SH data. Note that, in contrast with other formulas, in Eq. 8, we use abs in the numerator, because the sign of STDerr is not relevant. This is because the best value of STDerr is 0%, and therefore the same level of error is indicated regardless of the sign. A minus sign only indicates that HAM/CAM series has smaller STD than GAO series, whereas the absence of a minus means the opposite situation. However, both underestimation and overestimation of GAO variability by HAM/CAM is unwelcome. In Eqs. 9 and 10, the sign of Corr and Var exp is of great importance because both parameters have the maximum value (1 for Corr and 100% for Var exp ), which is also the most desirable outcome. The more negative the value of the Corr and Var exp , the less favorable it is. Therefore, abs is not used in the numerator of Eqs. 7 and 8. Moreover, abs is not considered in Eq. 11 because RMSD has only non-negative values (the best value is 0).

Analysis of overall time series and trends
The analysis was initiated with a comparison of trends in GAO, GRACE-based HAM/CAM, and SLR-based HAM/CAM (Table 1). In the case of GAO, we consider both the trend of the series in which the GIA signal is retained and the trend after removing the post-glacial rebound effect from the series. In this study, we removed the GIA contribution from GAO trend using the same model that was used in GRACE TWS and MAS data, provided by Peltier et al. (2018). Preliminary comparison of trend values indicates differences between HAM/CAM obtained from various data types (Table 1). This is not surprising, because the effects of post-glacial rebound are removed in all GRACE Level-3 datasets (TWS, MAS land, and MAS global), but are retained in GRACE Level-2 (GSM) and SLR. This naturally affects the consistency of the trend between series analyzed. It was expected that trends for SLR-based HAM/CAM are similar to those obtained for GSM-based HAM/CAM as both estimates are based on ΔC 21 and ΔS 21 coefficients in which the effects from post-glacial rebound are not eliminated.
For χ 1 , all the series agree with each other in terms of the trend sign, but there are visible discrepancies in their magnitude. In particular, TWS-based HAM/CAM and land MAS-based HAM/CAM are characterized by smaller trends than the other time series, whereas GSM-based HAM/CAM and SLR-based HAM/CAM provide high trends that are even bigger than that for GAO (both with and without GIA trend removed). Although the GIA signals are removed from global MAS-based HAM/CAM, it gives a trend as high as the GAO with the GIA retained. This may indicate that the trend in the χ 1 component of HAM/CAM is influenced not only by post-glacial rebound, but also factors such as residual ocean signals and ice signatures. This could be partly because χ 1 is more sensitive to mass variations over oceans and the Greenland ice sheet than χ 2 . Since the beginning of this century, a noticeable ice mass loss from Greenland and Antarctica has been observed by GRACE satellites (Luthcke et al. 2013;Velicogna 2009;Velicogna et al. 2014), which might affect the trend, especially in the χ 1 component of HAM/CAM. A comparison of χ 1 trends for global and land MAS-based HAM/CAM shows that using a mask that removes total signal from the ocean reduces the trend by 1.53 mas/year. Moreover, in the case of the χ 1 component of GAO, the removal of effects from post-glacial rebound reduces the GAO trend by 1.02 mas/year. The analysis of trend values for the χ 2 component indicates that there is no consistency among all series in terms of trend signs. The χ 2 trend signs for GAO (without GIA trend removed), GSM-based HAM/ CAM, and SLR-based HAM/CAM are consistent with each other (negative), because in each series the impact of post-glacial rebound is retained. However, similar to the χ 1 component, GSM-based HAM/CAM and SLRbased HAM/CAM overestimate the trend of the χ 2 component of GAO (without GIA trend removed). All other series have positive trends, because the GIA signal has been removed in all cases. There is also a consistency between trend values for TWS-based HAM/ CAM, land MAS-based HAM/CAM, and global MASbased HAM/CAM. However, all of them underestimate the GAO (with GIA trend removed) trend. Nevertheless, global MAS-based HAM/CAM provides the trend that is the most consistent with that for GAO (with GIA trend removed).
It is noticeable that for the χ 2 component, the signal from GIA has a higher impact on observed PM excitation trend than that for χ 1 , as the difference in the trend between GAO with GIA removed and GAO without GIA removed is as high as 4.60 mas/year (Table 1).
Moreover, the removal of post-glacial rebound effects causes the GAO trend sign to change from negative to positive. This is not surprising, because, on the one hand, post-glacial rebound has the greatest impact on areas of the Northern Hemisphere-Fennoscandia and the northern part of North America. On the other hand, the χ 2 component is more sensitive than the χ 1 component to mass changes over these areas (Śliwińska et al. 2020a;Youm et al. 2017). In general, removing signals from post-glacial rebound reduces the GAO trend for the χ 1 component, but increases the trend for χ 2 .
Masking ocean areas has an almost two times smaller impact on the χ 2 trend than on the χ 1 trend, which is not surprising as χ 1 is more sensitive to ocean mass changes. For both χ 1 and χ 2 components of HAM/CAM, the signals from the ocean increase HAM/CAM trends. We now compare an overall variability of GAO, GRACE-based HAM/CAM, and SLR-based HAM/ CAM by comparing their detrended time series (Fig. 1) and Taylor diagrams (Figs. 2, 3). The values of STD, STDerr, Corr, Var exp , and RMSD are also given in Table 3 in Appendix. We do not distinguish between GAO series with GIA removed and GAO series without GIA removed as post-glacial rebound mainly affects trends. The time series plots show that the χ 2 component is characterized by bigger amplitudes than χ 1 (note that there are different scales on y axis for χ 1 and χ 2 ), which has been already demonstrated in other works (e.g., Nastula et al. 2019; Seoane et al. 2011;, 2020aWińska et al. 2016Wińska et al. , 2017 and arises from continent-ocean distribution (Fig. 1). In general, the phase consistency between various time series is good for the χ 2 component, and similar minimum and maximum peaks appear for each time series. However, this consistency is lower for the χ 1 component. Figures 2 and 3 show that not removing the residual ocean signals in global MAS-based HAM/CAM increases the STD of series by 2.94 mas for χ 1 and 4.27 mas for χ 2 compared with land MAS-based HAM/CAM. Moreover, of all HAM/CAM series, global MAS-based HAM/CAM provides the highest consistency with GAO in terms of time series variability (the lowest STDerr) ( Table 3). The rest of HAM/CAM underestimate STD of GAO, especially for the χ 2 component (negative values of STDerr) (Table 3). In addition, the variability of the χ 2 component is higher than variability of χ 1 , as shown also in Fig. 1. This is as expected because higher variability and amplitudes of the χ 2 component of HAM/CAM have previously been reported (e.g., Nastula et al. 2019;Seoane et al. 2011;, 2020aWińska et al. 2016Wińska et al. , 2017. The comparison of GSM-based HAM/CAM and SLR-based HAM/CAM indicates that both series present similar STD consistency with GAO for χ 2 ; however, series from GRACE GSM are more compatible with GAO for χ 1 (Figs. 2, 3, Table 3).
The highest consistency between GAO and HAM/ CAM as measured by Corr, Var exp , and RMSD is observed for global MAS-based HAM/CAM (Corr of 0.72 and 0.84 for χ 1 and χ 2 , respectively; Var exp of 36% and 69% for χ 1 and χ 2 , respectively; and RMSD of 6.39 and 7.66 mas for χ 1 and χ 2 , respectively) (Figs. 2, 3, Table 3). GSM-based HAM/CAM, however, in which residual ocean signatures also remain, is characterized by the lowest agreement with GAO of all HAM/CAM series (Corr of 0.52 and 0.60 for χ 1 and χ 2 , respectively; Var exp of 7% and 36% for χ 1 and χ 2 , respectively; and RMSD of 7.71 and 10.94 mas for χ 1 and χ 2 , respectively). This might indicate that incomplete removal of the effects of ocean bottom pressure and currents from GAO and GSM is not a reason for these inconsistencies, but the SH approach might be less appropriate for estimating mass-related PM excitation than the MAS approach.
Of the two HAM/CAM series in which the mask is put on ocean areas, namely TWS-based HAM/CAM and land MAS-based HAM/CAM, the former provides better results for χ 1 , while the latter ensures a higher consistency with GAO for χ 2 . Overall, a higher consistency between GAO and HAM/CAM is observed for χ 2 than for χ 1 , and this has already been reported in previous works (Meyrath and van Dam 2016;Seoane et al. 2011;, 2020a and result from the spatial distribution of the main continents and oceans. Interestingly, in the HAM/CAM series obtained from ΔC 21 and ΔS 21 coefficients provided by the SLR technique are distinguished by higher consistency with GAO than corresponding series computed using ΔC 21 and ΔS 21 coefficients delivered by GRACE (GSM-based HAM/ CAM). However, previous studies have shown that SLR data are less appropriate for determination of HAM/ CAM than GRACE data and are mainly useful for determination of Earth's oblateness (Chen and Wilson 2008). In previous studies concerning PM excitation, SLR data were usually used in combination with GRACE or other datasets . Nevertheless, analysis of the newly released SLR series provided by CSR shows a clear improvement in SLR-based HAM/CAM.

Analysis of seasonal and non-seasonal variations
It is well known that seasonal signals are one of the strongest in PM excitation (Brzeziński et al. 2009;Dobslaw et al. 2010;Nastula et al. 2019;;  Wińska et al. 2017), therefore we now divide the analysis into seasonal and non-seasonal parts. We begin with the analysis of seasonal variations, which are presented as time series plots of χ 1 and χ 2 components (Fig. 4). As with the overall time series analysis, a validation of seasonal variations in HAM/CAM is based on Taylor diagrams (Figs. 5, 6). The corresponding values of STD, STDerr, Corr, Var exp , and RMSD are also given in Table 4 in Appendix.
An overall comparison of the seasonal time series indicates that the series agree quite well with each other, but some differences in amplitudes occur (Fig. 4). In particular, of all GRACE-based estimates, GSM-based HAM/ CAM have the smallest amplitudes while MAS-based  HAM/CAM have the biggest amplitudes of seasonal variation. This indicates that the SH approach may reduce amplitudes of seasonal oscillation. TWS-based HAM/ CAM, however, which is also based on SH coefficients, provides slightly bigger seasonal amplitudes than amplitudes of GSM-based HAM/CAM. These findings are supported by STD values of seasonal series (Figs. 5, 6, Table 4); GSM-based HAM/CAM are characterized by the lowest variability, while MAS-based HAM/CAM series have the highest STD values (global MAS-based HAM/CAM for χ 1 and land MAS-based HAM/CAM for χ 2 ). Notably, land MAS-based HAM/CAM series are characterized by small STDerr (the second lowest value for χ 1 and the lowest value for χ 2 ), which indicates a high agreement with GAO in terms of time series variability. Analysis of SLR-based HAM/CAM series shows that SLR data produce series with the smallest amplitudes for χ 1 . For χ 2 the amplitudes of SLR-based HAM/CAM are similar to those obtained from GRACE GSM. This is primarily because both SLR-and GSM-based HAM/CAM series are determined using ΔC 21 and ΔS 21 coefficients of geopotential. The use of SH coefficients instead of mascon approach can be a reason of lower amplitudes in HAM/ CAM. Notably, for the seasonal χ 1 component, SLRbased HAM/CAM exhibits an additional peak that precedes the annual maximum and has an amplitude about two times smaller than the amplitude of annual variation. In general, it can be stated that for seasonal oscillations, the choice of data type has a higher impact on amplitude consistency between HAM/CAM and GAO than on phase agreement.
Taylor diagrams (Figs. 5, 6) and parameters shown in Table 4 indicate that the MAS approach better determines seasonal oscillations in HAM/CAM than the SH method. In particular, global MAS-based HAM/CAM is best for determination of the χ 1 component (Corr of 0.91; Var exp of 69%; and RMSD of 2.25 mas), whereas land MAS-based HAM/CAM provides the most satisfactory results for χ 2 (Corr of 0.98; Var exp of 93%; and RMSD of 1.90 mas). This is unsurprising because χ 2 is more sensitive to mass changes over land, whereas χ 1 is more responsive to variations over oceans. In GAO and global MAS-based HAM/CAM some residual signals from ocean errors and leakage are retained, which results in a higher agreement of these series in the case of the χ 1 component. In general, to obtain the highest possible consistency between seasonal variation in HAM/ CAM and seasonal variation in GAO, one should use global MAS data (no mask) for determination of χ 1 and land MAS data (with mask) for χ 2 . However, to receive a "pure" hydrological plus cryospheric signal in PM excitation, without any residual ocean signatures, only the data from continents should be included. It should also be noted that, for seasonal HAM/CAM oscillations, SLRbased HAM/CAM series are less consistent with GAO than HAM/CAM obtained from any type of GRACE data.
We now separate seasonal oscillations into annual and semiannual variations and present their phasor diagrams separately for prograde and retrograde terms (Fig. 7). We do not consider terannual oscillations as Śliwińska et al. (2020a) showed that the amplitudes and phases of these oscillations are smaller than the errors in their determination. For annual oscillations, among all GRACEderived series, land MAS-based HAM/CAM has the best phase agreement with GAO, especially in the retrograde term. The annual amplitude of GAO is also well captured by these series. The GRACE series based on the SH method (TWS-based HAM/CAM and GSM-based HAM/CAM) visibly underestimates annual amplitudes of GAO for both prograde and retrograde terms. This indicates that methods based on SH cannot reconstruct the full observed amplitude of the annual oscillation in PM excitation. For semiannual variation, the amplitudes of all HAM/CAM series are smaller than amplitudes of GAO. Among all GRACE estimates, global MAS-based HAM/CAM appears to be the most consistent with GAO in terms of phases. However, it should be kept in mind that the annual variation has definitely higher contribution to seasonal PM excitation change than semiannual oscillation. In general, the highest consistency between seasonal HAM/CAM and GAO is obtained for annual retrograde oscillation. The phasor diagrams also show that not removing residual ocean signals in global MASbased HAM does not have a noticeable impact on annual phases. However, such signals increase annual prograde amplitudes and reduce annual retrograde amplitudes. Notably, SLR-based HAM/CAM, although it definitely underestimates the amplitudes of seasonal oscillations, provides high consistency with GAO phases. This substantial correspondence is particularly evident for prograde oscillation, for which the SLR results for phases are even better than those obtained for the GRACE MAS data.
We now analyze time series of non-seasonal variations in HAM/CAM and GAO (Fig. 8), together with their Taylor diagrams (Figs. 9, 10). The corresponding values of STD, STDerr, Corr, Var exp , and RMSD are also given in Table 5 in Appendix. Both TWS-based HAM/ CAM and land MAS-based HAM/CAM have the smallest amplitudes, which is also supported by the STD values (Figs. 8-10, Table 5). In contrast, the highest STD was found for global MAS-based HAM/CAM, which was also the most consistent with STD of GAO (the lowest STDerr). This may indicate that residual ocean signatures have an impact on non-seasonal changes in HAM/CAM especially. It is also noticeable that for non-seasonal oscillations, HAM/CAM series from SLR and GRACE GSM have good phase and amplitude agreement, which results from the same method of computing HAM/CAM. The SLR-based and GSM-based HAM/CAM series have also similar STD (Figs. 9, 10).  Table 5) indicate that the highest consistency between non-seasonal variations in GAO and non-seasonal variations in HAM/CAM can be obtained when global MAS data are used (Corr of 0.64 and 0.83 for χ 1 and χ 2 , respectively; Var exp of 25% and 64% for χ 1 and χ 2 , respectively; and RMSD of 5.98 and 6.90 mas for χ 1 and χ 2 , respectively). In contrast, GSM-based HAM/CAM, despite having residual ocean signals that have not been removed, provides the lowest level of agreement with GAO. It is also noticeable that TWS-based HAM/CAM and land MAS-based HAM/CAM give almost the same level of consistency with GAO, which mainly result from masking oceans in both datasets. Taylor diagrams and parameters listed in Table 5 also indicate that for the non-seasonal spectral band, SLR observations provide better consistency of HAM/CAM with GAO than GRACE GSM data. This indicates that SLR observations are especially appropriate for determining non-seasonal variations in PM excitation. However, these results are not as satisfactory as for GRACE global MAS data, which, as noticed earlier in this study, are the most appropriate for the study of PM excitation due to land hydrosphere and cryosphere in the non-seasonal spectral band.

Comparison between MAS and SH solutions
The analyses presented in the previous sections showed that, in general, MAS datasets provide higher consistency between HAM/CAM and GAO than solutions based on SH coefficients. We now examine the magnitude of change in consistency between HAM/CAM and GAO when MAS data are used instead of SH solutions, though analysis of the relative change of STDerr, Corr, Var exp , and RMSD. The changes in STDerr, Corr, Var exp , and RMSD are computed here for two cases: first, for land data (land MAS-based HAM/CAM compared to TWS-based HAM/CAM) and second, for global data (global MASbased HAM/CAM compared to GSM-based HAM/ CAM). We consider overall time series (without separating specific oscillations), seasonal oscillations, and nonseasonal changes. Table 2 presents STDerr, Corr, Var exp , and RMSD changes computed according to Eqs. 8-11. Table 2 shows that for land data, the use of the MAS solution improves the agreement between overall HAM/ CAM and GAO series mostly for the χ 2 component (improvement in STDerr by 29%, Corr by 5%, Var exp by 18%, and RMSD by 9%), while for the χ 1 component only STDerr is improved (by 30%). The highest change in consistency between HAM/CAM and GAO is observed for seasonal variation in the χ 2 component. The improvement in STDerr by 65%, Var exp by 43%, and RMSD by 53%, without any change in Corr, indicates that the use of land MAS data especially improves the amplitudes of seasonal oscillation. Indeed, Figs. 3 and 4 show that land MAS-based HAM/CAM has amplitudes are more consistent with GAO than TWS-based HAM/CAM. For non-seasonal variation, there is no noticeable change in results when land MAS-based HAM/CAM is considered instead of TWS-based HAM/CAM.
Analysis of results obtained for global data shows that the increase in overall consistency between HAM/CAM and GAO is visible for both equatorial components, with values ranging between 17 and 414%, and with the highest value for Var exp in the χ 1 component. Both equatorial components are also improved for non-seasonal oscillations (the values of improvement range between 19 and 267%, with the highest one for Var exp in χ 1 component). Similar to land data, the improvement in agreement between seasonal variation in HAM/CAM and seasonal variation in GAO is more noticeable for χ 2 . This is mainly due to underestimation of seasonal amplitudes of GAO by GSM-based HAM/CAM.

Preliminary analysis of HAM/CAM series from GRACE-FO
This section presents preliminary estimates of GSMbased, TWS-based, and global and land MAS-based HAM/CAM from the new GRACE-FO mission. A comparison with data from GAO and SLR is also provided. We analyze overall detrended time series from first 2 years of the mission activity (June 2018-July 2020). A comparison plot between GAO and various HAM/CAM series is presented in Fig. 11; Taylor diagrams are given in Figs. 12,13; and the values of STD, STDerr, Corr, Var exp , and RMSD are provided in Table 6 in Appendix. Figure 11 shows that for χ 1 , land MAS-based HAM/ CAM and TWS-based HAM/CAM series have quite similar amplitudes to that of GAO, but there is some phase shift between series. GSM-based HAM/CAM, global MAS-based HAM/CAM and SLR-based HAM/ CAM are in opposite phases to GAO. They also have visibly higher amplitudes than the other time series. This is probably because of the non-removal of residual ocean signals in GSM-based HAM/CAM, global MAS-based HAM/CAM and SLR-based HAM/CAM. In terms of χ 2 , the phase consistency between GAO and various HAM/CAM is higher than that for χ 1 , but discrepancies in amplitudes still exist. Similar to χ 1 , data from land only (in land MAS-based HAM/CAM and TWS-based HAM/CAM) enables us to obtain lower amplitudes than amplitudes of HAM/CAM based on other data types. Unsurprisingly, the consistency between various HAM/ CAM estimates is higher for χ 2 . In particular, each series shows a maximum in mid-2019, which is also apparent for GAO. Both land data (land MAS-based HAM/CAM compared to TWS-based HAM/CAM) and global data (global MAS-based HAM/CAM compared to GSM-based HAM/CAM) are considered in this table. We examine overall detrended time series, seasonal variations, and non-seasonal variations. The values showing an improvement in results after using MAS data are italicized. Note that for STDerr change and RMSD change the improvement is when variation is negative (STDerr and RMSD decreased after using MAS data), whereas for Corr change and Var exp change, the improvement occurs when the parameter value is positive (Corr and Var exp increased after using MAS data)

STDerr change (%)
Corr change (%) Var exp change (%) RMSD change (%) Overall ( Table 6 confirm that MAS data from land provide better consistency between HAM/CAM and GAO than global data. However, when analyzing the GRACE solutions, we found that for overall time series of GRACE-based HAM/CAM, global MAS datasets give slightly better results than land data. Nevertheless, it should be kept in mind that the length of currently available GRACE-FO series (24 months at the time of conducting this study) may affect objective assessment of HAM/CAM. TWS-based HAM/CAM presents a similar correspondence with GAO as land MAS-based HAM/ CAM, which was also observed for GRACE.
All in all, the actual accuracy of HAM/CAM computed using data from new GRACE-FO mission meets expectations (Corr values vary from − 0.12 to 0.54, Var exp values vary from − 364 to 28%, and RMSD values vary from 5.34 to 9.83 mas, depending on data type and equatorial component considered). The high consistency between HAM/CAM and GAO observed in the period of best GRACE performance has not yet been reached. However, the regular appearance of data from the following months of GRACE-FO operation will improve this consistency.
A detailed analysis of series obtained from SLR data indicates that for the period of GRACE-FO activity, SLRbased HAM/CAM series are less consistent with GAO than HAM/CAM obtained from any type of GRACE-FO data. This is partly because GRACE-FO has improved hardware, such as additional laser ranging interferometer for better measurement of distance between satellites or an additional star camera to improve satellite orientation (Kornfeld et al. 2019). The above upgrades undoubtedly contributed to the ameliorated quality of GRACE-FObased HAM/CAM.   Śliwińska et al. Earth, Planets and Space (2021) 73:71 Conclusions A comparison of HAM/CAM computed from different GRACE data products indicated that the use of MAS solutions provided the highest consistency between HAM/CAM and GAO for all considered spectral bands. However, depending on the oscillations and equatorial components considered, global or land MAS data were more appropriate for this purpose. The advantage of MAS data over more commonly used SH solutions is particularly evident in the seasonal spectral band (global MAS-based HAM/CAM provided the highest consistency with GAO for χ 1 , and land MAS-based HAM/ CAM provided the most satisfactory results for χ 2 ). For non-seasonal oscillations, MAS solutions (global for both χ 1 and χ 2 ) also performed better than the SH series. A similar situation was found for trends in HAM/CAM, which are most consistent with GAO trends when global MAS (for χ 1 ) or land MAS (for χ 2 ) are exploited. Analysis of overall HAM/CAM series from the first 2 years of the GRACE-FO mission data confirmed that the use of land MAS data enabled the highest consistency between HAM/CAM and GAO. The actual RMSD of HAM/ CAM based on GRACE-FO land MAS data is as high as 6.03 mas for χ 1 and 5.34 mas for χ 2 . With an increasing data length, GRACE-FO-based HAM/CAM quality will increase. We also noted that the newest release of the SLR solution can be used for determination of HAM/CAM to fill the one-year-long data gap between the end of GRACE activity and the start of the GRACE-FO measurements. In addition, the SLR series could be particularly useful in determination of HAM/CAM in the non-seasonal spectral band, as the performance of this technique is the best for these oscillations. Despite its low seasonal amplitudes, SLR-based HAM/CAM provided high phase consistency with GAO in this spectral band. The key advantage of SLR data is the length of the time series and the regularity of the subsequent monthly observations.
The main conclusion of this study is that the use of GRACE/GRACE-FO MAS solutions in determining PM excitation due to the land hydrosphere and cryosphere allows us to obtain higher consistency with observed excitation than exploiting data based on SH coefficients. We therefore recommend the use of the MAS series in studies concerned with analyses of variations in PM. The MAS method is more recent than the SH method, and has several advantages. First of all, it provides a better separation of land and ocean signals, which is particularly important for areas along the coastlines. Moreover, TWS maps obtained from MAS solutions, in contrast to TWS maps based on SH data, do not need to be filtered. This is because a priori constraints in MAS solutions help to filter out noise from the GRACE/GRACE-FO observations at the Level-2 processing step, which is a more rigorous approach than the empirical post-processing filtering applied to TWS grids derived from SH coefficients. It is known that filtering removes some part of real geophysical signal, which has to be restored using gain factors. However, as gain factors are usually determined using hydrological models (e.g., Landerer and Swenson 2012;Zhang et al. 2016), which are not free of errors, this operation could introduce additional uncertainties. The problem does not occur for the MAS data. Moreover, in MAS solutions, TWS distribution is obtained directly from GRACE/GRACE-FO observations using mascon functions, while in SH data, first the coefficients of geopotential are determined, and then TWS maps are developed on their basis. Additional processing step can increase the errors of resulting SH-based TWS.

STD (mas)
STDerr ( Table 6 Standard deviation (STD) of overall detrended time series of GAO, GRACE-FO-based HAM/CAM (GSM, TWS, MAS land, and MAS global), and SLR-based HAM/CAM series, error of STD (STDerr), correlation coefficients (Corr) between HAM/CAM and GAO, percentage of variance in GAO explained by HAM/CAM (Var exp ), and root mean square error (RMSD) of HAM/CAM Note that the critical value of correlation coefficient for 10 independent points and a 95% confidence level is 0.55, while the standard error of the difference between two correlation coefficients is 0.53 STD (mas) STDerr (%) Corr Var exp (%) RMSD (mas)