Validation of GRACE and GRACE-FO Mascon Data for the Study of Polar Motion Excitation

In this study, we calculate the hydrological plus cryospheric excitation of polar motion (hydrological plus cryospheric angular momentum, HAM/CAM) using mascon solutions based on observations from the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GRACE-FO) missions. We compare and evaluate HAM/CAM computed from GRACE and GRACEFO mascon data provided by the Jet Propulsion Laboratory (JPL), the Center for Space Research (CSR), and the Goddard Space Flight Center (GSFC). A comparison with HAM obtained from the Land Surface Discharge Model is also provided. An analysis of HAM/CAM and HAM is performed for overall variability, trends, and seasonal and non-seasonal variations. The HAM/CAM and HAM estimates are validated using the geodetic residual time series (GAO), which is an estimation of the hydrological plus cryospheric signal in geodetically observed polar motion excitation. In general, all mascon datasets are found to be equally suitable for the determination of overall, seasonal, and non-seasonal HAM/CAM oscillations, but some differences in trends remain. The use of an ellipsoidal correction, implemented in the newest solution from CSR, does not noticeably affect the consistency between HAM/CAM and GAO. Analysis of the data from the first two years of the GRACE-FO mission indicates that the current accuracy of HAM/CAM from GRACE-FO mascon data meets expectations, and the root mean square deviation of HAM/CAM components are between 5 and 6 milliarcseconds. The findings from this study can be helpful in assessing the role of satellite gravimetry in polar motion studies and may contribute towards future improvements to GRACE-FO


