Evaluating Gravimetric Polar Motion Excitation Estimates from the RL06 GRACE Monthly-Mean Gravity Field Models

Over the last 15 years, the Gravity Recovery and Climate Experiment (GRACE) mission has provided measurements of temporal changes in mass redistribution at and within the Earth that affect polar motion. The newest generation of GRACE temporal models, are evaluated by conversion into the equatorial components of hydrological polar motion excitation and compared with the residuals of observed polar motion excitation derived from geodetic measurements of the pole coordinates. We analyze temporal variations of hydrological excitation series and decompose them into linear trends and seasonal and non-seasonal changes, with a particular focus on the spectral bands with periods of 1000–3000, 450–1000, 100–450, and 60–100 days. Hydrological and reduced geodetic excitation series are also analyzed in four separated time periods which are characterized by different accuracy of GRACE measurements. The level of agreement between hydrological and reduced geodetic excitation depends on the frequency band considered and is highest for interannual changes with periods of 1000–3000 days. We find that the CSR RL06, ITSG 2018 and CNES RL04 GRACE solutions provide the best agreement with reduced geodetic excitation for most of the oscillations investigated.


Introduction
Polar motion (PM) is affected by a wide range of processes with different temporal variability ranging from several days to many decades [1]. Such disturbances include the gravitational influence of other celestial bodies, continuously changing mass distribution in the Earth's surficial fluids (atmosphere, oceans and land hydrosphere), the effects of core-mantle coupling, and also groundwater depletion and ice mass loss resulting from recent climate changes [2][3][4].
It is commonly known that for time scales of a few years or less, the main contributors to variation in Earth's rotation are angular momentum changes of the solid Earth caused by mass redistribution in the Earth's surficial fluids. The excitation of polar motion due to changes in mass redistribution of atmosphere, ocean and land hydrosphere can be represented with effective angular momentum (EAM) functions: atmospheric, oceanic and hydrological. Information about mass transport in the Earth's fluid layers can be obtained from observations of temporal changes of the gravity field. Such measurements were provided from 2002 to 2017 by the Gravity Recovery and Climate Experiment (GRACE) mission [5]. The variations of degree-2, order-1 Stokes coefficients (∆C 21 , ∆S 21 ) of geopotential are commonly used in determination of mass-related PM excitation since they are proportional to the χ 1 and χ 2 components of the PM excitation [6]. The leading GRACE data centres at GeoForschungsZentrum (GFZ) in Potsdam, Germany; Center for Space Research (CSR) in Austin, USA; and the Jet Propulsion

