Determining and Evaluating the Hydrological Signal in Polar Motion Excitation from Gravity Field Models Obtained from Kinematic Orbits of LEO Satellites

: This study evaluates the gravity ﬁeld solutions based on high-low satellite-to-satellite tracking (hl-SST) of low-Earth-orbit (LEO) satellites: GRACE, Swarm, TerraSAR-X, TanDEM-X, MetOp-A, MetOp-B, and Jason 2, by converting them into hydrological polar motion excitation functions (or hydrological angular momentum (HAM)). The resulting HAM series are compared with the residuals of observed polar motion excitation (geodetic residuals, GAO) derived from precise geodetic measurements, and the HAM obtained from the GRACE ITSG 2018 solution. The ﬁndings indicate a large impact of orbital altitude and inclination on the accuracy of derived HAM. The HAM series obtained from Swarm data are found to be the most consistent with GAO. Visible di ﬀ erences are found in HAM obtained from GRACE and Swarm orbits and provided by di ﬀ erent processing centres. The main reasons for such di ﬀ erences are likely to be di ﬀ erent processing approaches and background models. The ﬁndings of this study provide important information on alternative data sets that may be used to provide continuous polar motion excitation observations, of which the Swarm solution provided by the Astronomical Institute, Czech Academy of Sciences, is the most accurate. However, further analysis is needed to determine which processing algorithms are most appropriate to obtain the best correspondence with GAO.