Introduction
The Gravity Recovery and Climate Experiment (GRACE) mission, jointly managed by the American National Aeronautics and Space Administration (NASA) and the German Aerospace Center (Deutsches Zentrum für Luft-und Raumfahrt; DLR), was operational between 2002 and 2017 [1][2][3][4]. GRACE was launched to provide continuous and highquality observations of the temporal gravity field variations caused by mass changes of the atmosphere, oceans, hydrosphere, and cryosphere around the world. The mission concept was based on a constellation of twin satellites placed in the same low Earth orbit (about 500 km above the Earth surface), spaced approximately 220 km from each other. The payload of the GRACE mission allowed precise measurement of the variations of distance between two satellites, which reflected the Earth's gravity field changes. The unquestionable success of the mission, which terminated in 2017, contributed to the development of another mission based on the same concept, GRACE Follow-On (GRACE-FO) [5,6], which began its measurements in May 2018 and continues to the present. The currently operating GRACE-FO mission, in addition to the microwave K-band ranging system (KBR) previously used on GRACE satellites, is testing the laser ranging interferometer (LRI), a new instrument for inter-satellite ranging. of land water, snow, and ice as the area of the ocean is usually masked. In addition, signals from post-glacial rebound are always modelled and removed. Comparisons of HAM/CAM determined from Level-2 and Level-3 GRACE data have been provided in previous research (e.g., [42,47]).
Compared with TWS based on SH, TWS derived from mascon datasets has several advantages, such as better separation of land and ocean signals or no filtering required, and thus should be more appropriate for HAM/CAM determination. However, although many studies have used GRACE mascon datasets to study water and ice mass changes around the world (e.g., [11,[51][52][53][54][55]), to date, there were only few studies that have used mascon solution to analyze the PM changes due to the land hydrosphere and cryosphere. The example work is [47], where HAM/CAM obtained from mascon-based TWS was compared to HAM/CAM computed from SH-based TWS, and HAM/CAM calculated using ∆C 21 , ∆S 21 coefficients of geopotential. It was concluded that the use of mascon data provides the highest consistency between HAM/CAM and geodetic observations of PM excitation. However, only one mascon solution was tested in the mentioned study, while several solutions of this type are available. Therefore, in this paper, we determine HAM/CAM series from several available GRACE/GRACE-FO mascon solutions and analyze their variability. We also assess the quality of determined HAM/CAM and HAM by comparing them with the hydrological plus cryospheric signal in observed PM excitation (GAO), obtained from geodetic measurements and geophysical models. This could help decide which of the available mascon solutions is the most suitable for the study of PM excitation. Mascon-based HAM/CAM are also contrasted with HAM derived from the Land Surface Discharge Model (LSDM), a hydrological model that has demonstrated the highest consistency with GAO among the models tested to date [42,46]. The analyses are performed for the overall HAM/CAM series without considering specific oscillations, but also after decomposing them into trends and seasonal and non-seasonal variations. The mascon-based HAM/CAM series from the GRACE-FO mission are also considered, which may give a preliminary assessment of the accuracy of the new mission.

GRACE and GRACE-FO Mascon Data
In this paper, we use the newest mascon data provided by JPL, CSR, and the NASA Goddard Space Flight Center (GSFC). Recently, the CSR team has released the new CSR RL06M v02 solution which, in contrast to other solutions, defines mascons on an ellipsoid. Therefore, for CSR, we consider both the newest (CSR RL06M v02) and the previous solution (CSR RL06M v01). Basic information about mascon solutions used in this study, together with corrections included, data sources, and references, are given in Table 1.
In the JPL RL06M v02 solution (abbreviated here as JPL), TWS anomalies were computed for each equal-area 3 • × 3 • spherical cap mascon block [8,56]. In the final dataset, the TWS variations are sampled to 0.5 • × 0.5 • longitude-latitude grids; however, it should be noted that the real resolution of JPL RL06M v02 is 3 • and thus land-grid-scaling factors provided by JPL have to be used to undo the data smoothing resulting from the use of 3 • mascons as the basis. Compared with the previous version (JPL RL05M), JPL RL06M v02 data have improved separation of land and ocean mascons, provided through the use of a Coastline Resolution Improvement (CRI) filter, which reduces leakage errors across coastlines [56]. The processing teams at JPL offer mascon solutions with and without the CRI filter applied. In the current study, we use JPL mascon data with the CRI filter applied.
In CSR RL06M v01 (abbreviated here as CSR v01), TWS anomalies are given in 0.25 • × 0.25 • longitude-latitude grids, but they represent the equal-area grids of size 1 • [57]. However, it should be remembered that the neighboring grids are designed to be correlated with each other. In contrast to the previous CSR RL05M release, CSR RL06M solution uses a newly defined grid in which the hexagonal tiles that span the coastline are split into two (ocean and land tiles) to minimize the leakage between land and ocean signals. Because the CSR data are based on smaller mascons than those in the JPL solution and separate mascons into ocean and land along the coastlines, there is no need to apply scaling factors in this case.
The recently released CSR RL06M v02 solution (abbreviated here as CSR v02) has been processed similarly to CSR RL06M v01. There are only two differences between these releases. First, in contrast to the CSR RL06 v01 solution, CSR RL06 v02 mascon grids have been corrected for representation on ellipsoidal Earth according to methodology shown in [58]. This correction was applied separately for land and ocean areas to prevent any leakage of the land signal into the ocean [59]. Second, while computing GRACE-FO mascons in the new CSR solution, the ∆C 30 coefficient was replaced with a more accurate estimate from SLR. All other processing steps remain unchanged. GIA correction ICE6G-D model [66] ICE6G-D model [66] ICE6G-D model [66] ICE6G model [67]  In the GSFC v02.4 solution (abbreviated here as GSFC), TWS anomalies were computed in each 1 • × 1 • equal-area square mascon [53]. In contrast to the mascon solutions from JPL and CSR, TWS variations are not sampled to the regular longitude-latitude grids. Another difference is that the regional constraints were applied using seven different regions (Greenland Ice Sheet with elevation above 2000 m, Greenland Ice Sheet with elevation below 2000 m, Antarctic Ice Sheet with elevation above 2000 m, Antarctic Ice Sheet with elevation below 2000 m, Gulf of Alaska glacier area, land, and ocean) to define each mascon according to its geographical category. As a result, mascons defined for a specific region do not overlap with any other area. Metadata for the defined basins are provided in the mascon data files and allow the user to select a region of interest. There is no need to apply scaling factors to the data.
The goal of this study is to determine the equatorial components (χ 1 , χ 2 ) of HAM/CAM from the distribution of TWS provided by GRACE/GRACE-FO mascon solutions. To do so, we use a well-known formula [30,69]: where C and A are Earth's principal moments of inertia with C − A = 2.6398·10 35 kg m 2 and C = 8.0103·10 37 kg m 2 [70]; R 2 e is the square of the Earth's mean radius; ∆q is a change in TWS; ϕ, λ, and t are latitude, longitude, and time, respectively; and dS is the surface area. The factor 1.0966 reflects the yielding of the solid Earth to rotational deformation, surface load, and core-mantle decoupling [69].
It should be kept in mind that non-tidal atmospheric and ocean signals are mostly removed from all GRACE data products by using the AOD dataset [29]. However, there is always a residual signal above the ocean areas, which may be related to errors in ocean model used to determine AOD or leakage between land and ocean [71]. One method of separating the land signal from the residual ocean signal is to place a mask above the ocean area, and this is done in our study.

Land Surface Discharge Model (LSDM)
In this study, we also consider HAM series obtained from the LSDM provided by GFZ (ftp://esmdata.gfz-potsdam.de/EAM, accessed on 1 July 2020). LSDM is a hydrological model, elaborated in a regular 0.5 • × 0.5 • global latitude-longitude grid and with a daily time resolution [72]. LSDM provides global variations of water mass storage and water mass fluxes, and the HAM time series. The LSDM was processed consistently with the ECMWF (European Center for Medium-Range Weather Forecasts) [73] and MPIOM (Max Planck Institute Ocean Model) [74] models, which were also used in the GRACE/GRACE-FO AOD data product, therefore the mass balance is maintained when analyzing the results of LSDM, ECMWF and MPIOM. No sophisticated ice model is applied in LSDM, so it may not be possible to determine a reliable cryospheric angular momentum (CAM) series. However, LSDM includes the annual snow accumulation and melting in glaciated regions, which guarantees that the main part of the seasonal cycle over ice-sheets is captured. Nevertheless, the long-term ice mass is kept constant [72]. In contrast, total mass changes over polar ice sheets and mountain glaciers are measured by GRACE/GRACE-FO. Therefore, further in the paper, we use "HAM/CAM" for excitation series determined from GRACE/GRACE-FO solutions, and "HAM" from series taken from the LSDM. The LSDM also contains shallow groundwater layer that is not present in many other hydrological models.
In order to make a comparison of HAM from the hydrological model with HAM/CAM from GRACE and with GAO more reliable, the sea-level angular momentum (SLAM) should be considered [75]. SLAM reflects barystatic sea level changes resulting from the inflow of water from land into the oceans that may affect HAM. Considering SLAM is important to maintain the mass balance between the atmosphere, ocean, and land hydrosphere. The scientists from GFZ provide atmospheric angular momentum (AAM), oceanic angular momentum (OAM), HAM, and SLAM series, which were processed consistently in order to keep the global mass balance. Maintaining the mass balance means that the total mass in the atmosphere, oceans, and continental hydrosphere is constant at any given time. It was recommended by GFZ to add SLAM to HAM in order to obtain higher agreement with geodetic observations of PM. However, it has been found that considering the mass balance is of particular importance in the case of seasonal variations in length-of-day, while the impact on PM is smaller [75]. Nevertheless, in this study, we consider only the sum of HAM and SLAM (HAM+SLAM). The time series of χ 1 and χ 2 components of both HAM and SLAM, which are consistent with each other, are provided by GFZ (available at ftp://esmdata.gfz-potsdam.de/EAM, accessed on 1 July 2020) and used in this study.

Hydrological Plus Cryospheric Signal in Observed PM Excitation
The HAM/CAM series, determined from monthly GRACE/GRACE-FO gravity field estimates, and HAM series obtained from hydrological model LSDM, are evaluated here by comparison with the hydrological plus cryospheric signal in observed (or geodetic) PM excitation. The geodetic excitation, expressed as geodetic angular momentum (GAM), can be computed from measured variations in pole coordinates (x, y). The precise observations of (x, y) are provided by space geodesy techniques, such as SLR, VLBI (very long baseline interferometry), and GNSS (global navigation satellite systems). Pole coordinates, together with other Earth orientation parameters (EOP), namely length-of day and precession and nutation, are made available by the International Earth Rotation and Reference Systems Service (IERS) (https://www.iers.org/IERS/EN/Home/home_node.html, accessed on 1 March 2020). The calculation of (χ 1 , χ 2 ) equatorial components of GAM from (x, y) pole coordinates can be performed using the algorithms shown in [69,76]. The resulting GAM series reflect PM variations due to all external effects, including those from land hydrosphere, cryosphere, atmosphere, and oceans.
In order to separate the hydrological and cryospheric signals from GAM, as needed for assessment of HAM/CAM and HAM, one should remove the impact of atmospheric pressure and winds (described by AAM), and ocean bottom pressure and currents (OAM). Usually, AAM and OAM are determined from geophysical models of the atmosphere and ocean. The resulting series are described as geodetic residuals, GAM−AAM−OAM or simply GAO. The choice of AAM and OAM models should not affect phases of GAO, but may have some influence on the amplitudes [42]. Details on GAO computation can be found in previous research (e.g., [42,46]).
The GAO series reflects mainly the impact of the land hydrosphere and cryosphere on PM excitation; however, similar to HAM/CAM from GRACE/GRACE-FO, they also contain effects from SLAM and tectonic events. In addition, signals from post-glacial rebound, which are modeled and removed in GRACE/GRACE-FO Level-3 data, remain in GAO. Nevertheless, the effects from post-glacial rebound mainly affect the trend in GAO and can be eliminated using a GIA model that applies a readjustment in response to the last deglaciation and current changes in ice mass [67,68].
To determine GAO, we use the time series of (x, y) coordinates of the pole, obtained from EOP 14 C04 combined solution [77], provided by the IERS (https://www.iers.org/ IERS/EN/DataProducts/EarthOrientationData/eop.html, accessed on 1 March 2020). The (x, y) coordinates are used for computation of (χ 1 , χ 2 ) components of GAM, according to the algorithm shown at http://hpiers.obspm.fr/eop-pc/index.php?index=excitactive& lang=en (data accessed on 1 March 2020). For AAM and OAM, we used ESMGFZ (Earth System Model GFZ) datasets, made available by GFZ (http://rz-vm115.gfz-potsdam.de: 8080/repository, accessed on 1 March 2020). The ESMGFZ data provide time series of the χ 1 and χ 2 components of AAM and OAM (both mass and motion terms). AAM is based on the ECMWF model [73], while OAM is based on MPIOM [74]. These models were selected because they were used by the GRACE/GRACE-FO team to process the AOD data product, which is essential to remove atmospheric and oceanic non-tidal contributions from GRACE/GRACE-FO gravity fields. HAM from LSDM is also consistent with AAM from ECMWF and OAM from MPIOM (the mass balance is maintained).

Time Series Processing
It should be noted that GAM and HAM based on LSDM model are given at daily intervals, AAM and OAM are given at 3-hour intervals, while GRACE/GRACE-FO datasets are monthly. Therefore, some short-period signals present in GAM, OAM, AAM, and HAM are not included in GRACE/GRACE-FO estimates of HAM/CAM, which certainly affects the objective comparison of the time series. For this reason, we apply a Gaussian filter to the series to remove oscillations with periods shorter than one month. Our analyses are divided into those performed for the period of GRACE mission activity (between January 2003 and July 2016, because data for GSFC is not available before and after this period, Section 3.1) and analyses for the first two years of GRACE-FO operation (between June 2018 and July 2020, Section 3.2).
For the GRACE period, we analyze the following oscillations: Trends, overall detrended time series (without separating specific oscillations, but with trends removed), seasonal variations, and non-seasonal variations. In previous works, the seasonal oscillations in PM excitation have been of particular interest [27,28,73], because they are one of the strongest in PM. Trends and seasonal oscillations are computed together by fitting a seasonal model to the series using the least-squares method. This model consists of a first-degree polynomial (trend) and a sum of three sinusoids with the periods of 1 year (annual oscillation), 1/2 year (semiannual oscillation), and 1/3 year (terannual oscillation). The trend and seasonal fluctuations are then separated from each other and analyzed independently. Non-seasonal oscillations are obtained here after removing trends and seasonal oscillations from the overall time series. For the period of GRACE-FO activity, we consider only the overall detrended time series because of the short length of data (2 years).
To examine the amplitude distribution for different spectral bands, we use amplitude spectra of the complex-valued components of χ 1 + iχ 2 , which is a common method used in PM studies [41,46,56,78].
The motion of the Earth's pole is two-dimensional and is usually described by its coordinates (x, y) given in the terrestrial reference frame, while PM excitation is described by (χ 1 , χ 2 ) components. Both (x, y) and (χ 1 , χ 2 ) are expressed in a Cartesian plane coordinate system. However, it is known that the motion of the pole is a combination of two circular movements, namely prograde, which is a circular motion in a counter-clockwise direction, and retrograde, which is a circular motion in a clockwise direction [79][80][81]. The sum of prograde and retrograde terms gives an elliptic motion. The motion of the pole, P(t), can therefore be described by its amplitude and phase of prograde and retrograde components as follows [70]: where A is amplitude of oscillation, α is a phase of oscillation, σ is frequency of the motion, t 0 is reference date, subscript p denotes prograde, and subscript r denotes retrograde. However, such representation is useful only for specific frequency, like, for example, annual (with frequency of 1 cycle per year), semiannual (with frequency of 2 cycles per year), or terannual (with frequency of 3 cycles per year). It is common in the study of PM excitation that the seasonal elliptical wave is represented as a combination of prograde and retrograde circular motion rather than equatorial χ 1 and χ 2 components related to the Earth-fixed system. Prograde and retrograde terms of seasonal variations are usually represented by their amplitudes and phases separately for annual, semiannual, and terannual oscillations [27,28,41,43,45,46,48,49,73]. Such amplitudes and phases are often shown as vectors plotted on phasor diagrams. Interpreting phasor diagrams is much easier than interpreting time series plots, which in the case of seasonal oscillations have the form of sinusoids. In phasor diagram, the length of a vector reflects a magnitude of amplitude, and the vector direction represents a phase, facilitating the simultaneous comparison of several seasonal series.
To interpret visually seasonal variations, apart from time series of equatorial components (χ 1 and χ 2 ) of seasonal oscillation which is the sum of annual, semiannual, and terannual variation, we also use phasor diagrams which display the amplitudes and phases of seasonal oscillation separately for prograde annual, retrograde annual, prograde semiannual, and retrograde semiannual terms. Such analysis is consistent with the approach taken in previous PM excitation research [27,28,41,43,45,46,48,49,73].
Finally, a detailed evaluation of each oscillation is performed based on comparison of HAM/CAM (or HAM) with GAO. The quality of HAM/CAM and HAM data is assessed by analysis of the following parameters: Error of standard deviation (STDerr), correlation coefficient (Corr) between HAM/CAM (or HAM) and GAO, percent of variance of GAO explained by HAM/CAM (or HAM) (relative explained variance, Var exp ), and root mean square deviation (RMSD). To statistically analyze the significance of Corr values and the significance of the differences between them, we compute the critical value of Corr and the standard error of the difference between the two Corr values for assumed confidence level (see details in [40]).
STDerr is obtained here from the following equation: The most favorable value of this parameter is 0%. The positive values indicate higher STD for HAM/CAM (or HAM) series than for GAO, whereas negative values indicate the opposite situation.
Var exp is computed using the following formula: and the most favorable value is equal to 100%. The centered RMSD is calculated according to the equation [82]: where f n and r n are modeled (HAM/CAM or HAM) and referenced (GAO) time series, respectively, and f and r are their mean values, respectively. The most favorable value of RMSD is 0.

Trends
We begin with a comparison of trends in GAO, GRACE-based HAM/CAM, and LSDM-based HAM ( Table 2). The signals from GIA are eliminated in all considered GRACE datasets and are also not present in LSDM; however, they remain in GAO. Therefore, for consistency of all trends, we remove the effect from post-glacial rebound from GAO trends using the ICE6G-D model described in [66]. Note that GAO trends with GIA signal retained are also reported for completeness; however, for objective comparison with HAM/CAM and HAM, only GAO trends with GIA removed should be considered.
Naturally, the removal of the GIA signal leads to higher consistency between GAO trends and HAM/CAM or HAM trends ( Table 2). The GIA assumption has a higher influence on the χ 2 trend than on the χ 1 trend. This is as expected, because signals from post-glacial rebound mainly affect continental areas of the Northern Hemisphere [83], especially Fennoscandia and the northern part of North America, which also contribute more to the χ 2 variability than the χ 1 variability. This has been previously shown in many papers concerning HAM/CAM analysis (e.g., [41,42,45,46,48,49]). All these studies concluded that χ 2 is more sensitive to mass changes over land areas, while χ 1 is more responsible to the variations over oceans and the Greenland Ice Sheet. Moreover, it has been shown that HAM/CAM is correspondingly better correlated with GAO for χ 2 than for χ 1 , and χ 2 is characterized by higher variability than χ 1 [39,42,46].  Table 2 shows that among all GRACE-derived HAM/CAM estimates, the series based on the GSFC solution have visibly smaller trends in both χ 1 and χ 2 than GAO and other HAM/CAM series. The trend difference between HAM/CAM from GSFC and HAM/CAM from other GRACE data centers may result from the use of different GIA models; CSR and JPL applied ICE6G-D [66], whereas GSFC used the previous version of this model, ICE6G [67]. However, it should be kept in mind that the GIA models are not perfect and some residual GIA errors are still anticipated. To test the influence of the GIA model on trends in χ 1 and χ 2 , we computed χ 1 and χ 2 trends due to GIA signal for both ICE6G and ICE6G-D models, and the differences between them (Table 3). It can be seen that the differences between GIA models do not explain entirely the differences in χ 1 , χ 2 trends observed for GSFC and other GRACE solutions. Thus, there are other factors influencing the trends in HAM/CAM. We suppose that the discrepancies in the determination of ice mass variations and their rates by various GRACE datasets may have some contribution to the trend difference in HAM/CAM. Table 3. Trends in the χ 1 and χ 2 components of HAM/CAM due to GIA signal obtained for ICE6G and ICE6G-D models, and the differences between the models. Among all HAM/CAM series, the most consistent trends with GAO are provided by JPL (for χ 1 ) and CSR v02 (for χ 2 ). However, for χ 2 , the trends provided by JPL, CSR v01, and CSR v02 are quite similar. All these solutions use the same GIA model, and some small differences in trends may result from different native resolution or different sizes of mascons. It is also clear that a correction due to ellipsoidal Earth, which is applied in CSR v02 data increases the HAM/CAM trend compared with trends obtained for CSR v01, especially in χ 1 . For HAM from LSDM, a trend in χ 1 is the lowest of all the considered series. This is probably because the LSDM does not contain signals from glaciated regions, especially Greenland, which is considered to noticeably contribute to the change in χ 1 . In recent years, an intensification of Greenland ice mass loss has been observed [15,16], which certainly has an impact on trends in PM excitation; however, this is not observed in LSDM-based HAM. In contrast, for χ 2 , the LSDM-based HAM trend is the most consistent with the trend observed for GAO. The finding that the χ 2 trend is better determined by hydrological model than the χ 1 trend is expected; however, it is quite surprising that the consistency in χ 2 trends is better for LSDM than for GRACE. Groundwater storage, which has more long-lasting variability than other TWS components [84], should have bigger impact on χ 2 trends than factors such as more rapidly changing soil moisture or water in biomass. However, while GRACE and GAO can reflect the whole contribution of groundwater, the LSDM models only a shallow layer of groundwater of a limited storage capacity [72].  (Figure 1). Given this time scale, this oscillation might be related to groundwater variation. The authors of [85] showed that groundwater variations affect mainly interannual TWS variability, so it should indirectly influence HAM/CAM as well. Other authors [86,87] pointed out that the climate changes affect long-term groundwater loss, which can contribute to the change in polar motion. The few-year change in HAM/CAM might also be related to a variation in continental water distribution due to different precipitation patterns, as the trend in continental water distribution between 2007-2013 was different than in other years [46]. A similar 4-5-year oscillation is also visible for χ 1 , but is not as apparent. This variation has also been reported in our previous research [45,46].  There is quite good phase consistency between HAM/CAM or HAM with GAO, but they slightly underestimate GAO amplitudes, as indicated by the negative values of STDerr ( Figure 2). The LSDM-based HAM reveals the highest STD consistency with GAO for χ 2 (STDerr of around 10%). This shows that LSDM can successfully model the water storage variability over land. For χ 1 , LSDM provides a similar level of STDerr as that obtained for the GRACE solutions.

Trend (Mas/Year)
The HAM/CAM data from the CSR v02 solution provide the highest consistency with GAO for χ 1 (Corr = 0.63, Var exp = 40% and RMSD = 6.21 mas), while HAM from LSDM is shown to be the least compatible with GAO for this component (Corr = 0.34, Var exp = 0% and RMSD = 8.02 mas). However, for χ 2 , all HAM/CAM and HAM results present a similar level of consistency with GAO (Corr ranges between 0.71 and 0.79, Var exp ranges between 53% and 58%, and RMSD ranges between 8.83 and 9.37 mas). It should be noted that the differences in Corr obtained for particular HAM/CAM and HAM are lower than errors of Corr determination, therefore, from statistical point of view, all series are correlated with GAO at the same level. Only LSDM-based HAM is distinguished by a lower Corr, but only for χ 1 . Overall, at this stage, it is difficult to choose which of the four mascon solutions is the most appropriate for analyzing the overall PM excitation induced by land hydrosphere and cryosphere as they present similar results.
The amplitude spectra of the series are dominated by annual oscillation, which is stronger for the retrograde term ( Figure 3). Figure 3 also indicates that the GAO series have clearly stronger semiannual and terannual oscillations compared with other series. This is also apparent from visual inspection of Figure 1.

Seasonal Variations
We now focus on analysis of seasonal oscillations visible in the spectra (Figure 3) of the abovementioned series. To present these oscillations, we show time series plots (Figure 4), which are supplemented with STDerr, Corr, Var exp , and RMSD diagrams ( Figure 5). A comparison of the time series indicates a good phase consistency between HAM/CAM and GAO and between HAM and GAO ( Figure 4) as indicated by high values of Corr ( Figure 5). Some phase shift relative to other time series is apparent only for LSDM-based HAM in χ 1 , which is supported by the lower Corr value. For χ 1 , HAM also exhibits the lowest amplitude of all series, probably because some of the ice mass is excluded in LSDM, and is reflected by the most negative STDerr. The parameters plotted in Figure 5 indicate that for seasonal variation in χ 1 , all mascon solutions provide a similar level of consistency with GAO, as they are characterized by congruous values of Corr (0.78-0.84), Var exp (51%-60%), and RMSD (2.58-2.66 mas). For χ 2 , however, HAM/CAM from the JPL solution has the highest consistency with GAO (Corr = 0.99, Var exp = 94%, RMSD = 1.90 mas). Nevertheless, the Corr values in χ 2 are statistically similar for all the series, because a maximum difference between two Corr values is lower than standard error of the difference between the two correlation coefficients.
It is worth noting that ellipsoidal correction, which is included in the CSR v02 solution, but is not applied in the previous CSR v01 release, increases seasonal consistency between HAM/CAM and GAO in χ 2 , but slightly deteriorates results for χ 1 . However, for Corr, this difference is below statistical significance (a difference in Corr for CSR v01 and CSR v02 is lower than the standard error of the difference between the two correlation coefficients).
In addition, LSDM-based HAM provides high consistency with GAO for χ 2 , which indicates that the LSDM can reasonably model the seasonal rainfall and soil moisture fluctuations that contribute to the seasonal variability of TWS and the resulting seasonal HAM/CAM variations. This might also indicate that the variations in deep groundwater, which are not modeled in LSDM, do not have a noticeable impact on seasonal HAM/CAM variation. This might be true, because changes in groundwater storage are more longlasting than seasonal changes in soil water and precipitation. Nevertheless, no strong agreement was found between HAM and GAO in χ 1 because LSDM does not include the entire cryosphere. More detailed information on amplitude and phase consistency among seasonal series are provided as phasor diagrams (Figure 6), which are plotted separately for annual and semiannual variations. Terannual oscillations are not considered in detail as it has been shown that their amplitudes are small compared with the annual and semiannual oscillations ( Figure 3). The annual variation dominates in the seasonal spectral band, which is mainly caused by the annual variability of precipitation and the resulting seasonality of soil moisture and snow cover changes ( Figure 6). For annual variation, the retrograde term is characterized by higher amplitudes than the prograde term, as reported in previous research [40,45,46] and is observed in amplitude spectra (Figure 3). For the annual retrograde oscillation, a higher consistency between GAO, HAM/CAM, and HAM is also observed. Of all GRACE-based HAM/CAM, the series obtained from JPL solution are most consistent with GAO in terms of amplitudes and phases for both prograde and retrograde annual variation. For semiannual oscillation, although HAM/CAM vectors remain quite consistent with each other, they are not in high agreement with the GAO indices. In particular, they visibly underestimate semiannual GAO amplitudes in the prograde term, while in the retrograde term they provide phases that are different than those of GAO. It is worth noting that LSDM enables a good amplitude consistency with GAO for annual retrograde (Figure 6b) and semiannual retrograde (Figure 6d) terms as the vectors are of similar length. Moreover, HAM based on LSDM is in the highest phase consistency with GAO for semiannual prograde term (Figure 6c). Nevertheless, the amplitude agreement for semiannual prograde oscillation is poor.  (a,b) is defined by the annual term as sin(2π(t − t 0 ) + φ), where t 0 is the reference epoch (1 January 2004). φ in (c,d) is defined by the semiannual term as sin(4π(t − t 0 ) + φ).

Non-Seasonal Variations
Finally, we analyze non-seasonal variations in GAO, HAM/CAM, and HAM. Time series of non-seasonal changes in GAO and HAM are presented in Figure 7, while STDerr, Corr, Var exp , and RMSD are shown in Figure 8. All HAM/CAM and HAM underestimate the variability of GAO, as indicated by the time series plots (Figure 7) and STDerr values ( Figure 8). This suggests that in the non-seasonal spectral band, there is an additional signal in GAO, which is not captured by GRACE and LSDM. This effect may be caused by the residual ocean signal resulting, for example, from the errors of ocean model used for GAO determination [47]. This effect is eliminated from GRACE by masking ocean areas, and is not present in LSDM, as the hydrological model delivers data from land only. The interannual variability is well captured by both HAM/CAM and HAM (Figure 7). In particular, a similar 4-5-year oscillation is observed for all series and is probably related to groundwater variation [45]. The best agreement with GAO is for HAM/CAM obtained from the CSR v02 mascon solution (Corr of 0.60 and 0.72 for χ 1 and χ 2 , respectively; Var exp of 34 and 50% for χ 1 and χ 2 , respectively; and RMSD of 5.62 and 8.17 mas for χ 1 and χ 2 , respectively). The inclusion of ellipsoidal correction in CSR v02, which is not present in CSR v01, improves the results more for χ 1 than for χ 2 . It is also apparent that HAM/CAM based on GSFC solution has the lowest consistency with GAO of all GRACE estimates. However, the difference compared with the results from the other mascon solutions is not large; all HAM/CAM are characterized by similar level of correlation with GAO. The reason for the slightly lower consistency with GAO obtained for GSFC-based HAM/CAM compared with the other mascon solutions might be related to the use of a different AOD dataset. In particular, GSFC uses a different ocean model than CSR and JPL, namely MOG2D (2 Dimensions Gravity Waves model, [88]) instead of MPIOM. The atmosphere model is the same for all three data centers (ECMWF model). It is also worth noting that, for χ 2 , HAM obtained from the LSDM provides almost the same level of Corr, Var exp , and RMSD as HAM/CAM from GRACE, while the STDerr is even more satisfactory than for the GRACE estimates. This shows that the LSDM is particularly useful for determining interannual water mass variations over continents, which is consistent with previous research [42,46].

Preliminary Estimates of HAM/CAM from GRACE-FO Mascon Data
We now present preliminary estimates of HAM/CAM processed using mascon data from the first two years of the new GRACE-FO mission and contrast them with GAO and HAM obtained from the LSDM. We examine mascon solutions only from JPL and CSR v02, because CSR v01 and GSFC data do not cover the GRACE-FO period (Figures 9 and 10). The time series plot indicates high consistency of HAM/CAM series obtained from both considered mascon solutions (Figure 9), which is supported by the similar parameter values shown in Figure 10. The HAM/CAM series is not consistent with GAO in χ 1 , which is reflected by the low Corr values. However, for χ 2 , the phase consistency between HAM/CAM and GAO is visibly better, including the same maximum peak around mid-2019 followed by a minimum in late 2019. There is no noticeable difference between results obtained from JPL and CSR v02 solutions. The latter solution only provides a slightly higher agreement between HAM/CAM and GAO, but this is not statistically significant. This might indicate that the use of ellipsoidal correction and ∆C 30 replacement, which are applied only in CSR v02, does not have a noticeable impact on HAM/CAM over short time periods. Overall, the current quality of HAM/CAM determined from GRACE-FO mascon data is as expected (Corr values vary from 0.17 to 0.57, Var exp values vary from −75% to 32%, RMSD values vary from 5.18 to 6.03 mas, depending on the solution and equatorial component considered). The consistency between HAM/CAM and GAO is expected to increase with increasing length of the GRACE-FO time series. However, even for the observations from the initial phase of the GRACE-FO mission, the gravimetric data performed better than the hydrological model in χ 2 . The χ 1 component of GRACE-FObased HAM/CAM currently shows lower agreement with GAO than χ 1 of LSDM-based HAM; however, as the GRACE-FO data lengthens, this compliance can be expected to improve.

Discussion
Analysis of HAM/CAM computed from mascon data provided by the three different processing centers indicated high consistency between estimates provided by JPL and CSR, and considerably lower agreement between the GSFC solution and other data. This is likely to be related to the strategy used for GSFC mascon computation, which is different from the methodology exploited by the CSR and JPL teams. The GSFC mascons, unlike JPL and CSR (both v01 and v02) mascons, are also not resampled to regular grids. The GSFC uses a different ocean model to determine de-aliasing data, and a different GIA model. Nevertheless, it was found that for each oscillation considered, there were no clear differences in consistency of GAO with HAM/CAM determined from various mascon solutions.
The ellipsoidal correction recently introduced in the newest mascon solution processed by the CSR team (CSR v02 solution) did not have a high impact on consistency between HAM/CAM and GAO, so it can be omitted in the studies related to PM variations. However, it may have a meaningful effect on length-of day variations since they are more closely related to the flattening of the Earth than PM. As GRACE and GRACE-FO datasets use a number of different corrections to improve the estimation of TWS distribution, future research should quantify the effect of each of them on PM excitation.
The main conclusion of this study is that, in general, all mascon solutions are equally suitable for determination of overall, seasonal, and non-seasonal HAM/CAM oscillations, and only trends are less accurately determined by GSFC than by other data centers. However, the inclusion of the newest GIA model in GSFC should partly improve the trend consistency between HAM/CAM and GAO. Therefore, we recommend all currently available mascon data for studies focusing on the interpretation of PM variations. At present, the GSFC team is working on their new mascon solution which, in contrast to the current one, will cover the time period of both GRACE and GRACE-FO missions and use the newest GIA model. Therefore, we might expect that this solution will perform better in PM excitation studies than the solution assessed in this paper. However, this hypothesis should be tested based on the validation with external PM data.

Data Availability Statement:
The data that support the findings of this study are available upon request from the corresponding author.