Hydrological Signals in Observed Polar Motion-Geodetic Residuals
Observed or geodetic polar motion excitation (geodetic angular momentum, GAM) can be obtained from observed coordinates (x, y) of the pole by solving Liouville's equation [38,39]. In this research, the equatorial components (χ 1 and χ 2 ) of GAM are computed from x, y, which are taken from the daily combined 14 C04 series of Earth Orientation Parameters (EOP) derived from GNSS (Global Navigation Satellite Systems), SLR (Satellite Laser Ranging), and VLBI (Very Long Baseline Interferometry). EOP 14 C04 series are fully consistent with the International Terrestrial Reference Frame 2014 (ITRF 2014) [40,41]. The EOP data are provided by the International Earth Rotation and Reference System Service (IERS) (https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html) and updated on a regular basis with typically 30 days latency. In this study, for computation of GAM, we used an algorithm given at http://hpiers.obspm.fr/eop-pc/index.php?index=excitactive&lang=en.
Observed GAM are commonly reduced for the effects of atmosphere and ocean to study only the residual signals, denoted as geodetic residuals (GAM − AAM − OAM or GAO), using the following method: where GAM, AAM and OAM are geodetic angular momentum, atmospheric angular momentum and oceanic angular momentum, respectively. Note that here, AAM and OAM are effective angular momentum functions of atmosphere and ocean, which are sometimes abbreviated as EAM or EAMF of atmosphere and ocean. The resulting series are expected to represent the impact of land hydrosphere and barystatic sea-level changes on PM excitation. Nevertheless, it should be kept in mind that GAO will also contain signals from changes in the mass of ice sheets, the effects of large earthquakes [42,43], and the signatures of geomagnetic jerks [44].
The following data sets were used for this study to calculate GAO: • χ 1 and χ 2 equatorial components of GAM: obtained from the IERS website (https://www.iers.org/ IERS/EN/Home/home_node.html), with temporal resolution of 24 h; • χ 1 and χ 2 equatorial components of the AAM: provided by the GFZ, and based on ECMWF (European Center for Medium-Range Weather Forecasts) model [45]. The series are available at temporal resolution of 3 h and can be accessed from http://rz-vm115.gfz-potsdam.de:8080/ repository (Earth System Model GFZ (ESMGFZ) data). Note that the current AAM version provided by GFZ is consistent with GRACE AOD1B RL06 Atmosphere and Ocean De-Aliasing Level-1B Release-6 (AOD1B RL06) [35,36]; • χ 1 and χ 2 equatorial components of the OAM: provided by the GFZ, and based on MPIOM (Max Planck Institute Ocean Model) [46]. The series are available at temporal resolution of 3 h and can be accessed from http://rz-vm115.gfz-potsdam.de:8080/repository (ESMGFZ data). Note that the current OAM version provided by GFZ is consistent with GRACE AOD1B RL06.
The use of ESMGFZ data for computing GAO provides consistent results with GRACE estimations, because AOD1B RL06 and ESMGFZ were based on the same input models. Any other AAM and OAM data set would complicate the analysis by introducing unknown effects into the comparison.
It should be kept in mind that C04 EOP daily series suffer some weaknesses connected with aliasing of sub-daily signals caused by mis-modelled ocean tides. It is known that pole coordinates contain diurnal and sub-diurnal signals reaching 1 mas (miliarcsecond, 1 mas 3,09 cm on the Earth's surface) [41]. However, these high-frequency tidal signals cannot be included in 24-h C04 EOP series dedicated to providing reference Earth orientation data on a daily time rate. Nevertheless, in our study we focus on monthly PM excitation changes, so this should not affect our results.
The C04 EOP data have recently been reprocessed to ensure consistency with the new ITRF 2014 [47] and we used these updated series. The new daily EOP combined 14 C04 series were recomputed from 1984 according to improved algorithms and replaced the previous 08 C04 series in February 2017. The main reason for this reprocessing was that the differences between the former 08 C04 series and the EOP produced jointly with the ITRF 2014, starting from 2010, produced a bias in the y pole coordinate that reached 50 µas (microarcsecond) above the mean formal uncertainty estimated by the IGS (International GNSS Service) [41]. The 14 C04 EOP has been tied to the ITRF 2014 EOP solution and IVS (International VLBI Service for Geodesy and Astrometry) combination to provide consistency with the conventional ITRF 2014 [47] and ICRF2 (Second International Celestial Reference Frame) [48] reference frames.

Polar Motion Excitation Derived from GRACE
The relationship between ∆C 21 , ∆S 21 coefficients of the geopotential and the equatorial components of the PM excitation (χ 1 and χ 2 ) can be expressed as follows (based on [39] and Chapter 3.09.5 in [6]): where, ∆C 21 and ∆S 21 are the fully normalized Stokes coefficients of the Earth's gravity field; R e and M are the Earth's mean Earth's radius (6378136.6 m) and mass (5.9737 × 1024 kg), respectively; A = 8.0101 × 1037 kg·m 2 , B = 8.0103 × 1037 kg·m2, and C = 8.0365 × 1037 kg·m 2 are the Earth's principal moments of inertia, and A = (A+B)/2 is the average of the equatorial Earth's principal moments of inertia ( Table 1 in [6]). The PM excitation series obtained from the GRACE observations as represented by the monthly GRACE satellite-only models (GSM) are sometimes denoted as gravimetric-hydrological excitation, or simply as GRACE-based HAM. As indicated above, we note that the Level-2 monthly-mean gravity fields from GRACE contain not only the mass effects from terrestrial water storage, but also other signals such as the barystatic sea-level changes due to inflow of continental water into the oceans, GIA effects, or tectonic signals related to large earthquakes. We thus denoted the GRACE-based angular momentum series as GSMAM to emphasize that they were derived from the Level-2 GSM fields of GRACE. During the processing of GRACE data, the signals from short-term non-tidal atmospheric and oceanic variations were reduced by applying AOD1B [35,36]. Therefore, GSMAM series do not contain these effects.
In this study, we compared the following monthly times series of time-variable Level-2 GRACE GSM fields provided by five different processing groups: . Detailed information about considered solutions, with background models and data sources are given in Table A1 in Appendix A. Table A1 also contains corresponding  information for previous GRACE RL05 solutions.  It should be noted that the CNES RL04 and ITSG 2018 solutions were designed to be compatible  with official CSR RL06, JPL RL06 and GFZ RL06 data and were published around the same time,  so, we used the RL06 indication for all five solutions throughout the remainder of this paper. Regarding to the series from GFZ, it is now officially recommended to replace C 21 /S 21 coefficients for GFZ RL06 with SLR estimates [49]. These coefficients are estimated from a GRACE and SLR combination on the level of normal equations. Each monthly estimate is fully independent from the rest of the time series. Additionally, C 20 and C 30 coefficients need to be replaced according to this recommendation, but in this work we focus only on GSMAM computed from C 21 , S 21 . In the current study, for GFZ estimates, we consider GSMAM obtained from both the non-combined GFZ RL06 solution and the combined GFZ RL06 + SLR series. The latter are denoted in this paper as the GRACE + SLR series.

HAM and SLAM Processed by GFZ
We also considered HAM based on the Land Surface Discharge Model (LSDM) and accessed from GFZ (ftp://esmdata.gfz-potsdam.de/EAM). The LSDM simulates both vertical and lateral water transport and storage on land surfaces [52,53]. It is given in a regular 0.5 • × 0.5 • global latitude-longitude grid and with a 24 h time step. LSDM is forced with precipitation, evaporation, and surface temperature from the ECMWF model. The χ 1 , χ 2 LSDM-based HAM excitation series were taken directly from the GFZ website.
It should be remembered that while comparing HAM obtained from hydrological models with GSMAM derived from GRACE data, it is crucial to add the contribution of barystatic sea level changes due to inflow of water from lands into the oceans (sea-level angular momentum, SLAM) to model-based HAM. It is important, because both GRACE GSM fields and GAM (and also GAO which was computed after removing AAM and OAM from GAM) include SLAM, but it is not modeled in hydrological models. Therefore, in this study, we added SLAM to LSDM-based HAM. The SLAM series were based on ECMWF and LSDM data, provided by the GFZ (ftp://esmdata.gfz-potsdam.de/EAM) and available in the form of χ 1 , χ 2 [54].
Our motivation for choosing HAM from GFZ was that the consistent SLAM part of the excitation is publicly available only for the GFZ data set. This SLAM series are also compatible with GRACE AOD1B dealiasing data. To show how large the error would be if mass balance was not considered, we can give the relation of SLAM/(HAM + SLAM). The most preferred value of this dimensionless ratio is 1 (which means that SLAM and HAM + SLAM are the same). For the time period considered here, this value ranges from −9.0 to +3.7 for χ 1 and from −10.3 to +3.0 for χ 2, with standard deviations of 1.4 and 1.3 for χ 1 and χ 2 , respectively. The variability of SLAM is smaller than variability of HAM + SLAM as standard deviation of SLAM series represents 28% of standard deviation of HAM + SLAM for χ 1 and 19% for χ 2 . Another motivation to choose HAM based on LSDM were our previous results [15], which showed that compared to other hydrological and climate models, LSDM provided higher consistency with GAO and GRACE-based excitation series. In this study, we consider only the sum of HAM and SLAM, using the term HAM + SLAM throughout the remainder of this paper.

Results
In this paper, we focus on the overall analysis of month-to-month variations in time series (Section 3.1), month-to-month variations and trends in time series for four time periods characterized by different accuracy of GRACE observations (Section 3.2), secular trends, seasonal and non-seasonal oscillations (Section 3.3), and variations in four different frequency bands with periods of 60-100, 100-450, 450-1000, and 1000-3000 days (Section 3.4).

Overall Comparison of GSMAM Series
Initially, we display time series of GSMAM computed using different GRACE solutions and GRACE + SLR, without applying any filtering (Figure 1a,b). The series have quite similar character other oscillations of periods both shorter and longer than annual. It is obvious that GRACE + SLR series are less noisy than GSMAM from GFZ RL06 solution. Notably, after 2016, the series became more noisy, as the accuracy of GRACE decreased due to a change in thermal control of the spacecraft due to the increasingly limited battery capacity. In 2017, one accelerometer was turned off and that caused a further degradation of the results, especially visible for GFZ and ITSG solutions. Figure 1 gives also formal errors of estimated χ1, χ2, which were computed from the errors of ΔC21, ΔS21 coefficients ( Figure 1 c-f) (errors for CNES RL04 were not available). It is visible that the best GRACE performance occurs in the years 2006-2013, and before and after this period the error values are higher and noisier. During the entire 15-year GRACE operational period, some peaks of visibly higher errors occurs typically related to occasional periods of short repeat orbits. Compared to other GRACE solutions, the CSR RL06 series has the lowest and the most stable errors of χ1 and χ2. However, the errors of GRACE + SLR series are even lower. Notably, after 2005, the errors of GFZbased excitation series start to increase almost linearly.  Due to some gaps in the GRACE observations, different temporal sampling intervals (3 h for ECMWF and MPIOM models, 24 h for GAM, HAM from LSDM and SLAM, and only monthly for GRACE), and different periods in which these data are available, all time series were downsampled to monthly time-steps, using a Gaussian filter with full width at half maximum (FWHM) parameter of 60 days. In the reminder of this paper, we considered the same period for all data sets, namely between September 2002 and June 2016, except for the tests reported in Section 3.2 Figure 2 presents time series of χ1, χ2 components of hydrological PM excitation obtained from all considered GSMAM solutions as well as GRACE + SLR and GFZ HAM + SLAM, together with the GAO. To assess the variability of time series, we also computed their standard deviation (Table 1). In general, the new RL06 GRACE solutions underestimate the standard deviation of GAO particularly for χ2. In order to quantify the agreement between GAO and different GSMAM, we computed correlation coefficients and relative explained variance between evaluated and reference series ( Table  1). The relative explained variance (Varexp) describes the variability agreement between two time series, reference series (r) and evaluated series (e): where Var (r) , Var (e) and Var (r˗e) are variance of reference series, variance of evaluated series and variance of a difference between the reference and evaluated series, respectively. In this study, the reference series are GAO and the evaluated series are GSMAM, GRACE + SLR and HAM + SLAM. There is a better consistency of GSMAM series with GAO for the case of χ2 component, which has been already reported in previous research (e.g., [7,12,[15][16][17]30,55]), and is tightly related to the continent-ocean distribution. The values given in Table 1 reveal that CSR RL06, ITSG 2018 and CNES RL04 solutions provide the highest correlation coefficients with GAO and the most satisfactory relative explained variance values for both χ1 and χ2. Notably, HAM + SLAM obtained from models processed by GFZ, exhibits results comparable with those from CSR RL06, ITSG 2018 and CNES RL04, but only in χ2. GRACE + SLR series are highly more consistent with GAO than GFZ-based GSMAM, however, the results are not as satisfactory as for GSMAM from CSR RL06 data. Due to some gaps in the GRACE observations, different temporal sampling intervals (3 h for ECMWF and MPIOM models, 24 h for GAM, HAM from LSDM and SLAM, and only monthly for GRACE), and different periods in which these data are available, all time series were downsampled to monthly time-steps, using a Gaussian filter with full width at half maximum (FWHM) parameter of 60 days. In the reminder of this paper, we considered the same period for all data sets, namely between September 2002 and June 2016, except for the tests reported in Section 3.2. Figure 2 presents time series of χ 1 , χ 2 components of hydrological PM excitation obtained from all considered GSMAM solutions as well as GRACE + SLR and GFZ HAM + SLAM, together with the GAO. To assess the variability of time series, we also computed their standard deviation (Table 1). In general, the new RL06 GRACE solutions underestimate the standard deviation of GAO particularly for χ 2 . In order to quantify the agreement between GAO and different GSMAM, we computed correlation coefficients and relative explained variance between evaluated and reference series ( Table 1). The relative explained variance (Var exp ) describes the variability agreement between two time series, reference series (r) and evaluated series (e): where Var (r) , Var (e) and Var (r−e) are variance of reference series, variance of evaluated series and variance of a difference between the reference and evaluated series, respectively. In this study, the reference series are GAO and the evaluated series are GSMAM, GRACE + SLR and HAM + SLAM. There is a better consistency of GSMAM series with GAO for the case of χ 2 component, which has been already reported in previous research (e.g., [7,12,[15][16][17]30,55]), and is tightly related to the continent-ocean distribution. The values given in Table 1 reveal that CSR RL06, ITSG 2018 and CNES RL04 solutions provide the highest correlation coefficients with GAO and the most satisfactory relative explained variance values for both χ 1 and χ 2 . Notably, HAM + SLAM obtained from models processed by GFZ, exhibits results comparable with those from CSR RL06, ITSG 2018 and CNES RL04, but only in χ 2 . GRACE + SLR series are highly more consistent with GAO than GFZ-based GSMAM, however, the results are not as satisfactory as for GSMAM from CSR RL06 data.  In general, polar motion excitation shows two circular motions with opposite directions: clockwise, usually called retrograde, with frequency marked as negative, and counter-clockwise, named prograde, with frequency marked as positive [7,56,57]. It is common in the study of polar motion that the seasonal elliptical wave is represented as a combination of prograde and retrograde circular motions, rather than the equatorial χ1 and χ2 components related to the Earth-fixed system. This method of presentation facilitates the physical interpretation of seasonal PM excitation because the transfer function between the polar motion and polar excitation function is asymmetrical with respect to zero frequency. Such prograde and retrograde seasonal variations in PM excitation are  Table 1. Standard deviation of time series, correlation coefficients with GAO, and percentage of variance in GAO explained by GSMAM, GRACE + SLR and HAM + SLAM. Before computing the values, linear trends were removed from the series. Note that the critical value of the correlation coefficient for 54 independent points and the confidence level of 0.95 was 0.23. The standard error of the difference between two correlation coefficients for 54 independent points was 0.20. The highest values of correlations and relative explained variances are italicized.

Correlation Coefficient with GAO
Relative Explained Variance (%) In general, polar motion excitation shows two circular motions with opposite directions: clockwise, usually called retrograde, with frequency marked as negative, and counter-clockwise, named prograde, with frequency marked as positive [7,56,57]. It is common in the study of polar motion that the seasonal elliptical wave is represented as a combination of prograde and retrograde circular motions, rather than the equatorial χ 1 and χ 2 components related to the Earth-fixed system. This method of presentation facilitates the physical interpretation of seasonal PM excitation because the transfer function between the polar motion and polar excitation function is asymmetrical with respect to zero frequency. Such prograde and retrograde seasonal variations in PM excitation are usually represented by amplitudes and phases of annual, semi-annual and ter-annual oscillations [7,8,[16][17][18]30,45]. However, it is also possible to determine prograde and retrograde terms of PM excitation without focusing on specific frequency [56,57].
It is known that seasonal oscillation is the strongest signal in hydrological excitation of PM but it is possible to check the magnitude of other spectral bands by plotting amplitude spectra of the complex-valued components of χ 1 + iχ 2 . To check the oscillation amplitude distribution in GSMAM, GAO, GRACE + SLR and HAM + SLAM, we present such kind of spectra for all considered excitation series in Figure 3. The spectra are computed using a Fourier Transform Band Pass Filter (FTBPF) [58] in the broadband with a 60-and 3000-day cut-off ( Figure 3a) and in the broadband with 60-and 450-day cut-off ( Figure 3b). It should be kept in mind that the oscillations with periods shorter than 60 days cannot be considered, because all series were monthly sampled and their Nyquist frequency is about 6 cycles per year (cpy) that corresponds to the periods of about 60 days. The spectra show that there are several amplitude peaks in a band of 60-450 days with annual signal being the strongest (Figure 3a,b). For GAO and most of GSMAM series, some longer oscillations are also visible, namely 450-1000 days in prograde band, 450-700 days in retrograde band and 1000-3000 days in both prograde and retrograde bands ( Figure 3a). Figure 3b shows that the spectra of all considered series exhibit annual and semi-annual oscillations, both in the prograde and retrograde bands. The retrograde oscillations of the GAO have greater amplitudes than any of the considered GSMAM series (Figure 3a). In the prograde band, for the periods between 60 and 1000 days, GAO series is the most powerful, but for longer oscillations, GSMAM from GFZ RL06 reveals to have the highest amplitudes. The spectra obtained from HAM + SLAM has almost the weakest amplitudes for the whole prograde spectral band and also the retrograde semi-annual oscillations are the least energetic. However, HAM + SLAM captures the strong annual retrograde peak of GAO the best.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 29 usually represented by amplitudes and phases of annual, semi-annual and ter-annual oscillations [7,8,[16][17][18]30,45]. However, it is also possible to determine prograde and retrograde terms of PM excitation without focusing on specific frequency [56,57]. It is known that seasonal oscillation is the strongest signal in hydrological excitation of PM but it is possible to check the magnitude of other spectral bands by plotting amplitude spectra of the complex-valued components of χ1 + iχ2. To check the oscillation amplitude distribution in GSMAM, GAO, GRACE + SLR and HAM + SLAM, we present such kind of spectra for all considered excitation series in Figure 3. The spectra are computed using a Fourier Transform Band Pass Filter (FTBPF) [58] in the broadband with a 60-and 3000-day cut-off ( Figure 3a) and in the broadband with 60-and 450day cut-off ( Figure 3b). It should be kept in mind that the oscillations with periods shorter than 60 days cannot be considered, because all series were monthly sampled and their Nyquist frequency is about 6 cycles per year (cpy) that corresponds to the periods of about 60 days. The spectra show that there are several amplitude peaks in a band of 60-450 days with annual signal being the strongest (Figure 3a,b). For GAO and most of GSMAM series, some longer oscillations are also visible, namely 450-1000 days in prograde band, 450-700 days in retrograde band and 1000-3000 days in both prograde and retrograde bands ( Figure 3a). Figure 3b shows that the spectra of all considered series exhibit annual and semi-annual oscillations, both in the prograde and retrograde bands. The retrograde oscillations of the GAO have greater amplitudes than any of the considered GSMAM series ( Figure 3a). In the prograde band, for the periods between 60 and 1000 days, GAO series is the most powerful, but for longer oscillations, GSMAM from GFZ RL06 reveals to have the highest amplitudes. The spectra obtained from HAM + SLAM has almost the weakest amplitudes for the whole prograde spectral band and also the retrograde semi-annual oscillations are the least energetic. However, HAM + SLAM captures the strong annual retrograde peak of GAO the best.

Analysis of GSMAM Series for Different Time Periods
As it was noticed in Sect. 3.1, throughout the duration of the GRACE mission, the accuracy of its observations has not been stable, and this could affect the GSMAM quality assessment (see Figure 1). Moreover, in our previous work [59] we demonstrated that the trends in GAO depended strongly on considered time period. Therefore, this section evaluates GSMAM, GRACE + SLR and HAM + SLAM in several different time periods. Based on the information about the errors of GSMAM (Figure 1),

Analysis of GSMAM Series for Different Time Periods
As it was noticed in Section 3.1, throughout the duration of the GRACE mission, the accuracy of its observations has not been stable, and this could affect the GSMAM quality assessment (see Figure 1). Moreover, in our previous work [59] we demonstrated that the trends in GAO depended strongly on considered time period. Therefore, this section evaluates GSMAM, GRACE + SLR and HAM + SLAM in several different time periods. Based on the information about the errors of GSMAM (Figure 1 Table 2); the test of GSMAM, GRACE + SLR and HAM + SLAM compatibility with GAO, conducted using correlation coefficients (Table 3) and relative explained variances (Table 4).    The values of STD of time series given in Table 2   In general, the consistency between GSMAM and GAO, described by correlation coefficients (Table 3)       In terms of χ 1 , apart from 2002-2006 period, GAO is characterized by positive trend (Figure 4). As it was expected, for the period of the highest GRACE accuracy, the GSMAM trends in χ 1 are very compatible with each other and with GAO trend (Figure 1b). For this period, the maximum difference between GAO trend and GSMAM trend reaches 1.75 mas/year (for CSR RL06), while the minimum difference is equal to 0.79 mas/year (for GFZ RL06). For each considered time interval, GSMAM overestimates χ 1 GAO trend while HAM + SLAM underestimates it. For χ 2 , the GAO trend signs are less steady ( Figure 5). As it can be seen from Figure 5, χ 2 GSMAM trends are the most consistent with GAO trends and with each other for the initial period of the GRACE mission (2002)(2003)(2004)(2005)(2006). For this time interval, the maximum difference between GAO trend and GSMAM trend reaches 1.48 mas/year (for GFZ RL06), while the minimum difference is equal to 0.46 mas/year (for CSR RL06). Notably, in the period of the lowest GRACE observations errors (2006-2014, see Figure 1c,d) GSMAM trend signs in χ 2 are inconsistent with GAO trend sign (Figure 5b). We can generally conclude that the trend of HAM + SLAM series processed by GFZ is rather incompatible with trend of GAO for both χ 1 and χ 2 . Note that trends for the period 2017 until 2017.5 are reported for completeness, but not discussed due to the limited data coverage in that period.
The values of STD of time series given in Table 2  In general, the consistency between GSMAM and GAO, described by correlation coefficients (Table 3), is the highest for 2006-2014 and the lowest for 2017-2017.5. The exception are GSMAM series obtained from CSR RL06 and JPL RL06 solutions, for which we observe very high correlations with GAO during the last operational period of GRACE satellites, but only for χ 2 . However, only a six-month period of data cannot allow objective conclusions to be drawn. In the initial period of GRACE observations (2002)(2003)(2004)(2005)(2006), the biggest consistency with GAO for both χ 1 and χ 2 is provided by GSMAM from CSR RL06, while JPL RL06 gives satisfactory results for χ 1 only and ITSG 2018 provides high correlation for χ 2 only. The period of 2006-2014 is characterized by the highest correlation agreement between GSMAM and GAO, with the best results for CSR RL06 (for both χ 1 and χ 2 ), ITSG 2018 (for both χ 1 and χ 2 ) and CNES RL04 (for χ 2 only). After 2014, the correlation values start dropping and they reach their maximum for CSR RL06 (for both χ 1 and χ 2 ) and ITSG 2018 (for χ 1 only) in the period of 2014-2017. After turning off the accelerometer on one of the GRACE satellites in 2017, the correlations with GAO remain positive for both χ 1 and χ 2 only for GSMAM from CSR RL06. This is confirmed by the lowest errors for CSR data compared to other GRACE solutions (see Figure 1c,d). It is worth mentioning that HAM + SLAM series computed from numerical models processed by GFZ provides high correlation with GAO for all considered time periods particularly for χ 2 . This indicates that residual error of the ocean circulation remains both in GSM and GAO that are not explained by SLAM. This means that it requires further efforts to improve in particular the oceanic non-tidal background model for both GRACE and the EOP.
The analysis of relative explained variances (Table 4)

Analysis of Trends and Seasonal and Non-Seasonal Oscillation
We now focus on several analyses of GAO, GSMAM, GRACE + SLR and HAM + SLAM conducted for the same period, namely between September 2002 and June 2016. The most common methods to compare PM excitation series from different data sources involve either comparing full time series [12,14,30,37,[60][61][62] or decomposing them into seasonal and non-seasonal oscillations [8,15,17,18,45]. We now decompose all series into linear trends, and seasonal and non-seasonal changes. To calculate linear trends and seasonal oscillations, we use a least-squares method to fit a seasonal model to the series. The model comprises the first-degree polynomial and the sum of sinusoids with the periods of 365.25, 182.625, and 121.75 days [8]. In order to analyze seasonal variations, we remove the trends. Non-seasonal changes are obtained after removing trends and annual, semi-annual, and ter-annual signals. Figure 6 presents a comparison of linear trends in all series given in mas/year. We can observe that GSMAM trends for both χ 1 and χ 2 agree well with each other. They are also consistent with GAO trends in terms of their signs. The maximum difference in relation to GAO χ 1 trend is obtained for CSR RL06 (1.41 mas/year) and GRACE + SLR (1.46 mas/year), while the smallest one is found for GFZ RL06 (0.13 mas/year). The maximum difference in relation to GAO χ 2 trend is provided by GFZ RL06 (3.08 mas/year), and the smallest one is observed for JPL RL06 (1.60 mas/year). In general, GSMAM series overestimate negative trends of GAO in χ 2 , but for χ 1 they provide trends more consistent with GAO in terms of trend magnitude. The χ 1 trend detected in the HAM + SLAM is very low compared to GSMAM and GAO trends, certainly because the LSDM contains signals from nonglaciated regions only. We can therefore assume that the trend in GAO χ 1 is mainly due to Greenland ice mass loss, which is not included in the LSDM. Conversely, the trend in GAO χ 2 might be more closely related to the mass redistributions on land of the Northern Hemisphere. Nevertheless, χ 2 trends in HAM + SLAM are highly inconsistent with both GAO and GSMAM. These differences might result from the different treatment of groundwater changes by GAO, GRACE observations and LSDM. The whole contribution of groundwater variations is retained in GAO and GRACE data but LSDM models only a single shallow groundwater aquifer of rather limited storage capacity [54]. that GSMAM trends for both χ1 and χ2 agree well with each other. They are also consistent with GAO trends in terms of their signs. The maximum difference in relation to GAO χ1 trend is obtained for CSR RL06 (1.41 mas/year) and GRACE + SLR (1.46 mas/year), while the smallest one is found for GFZ RL06 (0.13 mas/year). The maximum difference in relation to GAO χ2 trend is provided by GFZ RL06 (3.08 mas/year), and the smallest one is observed for JPL RL06 (1.60 mas/year). In general, GSMAM series overestimate negative trends of GAO in χ2, but for χ1 they provide trends more consistent with GAO in terms of trend magnitude. The χ1 trend detected in the HAM + SLAM is very low compared to GSMAM and GAO trends, certainly because the LSDM contains signals from nonglaciated regions only. We can therefore assume that the trend in GAO χ1 is mainly due to Greenland ice mass loss, which is not included in the LSDM. Conversely, the trend in GAO χ2 might be more closely related to the mass redistributions on land of the Northern Hemisphere. Nevertheless, χ2 trends in HAM + SLAM are highly inconsistent with both GAO and GSMAM. These differences might result from the different treatment of groundwater changes by GAO, GRACE observations and LSDM. The whole contribution of groundwater variations is retained in GAO and GRACE data but LSDM models only a single shallow groundwater aquifer of rather limited storage capacity [54].  Figure 7. For χ1 seasonal oscillations, GSMAM has quite smaller amplitudes than GAO and HAM + SLAM. In the case of χ2, amplitude differences between all GSMAM and GAO are even more visible. Interestingly, the χ2 HAM + SLAM series from the numerical models has almost the same magnitude as GAO.  Figure 7. For χ 1 seasonal oscillations, GSMAM has quite smaller amplitudes than GAO and HAM + SLAM. In the case of χ 2 , amplitude differences between all GSMAM and GAO are even more visible. Interestingly, the χ 2 HAM + SLAM series from the numerical models has almost the same magnitude as GAO. A common way to analyze seasonal oscillations is to present their phasor diagrams showing amplitudes and phases of prograde and retrograde terms of annual, semi-annual, and ter-annual oscillations related to reference epoch (in our study the reference epoch is January 2004). These phasor diagrams for annual and semi-annual variations are shown in Figures 8 and 9, respectively, while the values of amplitudes and phases of seasonal (annual, semi-annual and ter-annual) oscillations and their errors are given in Tables A2-A4 (Appendix A). In general, the length of each vector represents A common way to analyze seasonal oscillations is to present their phasor diagrams showing amplitudes and phases of prograde and retrograde terms of annual, semi-annual, and ter-annual oscillations related to reference epoch (in our study the reference epoch is January 2004). These phasor diagrams for annual and semi-annual variations are shown in Figures 8 and 9, respectively, while the values of amplitudes and phases of seasonal (annual, semi-annual and ter-annual) oscillations and their errors are given in Tables A2-A4 (Appendix A). In general, the length of each vector represents the magnitude of amplitude, while the vector direction shows a phase. The greater a difference between two vector directions, the more the time series are shifted relative to each other. The phase shift between two time series can affect the size of the correlation. The congruous length of two vectors indicates that these two time series are characterized by consistent magnitude of amplitudes.  To analyze the agreement between non-seasonal GSMAM, GRACE + SLR, HAM + SLAM and GAO changes, we show time series of these oscillations (Figure 10), STD of the series, correlation coefficients and relative explained variance (Table 5). Figure 10 shows that there are some discrepancies between GSMAM from individual GRACE solutions, but the same amplitude peaks can be observed for all of them. Figure 10 also highlights a visible quasi-5-year oscillation in the χ2 component. Noticeably, the non-seasonal HAM + SLAM variations obtained from the geophysical models are rather smooth when compared to GRACE and GAO. However, the HAM + SLAM series represent the inter-annual variability in χ2 quite well. The values of STD given in Table 5 show that  To analyze the agreement between non-seasonal GSMAM, GRACE + SLR, HAM + SLAM and GAO changes, we show time series of these oscillations (Figure 10), STD of the series, correlation coefficients and relative explained variance (Table 5). Figure 10 shows that there are some discrepancies between GSMAM from individual GRACE solutions, but the same amplitude peaks can be observed for all of them. Figure 10 also highlights a visible quasi-5-year oscillation in the χ2 component. Noticeably, the non-seasonal HAM + SLAM variations obtained from the geophysical models are rather smooth when compared to GRACE and GAO. However, the HAM + SLAM series represent the inter-annual variability in χ2 quite well. The values of STD given in Table 5 show that GSMAM series underestimate a variability of GAO in χ2, while STD of HAM + SLAM for this component is quite consistent with STD of GAO. The highest correlations with GAO for both χ1 and Naturally the annual oscillations ( Figure 8) are stronger than semi-annual ones (Figure 9). The amplitude differences between GSMAM and GAO series shown in Figure 7 are reflected in inconsistencies in prograde and retrograde amplitudes (Figures 8 and 9). Notably, HAM + SLAM and each GSMAM series underestimate annual and semi-annual amplitudes of GAO. For annual prograde term, the highest amplitude agreement with GAO is provided by GSMAM from CSR RL06 and JPL RL06, while for annual retrograde part the highest amplitude consistency with GAO is obtained for HAM + SLAM. Apart from GFZ RL06 in prograde and CNES RL04 in retrograde band, GSMAM annual phases are rather consistent with each other. Notably, GSMAM from CNES RL04 solution provides high amplitude and phase consistency in semi-annual retrograde part, but for semi-annual prograde term this agreement is one of the lowest. GSMAM from ITSG 2018 data reveals to provide the highest compatibility with GAO for semi-annual prograde part. For ter-annual oscillations, in general, the amplitude and phase values are highly questionable as in many cases the errors exceed the amplitude values (see Table A4 in Appendix A). Therefore, we do not consider ter-annual variations. The above findings highlight that it is difficult to conclude which GRACE solution provides the highest seasonal consistency with GAO as the results depend on what oscillations are considered (annual, semi-annual, or ter-annual; prograde or retrograde; amplitude or phase).
To analyze the agreement between non-seasonal GSMAM, GRACE + SLR, HAM + SLAM and GAO changes, we show time series of these oscillations (Figure 10), STD of the series, correlation coefficients and relative explained variance (Table 5). Figure 10 shows that there are some discrepancies between GSMAM from individual GRACE solutions, but the same amplitude peaks can be observed for all of them. Figure 10 also highlights a visible quasi-5-year oscillation in the χ 2 component. Noticeably, the non-seasonal HAM + SLAM variations obtained from the geophysical models are rather smooth when compared to GRACE and GAO. However, the HAM + SLAM series represent the inter-annual variability in χ 2 quite well. The values of STD given in Table 5 show that GSMAM series underestimate a variability of GAO in χ 2 , while STD of HAM + SLAM for this component is quite consistent with STD of GAO. The highest correlations with GAO for both χ 1 and χ 2 components are obtained for CNES RL04, but ITSG 2018 and CSR RL06 data also provide significant correlations which exceed 0.56 for χ 1 and 0.64 for χ 2 . Accordingly, only these three solutions provide positive relative explained variances for both χ 1 and χ 2 and the best result is obtained for CNES RL04. Noticeably, the modelled HAM + SLAM series reveals sufficient correlation with GAO and high relative explained variance for χ 2 .

Analysis of GSMAM Series in Different Spectral Bands
We now extend our assessment to different spectral bands: (1) short-term changes with a period of 60-100 days, (2) a 100-450-day spectral band that includes periods of seasonal oscillations, and two longer-period spectral bands, namely (3) 450-1000 days, and (4) 1000-3000 days. This analysis aims to determine whether the level of agreement between GSMAM and GAO, between GRACE + SLAM and GAO or between HAM + SLAM and GAO depends on a specific spectral band in particular. The choice of spectral bands is determined by the length of available GRACE solutions, their monthly sampling, and the results of analysis of amplitude spectra of GAO, GSMAM and HAM + SLAM (Figure 3). These spectra show that, apart from strong seasonal oscillations, some longer-period variations occurs, namely with periods of around 450-1000 days and 1000-3000 days. Therefore, in the remainder of this paper, we consider these two spectral bands. Although we clearly indicated amplitude peaks for annual and semi-annual signals (see Figure 3), we decided to include them in one spectral band (100-450 days), since these changes were analyzed in detail already in Section 3.3. However, it should be noted that the 100-450-day range, apart from seasonal changes, also contains oscillations of both shorter and longer periods. Therefore, we cannot fully identify this spectral band with the seasonal changes analyzed in the previous chapter. To analyze short-term PM excitation, we introduce a spectral band of 60-100 days. It should be kept in mind that oscillations of higher frequencies cannot be considered due to monthly resolution of GRACE data. In this study, the oscillations in different spectral bands are calculated using a bandpass 4-pole Butterworth filter [63]. Figure 11a,b show long-term changes (with periods of 1000-3000 days) of GAO, GSMAM, GRACE + SLR and HAM + SLAM series, while the STD of time series, correlation coefficients with GAO and variance explained are presented in Table 6. Figure 11a shows that the χ 1 series of HAM + SLAM is almost flat but the HAM + SLAM χ 2 component shows the highest amplitude of all considered time series. Each GSMAM series underestimates STD of GAO for both χ 1 and χ 2 (Table 6). Nevertheless, there is good phase correspondence in χ 2 between all considered data. For GAO, GSMAM and HAM + SLAM, we observe the same oscillation with the minimum amplitude peak at the end of 2007 and maximum amplitude peak at the beginning of 2011. Similar to non-seasonal variations (Figure 10), a visible quasi-5-year oscillation in the χ 2 component is detected. In terms of correlation with GAO, the highest values for both χ 1 and χ 2 in the 1000-3000-day spectral band is obtained for ITSG 2018 and CSR RL06 data (Table 6). However, it should be noted that all GSMAM series perform very well in the case of χ 2 . As similar to the correlation results, the highest relative explained variances are observed for GSMAM from ITSG 2018 and CSR RL06 solutions (Table 6). Notably, the variance of the χ 1 GAO component is explained by GRACE less accurately than χ 2 . However, all variances are positive for χ 1 and only GFZ RL06 and CNES RL04 produce values that do not exceed 55%. the end of 2007 and maximum amplitude peak at the beginning of 2011. Similar to non-seasonal variations (Figure 10), a visible quasi-5-year oscillation in the χ2 component is detected. In terms of correlation with GAO, the highest values for both χ1 and χ2 in the 1000-3000-day spectral band is obtained for ITSG 2018 and CSR RL06 data (Table 6). However, it should be noted that all GSMAM series perform very well in the case of χ2. As similar to the correlation results, the highest relative explained variances are observed for GSMAM from ITSG 2018 and CSR RL06 solutions (Table 6). Notably, the variance of the χ1 GAO component is explained by GRACE less accurately than χ2. However, all variances are positive for χ1 and only GFZ RL06 and CNES RL04 produce values that do not exceed 55%. Figure 11. Changes with the periods of 1000-4000 (a,b), 450-1000 (c,d), 100-450 (e,f), and 60-100 (g,h) days in χ1 and χ2 components of GAO, GSMAM, GRACE + SLR, and HAM + SLAM. Table 6. Standard deviation of GSMAM, GAO, GRACE + SLR, and HAM + SLAM time series, correlation coefficients with GAO, and percentage of variance in GAO explained by GSMAM, GRACE + SLR and HAM + SLAM. Oscillations with periods 1000-3000 days are considered. The critical value Figure 11. Changes with the periods of 1000-4000 (a,b), 450-1000 (c,d), 100-450 (e,f), and 60-100 (g,h) days in χ 1 and χ 2 components of GAO, GSMAM, GRACE + SLR, and HAM + SLAM. Table 6. Standard deviation of GSMAM, GAO, GRACE + SLR, and HAM + SLAM time series, correlation coefficients with GAO, and percentage of variance in GAO explained by GSMAM, GRACE + SLR and HAM + SLAM. Oscillations with periods 1000-3000 days are considered. The critical value of the correlation coefficient for 9 independent points and a confidence level of 0.95 is equal to 0.58. The standard error of the difference between two correlation coefficients for the 9 independent points is equal to 0.58. The highest values of correlations and relative explained variances are italicized. For the periods between 450 and 1000 days, we found that all series are rather consistent with each other in terms of amplitudes (apart from χ 1 for HAM + SLAM), and the same amplitude peaks are visible for all of them (Figure 11c,d). Similar to the 1000-3000-day spectral band, χ 1 of HAM + SLAM is distinguished by lower amplitudes than χ 1 of GAO and GSMAM. Quasi-3-year changes in χ 2 are visible in all series, which might be the result of groundwater variations. Compared to correlations with GAO and relative explained variances obtained for the 1000-3000-day spectral band, the values for 450-1000-day-changes are lower ( Table 7). The best correlation and variance agreement with GAO is detected for GSMAM from CNES RL04 solution. Table 7. Standard deviation of GSMAM, GAO, GRACE + SLR, and HAM + SLAM time series, correlation coefficients with GAO, and percentage of variance in GAO explained by GSMAM, GRACE + SLR and HAM + SLAM. Oscillations with periods 450-1000 days are considered. The critical value of the correlation coefficient for 20 independent points and a confidence level of 0.95 is equal to 0.38. A standard error of a difference between two correlation coefficients for 20 independent points is equal to 0.34. The highest values of correlations and relative explained variances are italicized. Figure 11e,f show GAO, GSMAM, GRACE + SLR and HAM + SLAM changes with periods between 100 and 450 days. Because these variations include, among others, also annual, semi-annual, and ter-annual oscillations, the series resemble sinusoids but with amplitudes of various sizes. Compared with longer oscillations, 100-450-day changes are noticeably stronger (note that the scale on the y axis is twice as large as that in Figure 11a-d). This is not surprising as seasonal oscillations have previously been shown to be the most powerful in hydrological PM excitation [15,28]. GSMAM and HAM + SLAM series reveal to underestimate the GAO variability since their STD values are lower than for the reference series, which is particularly evident in the case of χ 2 ( Table 8). The comparison of correlation coefficients with GAO and relative explained variances given in Table 8 shows that the highest agreement with GAO, for both χ 1 and χ 2 components, is observed for GSMAM from CSR RL06 solution, but satisfactory results are also noticeable for GSMAM from CNES RL04 and ITSG 2018. Notably, in χ 2 , HAM + SLAM series reveals to provide the highest consistency with GAO of all considered time series, with correlation of 0.76 and relative explained variance of 57%. Table 8. Standard deviation of GSMAM, GAO, GRACE + SLR, and HAM + SLAM time series, correlation coefficients with GAO, and percentage of variance in GAO explained by GSMAM, GRACE + SLR and HAM + SLAM. Oscillations with periods 100-450 days are considered. The critical value of the correlation coefficient for 59 independent points and a confidence level of 0.95 is equal to 0.22. A standard error of a difference between two correlation coefficients for 59 independent points is equal to 0.19. The highest values of correlations and relative explained variances are italicized. Figure 11e,f show GAO, GSMAM, GRACE + SLR and HAM + SLAM changes with periods between 100 and 450 days. Because these variations include, among others, also annual, semi-annual, and ter-annual oscillations, the series resemble sinusoids but with amplitudes of various sizes. Compared with longer oscillations, 100-450-day changes are noticeably stronger (note that the scale on the y axis is twice as large as that in Figure 11a-d). This is not surprising as seasonal oscillations have previously been shown to be the most powerful in hydrological PM excitation [15,28]. GSMAM and HAM + SLAM series reveal to underestimate the GAO variability since their STD values are lower than for the reference series, which is particularly evident in the case of χ 2 ( Table 8). The comparison of correlation coefficients with GAO and relative explained variances given in Table 8 shows that the highest agreement with GAO, for both χ 1 and χ 2 components, is observed for GSMAM from CSR RL06 solution, but satisfactory results are also noticeable for GSMAM from CNES RL04 and ITSG 2018. Notably, in χ 2 , HAM + SLAM series reveals to provide the highest consistency with GAO of all considered time series, with correlation of 0.76 and relative explained variance of 57%.

Series
The short-term GAO and GSMAM oscillations (with periods of 60-100 days) are shown in Figure 11g,h and the STD of the series, correlations with GAO, and relative explained variances are given in Table 9. The short-term variations reveal the smallest STD of all considered oscillations. Neither GSMAM nor HAM + SLAM series provide amplitudes and STD as high as for GAO. Notably, for GAO and GSMAM series there are some periods of stronger and some periods of weaker amplitudes, which is not observed for HAM + SLAM. The latter is characterized by the smallest amplitudes and STD among all the data considered. Taking into consideration correlation coefficients with GAO and relative explained variances (Table 9), we found that the best agreement with GAO for both χ 1 and χ 2 components is obtained when ITSG 2018, CNES RL04, or CSR RL06 are used. In general, it can be concluded that short-period oscillations can be described by GRACE with lower accuracy than other variations as correlation and relative explained variances are smaller than for longer-period changes. One of the reasons for the worse consistency between GAO and GSMAM in the short-term spectral band than in other considered bands, may be an aliasing of sub-daily signals in C04 EOP daily series. Table 9. Standard deviation of GSMAM, GAO, GRACE + SLR, and HAM + SLAM time series, correlation coefficients with GAO, and percentage of variance in GAO explained by GSMAM, GRACE + SLR and HAM + SLAM. Oscillations with periods 60-100 days are considered. The critical value of the correlation coefficient for 112 independent points and a confidence level of 0.95 is equal to 0.15. A standard error of a difference between two correlation coefficients for 112 independent points is equal to 0.14. The highest values of correlations and relative explained variances are italicized.

Discussion
The study reveals that, although the previous researches reported an improvement in consistency between particular GRACE solutions for RL06 compared to former RL05 release, and an increase of the agreement between GSMAM and GAO [37,57], the differences between particular GSMAM series are still remarkable. One of the sources of these discrepancies are different background models used to develop particular GRACE solutions. The biggest update in RL06 data concerns the change of mean pole model from cubic to linear in all solutions, which mostly affected GSMAM trends, but also might have a minor impact on long-term variability. In fact, the GSMAM trends from RL06 are more consistent with each other and with GAO trends than GSMAM trends from RL05 [15]. All the considered GRACE solutions use the same solid Earth tides model based on IERS Conventions (2010) [34] (see Table A1 in Appendix A). Each GRACE RL06 temporal model, except CNES RL04, introduced updated dealiasing product, namely AOD1B RL06 [35,36] instead of AOD1B RL05, which certainly benefited GSMAM in the short-term spectral band. However, while GSMAM from CSR RL06 and ITSG 2018 provide similar level of agreement with GAO for 60-100-day spectral band, GFZ RL06 and JPL RL06 reveal visibly lower correlations and relative explained variances for these oscillations. Moreover, CNES RL04 solution, which is not based on AOD1B RL06, provides results similar to those shown from CSR RL06 and ITSG 2018 data. The CNES RL04 dealiasing product uses the ERA-Interim (ECMWF Reanalysis from January 1989 onward) model for the atmospheric part, and the TUGO (Toulouse Unstructured Grid Ocean) barotropic model for the oceanic response to the ERA-interim pressure and wind forcing [51]. Nevertheless, the work of [37] showed that the update of dealiasing model in RL06 solutions played a minor role in GSMAM improving, while the change of mean pole model from cubic to linear had a great influence. Other background models, including mean gravity field models, ocean tides models or planetary ephemerides for perturbations induced by the large planets of the solar system, vary from one GRACE solution to another.
It should be kept in mind that apart from background models, GRACE solutions developed by various data centers differ from each other in terms of data processing algorithms and various parametrization schemes of accelerometers and K-band ranging, whose impact on PM excitation determination is hard to quantify. Moreover, not all details on the processing of GRACE observations have been made available to users by the data centers. Therefore, it is difficult to determine which processing algorithm or model used most contributed to the discrepancies between different GSMAM.
It is worth mentioning that, although the GRACE mission has come to an end, the leading GRACE data centers are still working on new RL07 solution that will cover the whole period of the mission duration. The data processing methods of new series will be consistent with those used to process observations of the successor of the GRACE mission, GRACE Follow-On, launched in May 2018. Therefore, future research on the topic and evaluation of new GRACE RL07 solutions should be conducted. Recently, in December 2019, the specialists from CNES/GRGS have developed their new release of GRACE solution, namely CNES/GRGS RL05, available at monthly and 10-day intervals. Moreover, there are other available types of GRACE data, such as terrestrial water storage (TWS) gridded data (GRACE Level-3 data) or local mass concentration (mascon) solutions [64] which are not considered here.
Moving back to the GRACE solutions provided by CNES/GRGS, it should be mentioned that, apart from the data from GRACE, they use also the data from SLR satellites, namely Starlette, Stella, Lageos, Lageos-2 and Ajisai [51]. These additional observations were exploited in order to improve the quality of the coefficients of very low degrees of the gravity field. This had an impact on improved agreement between GAO and GSMAM from CNES RL04. Indeed, compared to other GRACE solutions, CNES RL04 provides better results in non-seasonal and 450-1000-day spectral bands. However, GSMAM series from previous release of the solution, CNES RL03, was characterized by poor consistency with GAO [15].
The analyses of GSMAM trends presented in Sections 3.2 and 3.3 reveal that, although the update of the mean pole model in GRACE RL06 data resulted in increased agreement between GSMAM and GAO trends, the full consistency is not achieved as GSMAM produced higher trends or both χ 1 and χ 2 . Moreover, these characteristics depends on time period considered. Apart from the adopted mean pole model, applying the GIA model might also have an impact on the size of the trends in the PM excitation series obtained from GRACE. However, it should be kept in mind that the post-glacial rebound is not modelled in GRACE Level-2 data [31,33,65,66], with the exception of the CNES solutions for which GIA is taken into account by the a priori gravity field model [51,67]. The GIA model is included in GRACE Level-3 products (available via the GRACE Tellus website https://grace.jpl.nasa.gov/), which represent monthly mass changes given in grids. The former studies by our group [15] showed visible differences between trends of excitation series computed from GRACE Level-2 and Level-3 data. We found that the trends in the latter were noticeably smaller. However, in this paper we focus only on hydrological excitations computed directly from GRACE Level-2 GSM coefficients. As GIA is a long-term linear process, it does not influence the results for the period bands considered in this study.
During the analysis of non-seasonal time series conducted in the previous section, we detect a visible quasi-5-year oscillation in the χ 2 ( Figure 10) which is also noticeable in 1000-3000-day spectral band (Figure 11b). This oscillation might show the effects from groundwater change that was previously shown to affect mostly inter-annual PM excitation variability [68,69]. In fact, noticeable groundwater depletion is currently observed in many major aquifers [70]. As indicated in other research, the main driver of long-term and non-seasonal PM variation is also ice mass loss (e.g., 4). Since the beginning of the 21st century, ice mass loss has been observed in the Greenland and Antarctic ice sheets, mainly as a result of the warming climate [2,71]. However, the quasi-5-year oscillation is also noticeable for HAM + SLAM series, but LSDM, which is used to obtain HAM + SLAM series, includes only shallow groundwater layer and does not model the long-term ice mass changes in Greenland and Antarctica. Therefore, we expect that there are other drivers for this variation in PM excitation. The five-year change might be related to a change in continental water distribution due to different precipitation patterns. During 2007 and 2013 the trend in continental water distribution was different from the years before and after. It is not obvious, however, whether the variations in land water patterns relate to overall water mass (distribution between land and ocean) or a land-internal redistribution. Nevertheless, the signal is more pronounced in χ 2 than in χ 1 , which more clearly indicates the direction of water exchange between land and ocean. Some authors suggested that electromagnetic core-mantle coupling also affects PM [72], but on decadal time scales, causing, for example, a Markowitz wobble with a period about 30 years [73].
It should be emphasized that for all considered oscillations, HAM + SLAM provides visibly higher consistency with GAO in χ 2 than in χ 1 . While for χ 1 HAM + SLAM reveals rather poor consistency with GAO compared to results from GRACE solutions, for χ 2 this agreement is similar or even higher than for GSMAM. It should be kept in mind that χ 2 might be more related to redistribution on the northern continents and between land and ocean, whereas χ 1 could be more sensitive to ice or glacier mass changes. Therefore, it is unsurprising that a higher agreement between considered series is found in the case of χ 2 . We can therefore conclude that HAM + SLAM is able to model observed signals in PM at longer periods, which can provide an indication of liquid water and snow redistributions; however, it cannot model ice and glacier mass change. In fact, LSDM includes only the annual snow accumulation and melting on top of the glaciated regions.

Conclusions
In this paper, we evaluate GSMAM from the newest GRACE RL06 series of global gravity field models by comparing them with GAO derived from precise measurements of the pole coordinates. We also consider the combined GRACE + SLR solution as well as HAM + SLAM series computed from numerical models processed at GFZ. The trends, seasonal and non-seasonal changes, and oscillations in a number of different spectral bands are considered. GAO, GSMAM, GRACE + SLR and HAM + SLAM series are also analyzed at different time intervals, which are determined on the basis of different quality of GRACE data during these periods.
The study finds that the variability of both GAO and GSMAM as well as trends of the series are variable for different time intervals. Moreover, the consistency between GSMAM and GAO also varies between periods and is the highest in 2006-2014. GSMAM trends are consistent with each other but they overestimate the trends in GAO for both χ 1 and χ 2 . The analyses show that for seasonal PM variations, it is difficult to identify the GRACE solutions that correspond best with GAO. For non-seasonal PM variations, the situation is clearer: the highest consistency with GAO are detected for the CNES RL04 solution.
We also identify the GRACE series that corresponds best with GAO in the chosen spectral bands and the results are as follows: (1) For 1000-3000-day oscillations, the highest correlation and variance agreement with GAO are observed for GSMAM from ITSG 2018 and CSR RL06 solutions. (2) For 450-1000-day changes, CNES RL04 is shown to provide the highest consistency between GSMAM and GAO. (3) In the 100-450-day spectral band, GSMAM computed from CSR RL06, ITSG 2018, and CNES RL04 solutions agree best with GAO. (4) For periods of of 60-100 days, the best correlation and variance agreement with GAO are observed for the ITSG 2018, CNES RL04, and CSR RL06 solutions. The level of agreement between GSMAM and GAO depend on the frequency band considered and is the highest for decadal changes with periods of 1000-3000 days and the lowest for short periods (60-100 days).
The results of this study demonstrate that there is no single solution that could give the highest consistency with GAO in all considered spectral bands. Nevertheless, the CSR RL06, ITSG 2018 and CNES RL04 appear to provide the best agreement with GAO for most of the oscillations investigated. However, it should be kept in mind that the latter use additional data from SLR satellites.
This study shows that the combination of GRACE and SLR, which is recommended during the use of C 21 , S 21 coefficients from GFZ RL06 data, improves the consistency between GSMAM and GAO. The formal errors of GRACE + SLR solution are smaller than errors of any GSMAM series. The combined series are also less noisy than GFZ RL06, especially in the last period of GRACE satellites activity. However, despite the fact that GRACE + SLR provides clearly better consistency with GAO than GSMAM from GFZ RL06, it does not perform as well as CSR RL06. GRACE + SLR results represents rather the mean of results obtained from considered GRACE solutions.
We note that the HAM + SLAM series from the numerical models processed at GFZ is unable to provide such high correlations with GAO as the GSMAM data, especially in χ 1 . Similar findings have been shown in previous studies by our group, where different hydrological and climate models were compared with gravimetric measurements and geodetic excitations [14,15,17,30]. LSDM-based hydrological excitations are found to be affected by weaknesses especially in χ 1 of 450-1000 and 60-100-day spectral bands. However, [15] showed that LSDM gives more satisfactory results than other hydrological and climate models. This study shows that the results obtained from HAM + SLAM are comparable with those provided by GRACE solutions at the 100-450-day and seasonal bands. Moreover, if only χ 2 would be considered, it will be found that HAM + SLAM provides one of the best results for all oscillations and time intervals considered.
Finally, it should be emphasized that none of the considered data sets are free of errors. The GAO series, which is taken as our reference for the quality assessment of GRACE-and model-based PM excitations, is affected by uncertainties of atmospheric and ocean models, while GRACE data are influenced by the errors in background models. Another method to assess uncertainties in hydrological excitations from different GRACE solutions, might be the estimating the noise in the series, using the three-cornered hat method, for example [74,75]. However, this could provide only internal agreement between different GRACE series since no external data are used. Nevertheless, such an assessment would be free of errors in atmosphere and ocean models used for GAO computation. Hopefully, the successor of the GRACE mission, GRACE Follow-On with the innovative Laser Ranging Interferometer to be used for the satellite-to-satellite tracking, will provide determinations of temporal Earth's gravity field changes with a higher quality than available so far, which also could improve the estimation of gravimetric excitation from satellite missions.