Introduction
The accurate determination of spatial and temporal changes of the Earth's gravity field is essential in many applications, including solid Earth science, hydrology, oceanography, glaciology, and geodesy. The gravity field varies in space and time due to disturbances in mass redistribution of Earth's surficial fluids, which include the atmosphere, oceans, and the land hydrosphere. These variations cause changes in Earth orientation parameters (EOPs), which describe the rotation of the planet. These parameters are: precession/nutation, polar motion (PM) and length-of-day (LOD) variations. EOPs are essential for a number of applications, including precise positioning and navigation in space and on the Earth's surface, pointing of astronomic instruments, and communication with deep space objects. Such great importance of these parameters results from the fact that they are necessary variables for transforming coordinates between a terrestrial reference frame (in which coordinates of ground stations are commonly available and used in all surveying tasks) and a celestial reference frame (in which the coordinates of ground stations are determined from space geodesy techniques such as Global Navigation Satellite Systems (GNSS)) because they provide the rotation between these two frames as a function of time.
Recent studies [10,12,13] have evaluated a few different hl-SST gravity field solutions obtained from orbits of Swarm satellites, using GRACE ll-SST observations for this purpose. The results from different institutes vary substantially, mainly due to different maximum degrees and orders of spherical harmonic representation of the geopotential and differences in gravity field estimation approaches and background models [12,13]. However, it has been stated that Swarm data could be able to cover the data gap between GRACE and GRACE-FO missions with satisfactory accuracy [9][10][11][12]. Moreover, some attempts to combine different Swarm solutions [12,13] have produced promising results and improved accuracy has been obtained. However, although the orbits of Swarm satellites have similar height and inclination as GRACE twin-satellites, they do not provide such high-quality data. The authors of Reference [13] showed that the satisfactory agreement between mass changes based on Swarm and GRACE data could be obtained when the coefficients of spherical harmonics representation of Earth's gravity field up to degree and order 15 is used. They demonstrated that the correlations between GRACE and Swarm estimations decrease sharply above degree 10, and the errors in Swarm estimates increase visibly above degree 15. They also showed that the higher consistency between GRACE and Swarm mass variations could be obtained in the periods of low solar activity. Indeed, due to their orbit inclination and altitude, the material from which they are built, and the characteristics of instruments aboard, Swarm satellites may be more sensitive to the solar activity. As a result, they produce noisier solutions than the GRACE-based ones.
The authors of Reference [5] evaluated monthly gravity field solutions based on kinematic orbits of GRACE A & B, GOCE [14], Swarm A, B & C [8], TanDEM-X [15], TerraSAR-X [16], MetOp-A, MetOp-B [17], and COSMIC [18] satellites and compared them with the results obtained from the GRACE CSR RL05 (Center for Space Research Release 5) solution [19]. The comparison of degree variances of the solutions based on observations from different satellite missions confirms they could all be used to estimate Earth's gravity field, but with different accuracies. The COSMIC, MetOp-A, and MetOp-B data provide results clearly inferior to the solutions obtained from other missions. Reference [5] suggests that this is mainly due to higher orbital altitude and worse GNSS positioning accuracy. For low degrees of geopotential, the TerraSAR-X, TanDEM-X, and Swarm satellites give comparable results to those obtained from GRACE and GOCE gravimetry satellite orbits. However, in the case of the Swarm mission, although two of the three satellites are orbiting at a lower altitude (430 km), the results are promising only if the data from all three satellites are combined [5].
Other studies have also shown that hl-SST can be useful for deriving long-wavelength gravity signals, related to several circumstances, e.g., mass changes over large river basins and ice sheets or glacial isostatic adjustment (GIA) signals [20]. References [21] and [22] show that combined hl-SST models based on either dedicated plus non-dedicated or only non-dedicated satellite data are able to estimate large-scale mass variations, annual signals, and trends. However, the solutions are limited because of a spatial resolution no better than 750 km and a high noise level [23]. Moreover, to achieve the best possible accuracy, several conditions should be taken into account during selection of the satellites, i.e., low altitude of the orbit (below 600 km), to provide an appropriate sensitivity of the satellite for gravity field variations, high inclination to avoid a polar gap in data, and good quality of GNSS observations. Some improvement in the accuracy of derived gravity field variations due to mass changes can be obtained through combining hl-SST data with satellite laser ranging (SLR) measurements [20,23].
In this study, we examined the usefulness of several global temporal gravity field models determined from kinematic orbits of selected LEO satellites (hl-SST technique) in designating PM variations. The variations of PM due to mass redistribution changes in the Earth's surficial fluids (atmosphere, ocean, and land hydrosphere) can be expressed as χ 1 and χ 2 components of the PM excitation function (atmospheric angular momentum (AAM), oceanic angular momentum (OAM), and hydrological angular momentum (HAM), respectively). The temporal variations of the χ 1 and χ 2 components are proportional to the changes of C 21 and S 21 (degree-2 order-one) coefficients of geopotential, respectively [24]. This relationship allows us to use global temporal gravity field models Remote Sens. 2019, 11, 1784 4 of 19 to quantify mass-related PM excitation. The GRACE monthly solutions, based on precise range-rate measurements between two satellites (ll-SST), have been widely used for this purpose following extensive research in recent years [25][26][27][28][29][30][31]. Special interest has been placed on HAM, as AAM and OAM are well established [32][33][34][35][36][37][38].
The assessment of gravity field models based on hl-SST measurements presented in this study was performed to find the most appropriate solution to bridge the gap between GRACE and GRACE-FO. Similar estimations have been performed by other authors [9,10,12]. However, all of these previous analyses are based mainly on comparison of degree variances of the solutions or regional water mass changes over ocean, land, and ice sheets. Such types of analyses require full spherical harmonic representation of the geopotential. Here, we focused our attention on PM variations that are described by the C 21 and S 21 coefficients only. However, it should be kept in mind that all five EOPs are equally important for transforming station coordinates between the celestial and terrestrial reference frames. Nevertheless, the mass variations of land hydrosphere, of which impacts on the Earth's rotation are analysed here, have a smaller impact on precession/nutation and LOD. Moreover, precession and nutation are well described by theoretical models.
In this study, we compared the hydrological part of mass-related PM excitation functions, expressed as HAM and derived from a number of hl-SST gravity field models. We not only considered the time series of these excitations, but also decomposed them into linear trends, seasonal, and non-seasonal oscillations. The study of the accuracy of these estimates was based on comparisons with reference data, the ll-SST GRACE ITSG 2018 solution [39], and hydrological signals in observed PM derived from precise geodetic measurements (also called geodetic residuals, GAO). Quantitative analyses were conducted based on comparison of correlation coefficients with reference data and the standard deviation (STD) of the series. This study aimed to identify the solutions that, on the one hand, could best match geodetic and gravimetric observations of PM and, on the other hand, cover the time period of the gap between GRACE and GRACE-FO. Finding such solutions is essential for providing consistent, uninterrupted, and high-quality PM data, which are needed in many geodetic applications.

Reference Data
To evaluate HAM functions obtained from models based on kinematic orbits of LEO satellites, we compared them with two reference series: (i) the hydrological signal in geodetically observed PM excitation (GAO), and (ii) the HAM function obtained from the latest monthly GRACE gravity field model ITSG 2018.

Hydrological Signal in Observed Polar Motion
The observed geodetic PM excitation (geodetic angular momentum, GAM) was computed from precise pole coordinates provided by space geodesy techniques such as GNSS, SLR, and very long baseline interferometry (VLBI). The equatorial components (χ 1 and χ 2 ) of GAM can be calculated from x and y pole coordinates by solving Liouville's equation [40,41]. The coordinates of the Earth's pole are available in the daily combined C04 series of EOPs [42,43] provided by the International Earth Rotation and Reference System Service (IERS) (https://www.iers.org/). To obtain the hydrological signal in observed excitation, GAM was reduced for the effects of atmosphere and ocean, using AAM and OAM: The expression GAO−AAM−OAM is usually shortened by taking only the first letters of GAM, AAM, and OAM to obtain G-A-O or simply GAO. Such residual series, denoted as geodetic residuals or GAO, reflect not only the impact of the land hydrosphere, but also barystatic sea-level changes due to the inflow of water from land into the oceans and some solid Earth-related signals, including tectonic signals and GIA, on PM excitation. While GAM is obtained from geodetic measurements, AAM and OAM are usually computed using the atmospheric pressure and wind speed, and ocean bottom pressure and currents. These variables can be obtained from geophysical models of atmosphere and ocean.
The following data sets were used to compute GAO:

GRACE ITSG 2018 Solution
The GRACE GSM data have the form of monthly time series of spherical harmonic coefficients of the Earth gravity field, with a specified maximum degree and order (d/o) (typically 60, 90, or 120). The non-tidal short-term atmospheric and oceanic signals were removed from the series, and therefore they reflect the impact of mass effects from land hydrosphere, the barystatic sea-level changes, GIA, and tectonic signals resulting from large earthquakes.
Recently, the leading GRACE data centres, including GFZ in Potsdam, Germany; Center for Space Research (CSR) in Austin, USA; the Jet Propulsion Laboratory (JPL) in Pasadena, USA; and the ITSG of the Graz University of Technology, Austria, have provided the newest temporal gravity field solutions, with updated processing algorithms and background models [39,[48][49][50].
In this study, we used the newest GRACE ITSG 2018 solution provided by the ITSG [39]. The motivation for this choice was that most of gravity field models from kinematic orbits of LEO satellites used in this study are also provided by this institute, so the processing methods and background models used are similar. Moreover, our previous researches [51][52][53] show that compared with series from other institutes, ITSG 2018 provides the highest correlation, variance, and trend agreement with GAO. However, it should be kept in mind that the consistency between GRACE-based HAM functions and GAO depends on specific oscillation and time period considered [51][52][53]. In general, the HAM trends from GRACE are usually stronger than GAO trends, especially for χ 1 component. In terms of amplitudes of oscillations, for some oscillations, they are stronger for GAO while for other they are weaker. Based on numerous comparisons of different oscillations in HAM series computed from different GRACE solutions, we indicated that ITSG 2018 provides the highest consistency with GAO in the largest number of cases [51][52][53]. The ITSG 2018 series have quasi-monthly time resolution, maximum d/o 60 and were accessed from http://ftp.tugraz.at/outgoing/ITSG/GRACE/ITSG-Grace2018/monthly/.

Gravity Field Models Obtained from Kinematic Orbits of the LEO Satellites
In this study, we used several temporal models of Earth's gravity field, developed based on the hl-SST method which uses kinematic orbits of LEO satellites, and provided by scientists from ITSG [5,6]. Among these models, there are series that use orbits of either gravimetric satellites (dedicated to gravity field determination) or satellites of which the main goals are not monitoring the gravity field (non-dedicated satellites). The models are based on orbits of one or few satellites of the same type. The institute also provides solutions that combine satellite data of different types. Summary information about satellites used by ITSG to develop gravity field models, including satellite type, mission duration, and orbit altitude and inclination, is included in Table 1. The data are publicly available to all users via http://ftp.tugraz.at/outgoing/ITSG/tvgogo/gravityFieldModels/.
Based on the data length and the considered time period, which should cover both time gap between GRACE and GRACE-FO (to potentially fill this gap) and at least the part of an operational period of the GRACE mission (to enable validation of models using the reference GRACE ITSG 2018 model), we selected twelve models from ITSG, Graz University of Technology and two models from ASU CAS, Prague ( Table 2). All GRACE solutions considered here were based on a combination of both GRACE A and GRACE B orbits, while all Swarm solutions were developed using a combination of Swarm A, Swarm B, and Swarm C orbits. We also considered two combined solutions, Combined v2 and Combined v3, which used data from several types of satellites, as listed in Table 3. The rest of the solutions considered were developed based on observations from a single satellite. To distinguish between GRACE-based and Swarm-based models from Graz and Prague, the former and the latter will be further abbreviated with ITSG and CAS, respectively. The models from other satellites and combined models were provided by ITSG only.
All models had the same form as the reference ITSG 2018 series, namely time series of spherical harmonic coefficients with the monthly time resolution and maximum degree and order 60 (with one exception for GRACE A, B ITSG v2 for which the maximum d/o is equal to 100).
The (χ 1 , χ 2 ) components of HAM functions (both from models based on kinematic orbits of LEO satellites and GRACE ITSG 2018 reference solution) were computed from (C 21 , S 21 ) coefficients of geopotential using the following formulas [24]: where ∆C 21 and ∆S 21 are the changes of coefficients of the Earth's gravity field; R e and M are the Earth's mean radius (6,378,136.6 m) and mass (5.9737 × 10 24 kg), respectively; A = 8.0101 × 10 37 kg·m 2 , B = 8.0103 × 10 37 kg·m 2 , and C = 8.0365 × 10 37 kg·m 2 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 Reference [24]).  Table 3. Satellites used by ITSG to develop combined models of Earth gravity field.

Time Series and Trends
Initially, we performed a comparative analysis of the HAM and GAO series, without considering specific oscillations. Because of the different time resolutions of analysed data (3 h for AAM and OAM, 24 h for GAM, and monthly for all gravity field models), we downsampled all series to monthly changes using a Gaussian filter with a full width at half maximum (FWHM) equal to 60 days.
The STD of the χ 1 and χ 2 components of GAO and different HAM functions (after removing linear trends) are provided in Table 4. The data in Table 4 were supplemented by the time period of each solution and its length. Notably, there were differences in data length between the considered solutions. However, all missions that are not dedicated to monitoring the gravity field are still operational (see Table 1), and therefore, the data from these missions may be useful in filling the gap between GRACE and GRACE-FO. Table 4 indicates that the STDs of χ 1 and χ 2 PM excitation functions from reference GAO and ITSG 2018 data were overestimated by values from all hl-SST estimates. Unsurprisingly, the most compatible with the reference series were GRACE AB (the best from CAS), Combined v2, and Combined v3 solutions. It should be noted that both combined solutions included orbits of dedicated satellites, which greatly contributed to the high consistency with GAO and ITSG 2018. However, the STD consistent with those obtained for reference PM excitation functions was also detected for the Swarm solution from CAS (Swarm ABC CAS). Surprisingly, these good results were not found for the Swarm ABC ITSG solutions. The reason for these differences may be the fact that both institutes used different approaches during computation of their gravity field models, namely the acceleration approach [7] for CAS solution and the short-arc approach [54] for ITSG solutions. The series computed from Jason 2, MetOp-A, and MetOp-B data revealed much higher STDs of χ 1 and χ 2 than reference data and other missions. This may result from visibly higher orbital altitude for these satellites (see Table 1). The Jason 2 solution was also affected by a lower inclination, which produced polar gaps in the data. Reference [5] shows that the lower accuracy of GNSS positioning for Jason 2, MetOp-A, and MetOp-B satellites could also contribute. In general, the higher the orbit of the satellite, the higher the STD of the obtained PM excitation function, and this could result from the fact that solutions from these higher orbits were noisier than series based on data from lower orbits.
The linear trends of χ 1 and χ 2 components of GAO and different HAM functions are shown in Table 5. Because the considered series were available in different time periods and in different lengths, to make them comparable with the reference series, the trends in GAO and HAM from ITSG 2018 were computed for the full period of their availability and also for the time period of each evaluated solution. It was noticeable that the trends of reference series were strongly dependent on the considered time period. For GAO, the trends varied from +1.83 to +5.62 milliarcseconds mas/year for χ 1 and from -0.60 to +3.39 mas/year for χ 2 . The trends of HAM functions computed from ITSG 2018 solutions had values between +5.76 and +9.74 mas/year for χ 1 and between −5.49 and +1.25 mas/year for χ 2 . These non-negligible differences suggested that GAO and ITSG 2018 excitation functions were characterized by non-linear trends. Table 5 shows that the best compatibility with χ 1 trends of the reference series can be provided by GRACE AB CAS, GRACE AB ITSG v2, Combined v2, Swarm ABC ITSG v2, and TerraSAR-X solutions. For χ 2 , this correspondence was noticeable for GRACE AB ITSG v3, Combined v2, Swarm ABC CAS, and Swarm ABC ITSG v3. The results showed that the Jason 2, MetOp-A and MetOp-B solutions, which were distinguished by the highest STDs in χ 1 and χ 2 , also visibly overestimated the trends in GAO and ITSG 2018 HAM series. The results from GRACE and Swarm orbits gave clearly superior results to the data from other satellites. However, TerraSAR-X provided satisfactory results for the χ 1 trend. It should be also emphasized that for χ 2 there were some periods where trends for GAO and ITSG 2018 HAM differed from each other in terms of their signs, such as for 2003. 12-2012.98, 2011.87-2016.47, 2013.87-2015.35, and 2013.12-2016.49. For the same periods, for χ 2 , HAM function from ITSG 2018 gave a visibly higher trend than GAO. Table 4. Standard deviations (STDs) of χ 1 and χ 2 components of geodetic residuals (GAO) and different hydrological angular momentum (HAM) series. Note that series were available in different time periods (given in decimal years) and with different period lengths (given in years).

Data
Time

Seasonal Variations
We now extended our assessment of HAM series into seasonal oscillations. To calculate these seasonal changes, we first removed linear trends from the time series. Then, the seasonal variations were computed using the least squares method [33]. The fitted model included a sum of sinusoids with the periods of 365.25, 182.625, and 121.75 days.
The time series of seasonal changes (sum of annual, semi-annual, and ter-annual oscillations) in GAO HAM are shown in Figure 1 (for all GRACE and combined solutions) and in Figure 2 (for Swarm, TerraSAR-X, and TanDEM-X solutions). The MetOp-A, MetOp-B, and Jason 2 solutions were excluded from further analyses because of the poor STD agreement and trend inconsistency with reference data.
As shown in Figure 1, for χ 1 , only ITSG 2018 and GRACE AB ITSG v3 provided good phase agreement with GAO, while the best amplitude agreement was detected for Combined v3 solution. Almost all GRACE models, both from ll-SST and hl-SST estimates (GRACE AB CAS, GRACE AB ITSG v3, and ITSG 2018), notably underestimated seasonal variation in GAO. For χ 2 , very good phase agreement was found between all PM excitation series. The most consistent with GAO in terms of χ 2 amplitudes was the GRACE AB ITSG v2 solution. Notably, the newest solutions from ITSG (GRACE AB ITSG v3 and Combined v3) visibly overestimated amplitudes of GAO in χ 2 , whereas GRACE AB CAS underestimated amplitudes of GAO in χ 2 . In general, better phase agreement with GAO was obtained for χ 2 component and this was also observed in previous works, e.g., in References [26,[28][29][30]32,37,55]. This resulted from spatial distribution of land and ocean that determine χ 2 to be more sensitive to mass changes over land, and χ 1 to be more sensitive to mass changes over ocean, ice, and glaciers. Figure 2 shows that, apart from Swarm ABC CAS and Swarm ABC ITSG v2, all hl-SST solutions from non-dedicated satellites visibly overestimated seasonal amplitudes of GAO for both χ 1 and χ 2 components. Moreover, for TanDEM-X and TerraSAR-X, the semi-annual χ 2 changes were almost as strong as annual ones. For χ 1 , this was apparent for the Swarm ABC ITSG v3 data. For geodetic observations of PM excitation, however, the seasonal signal was notably stronger than semi-annual and ter-annual signals. In general, for both χ 1 and χ 2 , the best amplitude and phase agreement with GAO and ITSG 2018 provided Swarm ABC CAS. Surprisingly, the former Swarm ABC ITSG v2 solution provided better amplitude agreement with GAO than the newer Swarm ABC ITSG v3.
The correlation coefficients between seasonal changes in GAO and HAM are shown in Figure 3. Because of different data lengths of the different series, we considered three periods: 2003-2013 (for GRACE AB CAS, GRACE AB ITSG v2, and Combined v2); 2011-2016 (for GRACE AB ITSG v3 and Combined v3); and 2014-2016 (for Swarm ABC CAS, Swarm ABC ITSG v2, Swarm ABC ITSG v3, TerraSAR-X, and TanDEM-X). To compare the results with those obtained for the reference GRACE ITSG 2018 solution, we also computed correlations between HAM functions from ITSG 2018 and GAO for the three considered time periods.
Keeping in mind that even high correlation between two time series can be statistically insignificant, for an objective assessment of the correlations, we determined the critical value of correlation coefficients. For this purpose, we used the autocorrelation function and the number of independent points. The autocorrelation function describes how a time series correlates with itself over different timescales. In other words, this function shows how rapidly time series changes or how rapidly it "forgets" about its previous values [56,57]. Therefore, the autocorrelation function measures for what time lag the correlation between series and between the same series but shifted with a lag will be close to zero (for what time shift time series become decorrelated). The decorrelation time are usually assumed as: (1) time required for autocorrelation function drop to the first zero-crossing, (2) twice the time required for autocorrelation function drop to 1/2, (3) the time required for autocorrelation function drop to 1/e, and (4) twice the time required for autocorrelation function drop to 1/e [56,57]. The autocorrelation function is symmetric about zero. To compute the critical value of correlation coefficient for the assumed significance level, the number of degrees of freedom (or number of independent points) has to be computed. This number can be calculated by dividing the number of series points by the decorrelation time. Then, the critical value of correlation coefficient for the computed number of independent points and assumed significance level can be read from statistical tables.
In this work, we estimated the number of independent points required for the autocorrelation function to drop to 1/e [54]. Next, based on Student's t-test and a chosen significance level equal to 0.95, we determined the critical value of the correlation coefficient for three considered periods separately. For the period 2003-2013, we obtained 47 independent points and a critical value of correlation coefficient equal to 0.24; for 2011-2016, 25 independent points and a correlation coefficient equal to 0.34; and for 2014-2016, 12 independent points and a correlation coefficient equal to 0.50.
The correlations between seasonal GAO and seasonal HAM from ITSG 2018 did not depend on the time period considered (Figure 3). Similar to amplitude and phase agreement shown in Figures 1 and 2 and results shown in previous works [26,[28][29][30]32,37,55], visibly better correlations were obtained for χ 2 . This was also observed for almost all hl-SST solutions. Notably, all Swarm solutions and TanDEM-X data revealed better consistency with GAO for χ 1 than for χ 2 . It was not surprising that the best correlation agreement with GAO was provided by the ll-SST solution from ITSG 2018. However, GRACE AB ITSG v3 and Swarm ABC CAS data gave correlations of a similar size. Nevertheless, taking into consideration critical values of correlation coefficient for different time periods, for χ 2 , HAM from Swarm ABC CAS was almost insignificant. The Swarm models provided by ITSG were apparently inferior to the Swarm ABC CAS in terms of their correlations with observed PM excitation. We also noted that combined and Swarm solutions from ITSG did not improve in the latest releases (v3) compared with former ones (v2). This may be caused by the different time periods considered. Newer combined solutions were also shorter and were obtained from orbits of slightly different satellites than the previous combined v2 model (see Table 3). The decrease of accuracy in Swarm ABC ITSG v3 and GRACE AB ITSG v3 was especially visible for χ 1 . Regarding HAM computed from SAR satellites, both TerraSAR-X and TanDEM-X solutions did not provide satisfactory correlations with GAO, and only χ 2 from the TerraSAR-X was well correlated with χ 2 from GAO. Taking into consideration all of the above findings, the Swarm ABC CAS model was the best candidate to provide seasonal PM excitation data for the period lacking observational data from GRACE/GRACE-FO.
periods, for χ2, HAM from Swarm ABC CAS was almost insignificant. The Swarm models provided by ITSG were apparently inferior to the Swarm ABC CAS in terms of their correlations with observed PM excitation. We also noted that combined and Swarm solutions from ITSG did not improve in the latest releases (v3) compared with former ones (v2). This may be caused by the different time periods considered. Newer combined solutions were also shorter and were obtained from orbits of slightly different satellites than the previous combined v2 model (see Table 3). The decrease of accuracy in Swarm ABC ITSG v3 and GRACE AB ITSG v3 was especially visible for χ1. Regarding HAM computed from SAR satellites, both TerraSAR-X and TanDEM-X solutions did not provide satisfactory correlations with GAO, and only χ2 from the TerraSAR-X was well correlated with χ2 from GAO. Taking into consideration all of the above findings, the Swarm ABC CAS model was the best candidate to provide seasonal PM excitation data for the period lacking observational data from GRACE/GRACE-FO.

Non-Seasonal Variations
This section addresses the analysis of non-seasonal variations in GAO and HAM, which were obtained after removing linear trends and seasonal oscillations from the time series. The series of non-seasonal changes in GAO and HAM are shown in Figure 4 (for all GRACE and combined solutions) and in Figure 5 (for Swarm, TerraSAR-X, and TanDEM-X solutions).

Non-Seasonal Variations
This section addresses the analysis of non-seasonal variations in GAO and HAM, which were obtained after removing linear trends and seasonal oscillations from the time series. The series of non-seasonal changes in GAO and HAM are shown in Figure 4 (for all GRACE and combined solutions) and in Figure 5 (for Swarm, TerraSAR-X, and TanDEM-X solutions).       Figure 4 shows that there is a very good consistency in terms of χ 1 and χ 2 amplitudes and phases between all hl-SST GRACE and combined solutions. However, all of them visibly overestimated amplitudes of reference PM excitation functions. Similar to the seasonal changes, HAM computed using GRACE AB CAS was distinguished by smaller amplitudes than the HAM obtained from the remaining hl-SST models. There were some periods in which evaluated models agreed in phase with the reference series, such as the same maximum in 2009 and minimum in 2005-2006 detected for χ 1 , or the same minimum in 2007 observed for χ 2 . Nevertheless, the considerable differences in the magnitude of amplitudes made this assessment difficult. Therefore, it is difficult to indicate the hl-SST solution that gives the best results. Even higher amplitudes were obtained for excitation functions computed using Swarm, TerraSAR-X and TanDEM-X data ( Figure 5). Among all models based on orbits of non-dedicated satellites, only Swarm ABC CAS provided good amplitude and phase agreement with the reference series.
The correlation coefficients between non-seasonal changes in GAO and HAM are shown in Figure 6. performed higher correlation consistency with GAO than ITSG 2018. However, it should be kept in mind that HAM functions from TanDEM-X data, although they provide high correlation coefficients with reference series, provide very big residual Root Mean Square RMS compared to GAO and HAM from ITSG 2018. Notably, HAM functions from GRACE, Swarm, and combined solutions provided by ITSG revealed decreased correlation agreement with GAO in the v3 release compared with in the former v2 release, and this might be an effect of different time periods and data lengths. Taking into consideration all of the above findings, the best candidate to provide non-seasonal PM excitation data for the period of lack of observational data from GRACE/GRACE-FO is the Swarm ABC CAS model.

Summary and Conclusions
In this study, we evaluated HAM computed using gravity field models determined from hl-SST of the orbits of LEO satellites (GRACE, Swarm, TerraSAR-X, TanDEM-X, MetOp-A, MetOp-B, and Jason 2) and some combined solutions. The assessment was based on the comparison with GAO derived from precise measurements of the pole coordinates and with the latest GRACE ITSG 2018 solution. The STD of the HAM series, trends, and seasonal and non-seasonal changes were considered [58].
The results indicated that HAM excitation functions obtained from hl-SST solutions visibly overestimated the magnitude of amplitudes and STD of the reference series. In general, use of a satellite with a higher orbital altitude led to a higher STD of the HAM series obtained from its orbit. However, an objective evaluation of hl-SST solutions was difficult because of the short period of common data.
The comparison of trends revealed a considerable dependence of trend values in the reference data from the period considered, which suggested a non-linear character of trends in PM excitation. Apart from hl-SST solutions obtained from GRACE orbits and from combined data, the most compatible with GAO and HAM ITSG 2018 in terms of trends were series computed from Swarm ABC ITSG v2 and TerraSAR-X for χ1 and Swarm ABC CAS and Swarm ABC ITSG v3 for χ2. For HAM series computed using MetOp-A, MetOp-B, and Jason 2 satellite data, we reported visibly stronger amplitudes than for other HAM series, as well as trends that did not match the trends in reference data. The reasons for this were probably due to the higher orbital altitude and less accurate GNSS data than for other considered missions.
Among the HAM series computed using solutions obtained from satellites not dedicated to

Summary and Conclusions
In this study, we evaluated HAM computed using gravity field models determined from hl-SST of the orbits of LEO satellites (GRACE, Swarm, TerraSAR-X, TanDEM-X, MetOp-A, MetOp-B, and Jason 2) and some combined solutions. The assessment was based on the comparison with GAO derived from precise measurements of the pole coordinates and with the latest GRACE ITSG 2018 solution. The STD of the HAM series, trends, and seasonal and non-seasonal changes were considered [58].
The results indicated that HAM excitation functions obtained from hl-SST solutions visibly overestimated the magnitude of amplitudes and STD of the reference series. In general, use of a satellite with a higher orbital altitude led to a higher STD of the HAM series obtained from its orbit. However, an objective evaluation of hl-SST solutions was difficult because of the short period of common data.
The comparison of trends revealed a considerable dependence of trend values in the reference data from the period considered, which suggested a non-linear character of trends in PM excitation. Apart from hl-SST solutions obtained from GRACE orbits and from combined data, the most compatible with GAO and HAM ITSG 2018 in terms of trends were series computed from Swarm ABC ITSG v2 and TerraSAR-X for χ 1 and Swarm ABC CAS and Swarm ABC ITSG v3 for χ 2 . For HAM series computed using MetOp-A, MetOp-B, and Jason 2 satellite data, we reported visibly stronger amplitudes than for other HAM series, as well as trends that did not match the trends in reference data. The reasons for this were probably due to the higher orbital altitude and less accurate GNSS data than for other considered missions.
Among the HAM series computed using solutions obtained from satellites not dedicated to gravity field monitoring, only HAM Swarm ABC CAS proved to be sufficiently consistent with seasonal signal in reference data. For non-seasonal change, the χ 1 and χ 2 correlations with reference data were found to be more consistent with each other than in the seasonal spectral band. For non-seasonal oscillations, Swarm ABC CAS and TanDEM-X models were found to be the most appropriate to replace PM observations of the GRACE/GRACE-FO missions. However, it should be kept in mind that for the short time period considered here, the correlation coefficients, although high, are almost insignificant statistically. Nevertheless, an analysis of correlation coefficients provides a general overview that indicates that HAM series are superior to the others. Our results also showed that the level of correlation between non-seasonal changes in GAO and HAM obtained from the latest ITSG 2018 solution was dependent on the time period considered, and was the lowest for the last three years of the GRACE mission duration.
We noted that the Swarm ABC CAS solution revealed visibly better HAM consistency of seasonal and non-seasonal changes with GAO than Swarm models developed by ITSG (Swarm ABC ITSG v2 and Swarm ABC ITSG v3 solutions). The reason for such discrepancies might be differences in gravity field estimation approaches or different background models used. Similarly, seasonal and non-seasonal HAM changes obtained from GRACE AB CAS were better correlated with GAO than the corresponding HAM computed using the GRACE AB ITSG v2 model.
The comparison of HAM obtained from GRACE AB v2 and v3, Swarm ABC v2 and v3, and Combined v2 and v3 solutions provided by ITSG demonstrated that most of the newest v3 releases did not considerably improve the accuracy of the determined PM excitation compared with the former v2 series. Notably, for seasonal variations, Swarm and combined solutions revealed decreased accuracy, while for non-seasonal changes, only the χ 2 component in HAM from the Swarm data had improved correlation with GAO. We suggest that the main reason why HAM from GRACE v3 solution is not clearly superior to the HAM obtained from GRACE v2 solution was the different time periods-in contrast to GRACE v2, the GRACE v3 model was developed for the period including the last few years of the mission when the accuracy of the GRACE measurements was lower. This decreased quality of GRACE data might also contribute to the lower accuracy of Combined v3 solution compared with that of Combined v2, as GRACE orbits were used in creating these models. Moreover, the Combined v2 and Combined v3 models were developed based on the orbits of different satellites. In particular, the Combined v2 solution was obtained from orbits of four satellites solely dedicated to the gravity field monitoring (CHAMP, GRACE A, GRACE B, and GOCE), whereas for Combined v3, only three such satellites were used (CHAMP was excluded). However, both Swarm ABC ITSG v2 and Swarm ABC ITSG v3 solutions were determined for the same time period but there were visible differences between them; seasonal, non-seasonal, and seasonal plus non-seasonal HAM variations revealed much bigger amplitudes and higher STD for Swarm v3 data, and differences in the trends were also visible. We assumed that the causes of such discrepancies may be updates to the processing algorithms or background models.
Overall, this study revealed that, as expected, neither of the considered hl-SST models provided the agreement with GAO at the same level of accuracy as ll-SST GRACE solutions. However, the Swarm data gave the most satisfactory results in both the seasonal and non-seasonal spectral bands. The problem of inconsistency between Swarm solutions from different processing centres can be overcome by combination of the different Swarm series. Some attempts to do so have been made by References [12,13]. However, they focused on analyses of geoid height [12] and global mass changes over ocean, continents, and ice sheets [13]. It is likely that such a combination of Swarm data will also lead to the increased consistency with geodetic observations of PM excitation. Nevertheless, detailed analyses have not been conducted so far, and they will be a subject of our future work.
Some high-quality information on PM excitation can be also obtained from the combination of data from different LEO satellites. In fact, the HAM functions computed from the Combined v2 and Combined v3 solutions considered here gave high correlation coefficients with GAO and satisfactory amplitude and STD agreement. However, these good results were mainly induced by including data from gravimetric satellites (CHAMP, GRACE, and GOCE). On the basis of these results, we concluded that the hydrological part of PM excitation can be derived from gravity field models obtained from kinematic orbits of LEO satellites. However, a notable impact on the accuracy of determined HAM has orbital parameters such as altitude and inclination, as well as computation approaches and algorithms and accuracy of GNSS observations. We identified that the best candidate to fill a PM data gap between GRACE and GRACE-FO are Swarm models provided by ASU CAS.