1 Introduction

The Earth’s rotation variables, comprising polar motion (X and Y) and change of rotation rate (length of day (LOD)), have been investigated for over four decades (Munk 1960). The Earth’s rotational changes are caused by movement of the Earth’s mass in mobile layers, which include atmosphere, oceans and the land-based hydrosphere.

There are many geophysical models of these fluid layers (the atmosphere, oceans and continental hydrosphere). However, because these geophysical processes are very complex, these models still suffer from uncertainties in the process descriptions, parametrizations and forcing (Güntner 2008). Consequently, estimates of the excitation of the Earth’s rotation derived from geophysical models are also affected by uncertainties (Gross et al. 2003; Seitz and Schmidt 2005; Chen and Wilson 2005; Nastula et al. 2007, 2011; Brzeziński et al. 2009; Dobslaw et al. 2010).

The most significant influence on Earth rotation changes, including polar motion at seasonal timescales, is the atmospheric angular momentum (AAM), whose orientation involves motion (wind) and mass (pressure) terms (Barnes et al. 1983; Chao and Andrew 1991; Gross et al. 2003).

Previous studies of Earth rotation changes, based on different ocean models containing important parameters, such as ocean bottom pressure and currents, demonstrate that the oceans also play a major role in the seasonal excitation of the polar motion (Wahr 1983; Dickey et al. 1993; Ponte et al. 1998; Gross et al. 2003).

Research by Gross et al. (2003), Brzeziński et al. (2005) demonstrated that adding the non-tidal oceanic angular momentum (OAM) estimate to that of AAM improves the agreement of model simulations with the seasonal excitation calculated from the geodetic determinations of Earth orientation. Nevertheless, the sum of AAM and OAM does not entirely explain the observed variations of polar motion (Ponte et al. 1998; Nastula and Ponte 1999; Brzeziński et al. 2009). The remaining power might be provided by the hydrological angular momentum (HAM).

Many studies on the impact of land hydrology on the polar motion excitation have been carried out in recent years (Chen and Wilson 2005; Seoane et al. 2011; Nastula et al. 2011; Jin et al. 2012; Adhikari and Ivins 2016). However, the global land hydrology models used in such studies each tend to provide significantly different amplitudes and phases for polar motion excitation (Nastula et al. 2011; Chen and Wilson 2005). These models may not represent the complete hydrological variation (Chen and Wilson 2005) because terrestrial water storage (TWS) has not yet been adequately measured at the continental scale (Lettenmaier and Famiglietti 2006).

The determination of equivalent water thickness (EWT) from the Earth’s gravity field observations represents an indirect approach for estimating land hydrology. Such geophysical information can be obtained from the Gravity Recovery and Climate Experiment (GRACE), and to a lesser extent from the Laser Geodynamics Satellite (LAGEOS) missions as time series of the spherical harmonic coefficient of the Earth’s gravity fields. These observations are often presented as monthly mass grids, which contain terrestrial water storage anomalies (i.e., in aquifers and river basins) from time-variable gravity data relative to a temporal mean. The storage anomalies are given as EWT values. The EWT patterns from GRACE are subject to the removal of artifacts by de-striping, filter smoothing methods and are also subject to aliasing errors (Jin et al. 2010). Several centers, such as the GeoForschungsZentrum (GFZ), the Center for Space Research (CSR) and the Jet Propulsion Laboratory (JPL) have computed the coefficients of series of time-variable gravitation and of the adequate layer of water storage (Brzeziński et al. 2009; Chen and Wilson 2005; Seoane et al. 2009). In this study, we used the GRACE data of RL05 series as computed by the CSR.

Our investigations are focused on the influence of different land hydrosphere surface parameters, especially soil moisture, on polar motion excitation functions at seasonal and non-seasonal timescales. Here, these different variables are obtained from the US-NASA-based Global Land Data Assimilation System (GLDAS) data products. In addition, this model of the land hydrosphere is based on the following different Land Surface Models: Common Land Model (CLM), MOSAIC Model, National Centers for Environmental Prediction/ Oregon State University/ Air Force/ Hydrologic Research Lab (NOAH) Model and the Variable Infiltration Capacity (VIC) Model (Fang 2008). The GLDAS generates optimal fields of land surface states and fluxes by integrating satellite and ground-based observational data products, using advanced land surface modeling and data assimilation techniques (Rodell et al. 2004). The water storage is the sum of soil moisture, snow water equivalent and canopy surface water (not counting changes in groundwater below the depth defined by the model). The data are available at http://ldas.gsfc.nasa.gov/gldas/GLDASdownload.php.

As a comparative dataset, we also use estimates of the EWT determined by the CSR from the GRACE satellite observations. Previous studies show a detailed analysis of continental scale water storage changes using GRACE and output from a high-quality global land hydrological modeling system (e.g., Chen et al. 2009; Landerer and Swenson 2012; Zotov et al. 2015).

Table 1 Global average values of annual precipitation, evaporation, total runoff, depth-integrated soil moisture and accumulated snow from each GLDAS model

External transfers of angular momentum, such as from global geophysical fluids to the solid Earth, are described using the equatorial components \(\chi _1\) and \(\chi _2\) of the polar motion excitation functions, where \(\chi _1\) is positive toward \(0^{\circ }\)E longitude and \(\chi _2\) is positive toward longitude \(90^{\circ }\)E (Barnes et al. 1983). Observationally, the geodetic excitation functions of polar motion can be determined on the basis of the equations of motion using the measured xy components of the pole (Eubanks 1993). The coefficients of the second degree and of the first order of the Earth gravity field are proportional to variations of the equatorial component \(\chi _1\), \(\chi _2\) of the series of the hydrological excitation function of polar motion. Satellite gravimetry data, in the form of the GRACE RL05 equivalent water heights (EWH) from Center for Space Research were used. The mass effects of the ocean and atmosphere and postglacial rebound have been removed, so in this way hydrological excitation of polar motion can be estimated from the gravimetric data, where EWH comprises terrestrial water in all of its components—biomass, surface water, ice, snow, soil moisture and groundwater, and is called here gravimetric excitation function. This gravimetric function can be compared with the mass term of the geophysical excitations of polar motion.

Our analysis includes two steps: first, the determination and comparisons of regional patterns of hydrological excitation functions of polar motion; and second, comparison of the global hydrological excitation functions determined from the GLDAS and GRACE data with a hydrological signal in the geodetic excitation function of polar motion.

2 Data descriptions

2.1 GLDAS data

Global and regional hydrological excitation functions of polar motion are determined using the different selected GLDAS realizations: GLDAS CLM, GLDAS MOSAIC, GLDAS NOAH and GLDAS VIC. Each of these models contain different numbers of soil moisture layers, which are taken into account when the different hydrological excitation functions are determined (see Table 1).

The following GLDAS dataset contains a series of land surface parameters simulated from the Common Land Model (CLM) V2.0 model in the Global Land Data Assimilation System (GLDAS) (Rodell and Beaudoing 2007a), the MOSAIC model in the Global Land Data Assimilation System (GLDAS) (Rodell and Beaudoing 2007b), the NOAH 2.7.1 model in the Global Land Data Assimilation System (GLDAS) (Rodell and Beaudoing 2013), the Variable Infiltration Capacity (VIC) model in the Global Land Data Assimilation System (GLDAS) (Rodell and Beaudoing 2007c). The number of vertical levels for soil moisture content is specified in a different number of layers and at different depths in each model (see Table 1). All of the described GLDAS models are at 1\(^\circ \) resolution in latitude and longitude and range from 1979 to the present at a monthly temporal resolution. These simulations were forced by a combination of NOAA/GDAS atmospheric analysis fields, spatially and temporally disaggregated NOAA Climate Prediction Center Merged Analysis of Precipitation (CMAP) fields and observation-based downward shortwave and longwave radiation fields derived using the method of the Air Force Weather Agency’s AGRicultural METeorological modeling system (AGRMET). The simulations were initialized on January 1, 1979 using soil moisture and other state fields from the GLDAS model climatology for that day of the year.

2.2 EWT changes from GRACE data

In this study, the CSR-based GRACE RL05 (200 km Gaussian smoothing) solution over the land is used from April 2002 to December 2014 in terms of EWT given in cm. The GRACE land monthly data are available at http://grace.jpl.nasa.gov, supported by the NASA MEaSUREs Program (Swenson 2012; Landerer and Swenson 2012; Swenson and Wahr 2006). The standard products of the GRACE data are sets of spherical harmonic coefficients describing the monthly variations in Earth’s gravity field, which were converted to estimate changes in mass at the surface (Wahr et al. 1998). The spatial sampling resolution of GRACE EWT grids is also 1 degree, and a spatial Gaussian averaging filter with a half-width of 200 km has been applied.

2.3 Geodetic residuals

To estimate the role global hydrological and gravimetric excitation functions of polar motion in polar motion, the geodetic residuals GAM-AAM-OAM (GAO) are used in our comparisons. They are computed as the differences between the geodetic angular momentum (GAM) and the sum of AAM and OAM. Similar comparisons were considered in studies by Nastula et al. (2011), Jin et al. (2010, 2012). It should be emphasized here that geodetic residuals are affected by uncertainties both in atmospheric and oceanic models (e.g., Naito et al. 2000; Zhou et al. 2005). In our study, we use the time series obtained from the International Earth Rotation and Reference System Service (IERS) C04 series of polar motion (Bizouard and Gambis 2009). The atmospheric excitation functions (AAM) are derived from the time series of NCEP/ NCAR reanalysis data (Salstein et al. 1993). The oceanic excitation functions (OAM), including bottom pressure and currents, are computed based on the ECCO-JPL ocean model (Gross et al. 2003).

3 Methods

In this study, the hydrological and gravimetric excitation functions of polar motion are computed using the following formulae given by Eubanks (1993):

$$\begin{aligned} \left[ \begin{array}{cc} \chi _1\\ \chi _2 \end{array}\right]= & {} -\frac{1.098R^{2}_{e}}{C-A}\int \!\!\!\int {\varDelta }{q}(\phi , \lambda , t) \nonumber \\&\cdot \sin (\phi )\cos (\phi ) \left[ \begin{array}{cc} \cos (\lambda )\\ \sin (\lambda ) \end{array}\right] \mathrm{d}S, \end{aligned}$$
(1)

where \({\varDelta }q(\phi ,\, \lambda ,\, t)\) represent the changes in water storage per unit area, \(R_e\) is the Earth’s mean radius, \(\mathrm{d}S\) is the surface element area, and C and A are the Earth’s principal moment of inertia. The factor 1.098 accounts for the combined effects of the yielding of the solid Earth to the surface load and rotational deformations (Eubanks 1993).

Here, changes in terrestrial water storage \({\varDelta }q\) for hydrological excitation functions are calculated from each GLDAS model, from state (snow water equivalent, SNW; total soil moisture in all layers, SM) and flux variables (precipitation, P; evaporation, E, total surface runoff, R):

$$\begin{aligned} ({\varDelta }\mathrm{TWS}_1)_n= & {} (P_{n}-E_{n}-R_{n}), \end{aligned}$$
(2)
$$\begin{aligned} ({\varDelta }\mathrm{TWS}_2)_n= & {} (\mathrm{SM}_{n+1}+\mathrm{SNW}_{n+1})-(\mathrm{SM}_{n}+\mathrm{SNW}_{n}),\nonumber \\ \end{aligned}$$
(3)

Subscripts (n) and (\(n+1\)) represent a given month and the next month, respectively.

In this study, the global and regional hydrological excitation functions are determined from TWS changes containing all flux variables (as shown in Eq. 2) and are called HAM functions, and from TWS changes assumed to be the sum of state variables, mainly soil moisture in all layers and snow water equivalent, and are referred to as SM HAM functions (Eq. 3).

Fig. 1
figure 1

Maps of the variances of the complex-valued components \(\chi _1+{i}\chi _2\) of the regional HAM excitation functions of polar motion computed from different GLDAS models: GLDAS CLM, GLDAS NOAH, GLDAS MOSAIC, GLDAS VIC and from the GRACE CSR RL05 gravimetric data. The units are mas. The solution for all excitations are in \(1{^\circ } \times 1{^\circ }\) grids

The changes in the terrestrial water storage \({\varDelta }q(\phi ,\, \lambda ,\, \hbox {t})\) based on the gravity data are computed from the spherical harmonic coefficients of the Earth’s gravity field (Wahr et al. 1998):

$$\begin{aligned}&{\varDelta }\overline{\sigma }(\lambda , \phi ) = \frac{2\pi a\rho _{ave}}{3} \sum \limits _{l=0}^{L} \frac{2l+1}{1+k_l}W_l \sum \limits _{m=0}^{l} \overline{P}_{lm}\nonumber \\&\qquad \cdot (\sin \phi ) ({\varDelta }\overline{c}_{lm}\cos m\lambda +{\varDelta }\overline{s}_{lm}\sin m\lambda ) \end{aligned}$$
(4)

and

$$\begin{aligned} {\varDelta }q(\lambda , \phi )=\frac{{\varDelta }\overline{\sigma }(\lambda , \phi )}{\rho _w}. \end{aligned}$$
(5)

In Eqs. 4 and 5, a is the radius of the Earth, \(\overline{P}_{lm}\) is the normalized associated Legendre function, \(\rho _{ave}\), \(\rho _{w}\) are the average density of the Earth and water, respectively; \(\overline{c}_{lm}\) and \(\overline{s}_{lm}\) are dimensionless coefficients, and \(W_l\) is the Gaussian averaging kernel defined in Wahr et al. (1998).

All of these data series are available at a resolution of one month. The common time span for all considered data covered 2002.3 to 2014.9.

4 Results

4.1 Regional patterns

Our goal in this study is to indicate how those hydrological land-based GLDAS models differ from each other in terms of the determination of the global HAM and SM HAM excitation functions. First of all, global simulated average values of annual precipitation, evaporation, total runoff, depth-integrated soil moisture and accumulated snow are determined by each GLDAS model and are presented in Table 1. These annual average values, computed from 2003 to 2014, are different for each GLDAS model and generally may influence the determination of TWS values.

As a first step, we compare the maps of variations of these regional hydrological excitation functions with each other as well as with the regional gravimetric excitation function determined from the GRACE CSR RL05 datasets. Figure 1 compares the geographic variability of the amplitudes of complex-valued components \(\chi _1+i\chi _2\) of the hydrological and gravimetric excitation functions estimated with the different realizations of the GLDAS models, using TWS calculated as given in Eq. 2, and the GRACE data. Figure 2 also compares the geographic variability of the amplitudes of complex-valued components \(\chi _1+i\chi _2\) of the hydrological and gravimetric excitation functions estimated from different realizations of the GLDAS models, using TWS calculated as the sum of soil moisture in all layers and snow water equivalent, and the GRACE data. These results suggest that the areas with the largest variances are generally similar for the different land hydrology models and for the GRACE data. These areas are the west coast of North America, the Brazilian Highlands, tropical southern Africa, Europe, the Mississippi basin, the Great Lakes region of North America and the southern region of the Asian continent. These results are similar to those of earlier investigations by Nastula and Salstein (2012), Winska et al. (2016), in which other land hydrology models were taken into account.

Fig. 2
figure 2

Maps of the variances of the complex-valued components \(\chi _1+i\chi _2\) of the regional SM HAM excitation functions of polar motion computed from the different layers of soil moisture and accumulated GLDAS models: GLDAS CLM, GLDAS NOAH, GLDAS MOSAIC, GLDAS VIC and from the GRACE CSR RL05 gravimetric data. The units are mas. The solution for all excitations are in \(1{^\circ }\times 1{^\circ }\) grids

The largest amplitudes of the regional excitation functions, shown in Fig. 1, are the results of two models: GLDAS VIC and GLDAS MOSAIC. The smallest amplitude is obtained from GLDAS CLM model and the GRACE CSR RL05 gravimetric data. The GLDAS NOAH model is associated with intermediate amplitudes. The amplitudes of the pattern based on the GRACE data are most similar to those of the GLDAS CLM. However, the GRACE-based pattern is large over the Himalayas, while in the GLDAS CL-based pattern the variability is most distinct in eastern and western North America.

Similarly, the amplitudes of the regional excitation functions, shown in Fig. 2, are similar to those shown in Fig. 1. There is a slight decrease in amplitude from regional HAM and SM HAM functions from the GLDAS NOAH, GLDAS MOSAIC and GLDAS VIC models. However, the amplitudes of regional HAM and SM HAM functions from the GLDAS CLM model are significantly decreased.

4.2 Time series comparisons

In this section, we compare global variations of hydrological excitation functions of polar motion with the geodetic residuals (GAO) in terms of time series variability, correlation coefficients, power spectrum and phasor diagrams.

Global HAM and SM HAM excitation functions, the GAO and global gravimetric excitation functions of polar motion were smoothed using the Gaussian method. Linear trends, resulting from the solid Earth Post Glacial Rebounds effects, and from climate changes (Chen et al. 2013), are removed from all time series. Here, the Greenland mass variations are not considered in HAM, SM HAM and gravimetric excitation functions, since the accumulated snow is poorly estimated by GLDAS models (Wenbin and Rong 2013). To investigate how the time series differ from each other, we also conducted interpolation of each time series in the step of 30 days from 2002.3 to 2014.9. The differences between regional patterns of HAM are also reflected in the time series of the global hydrological excitation functions of polar motion shown in Figs. 3 and 4.

Fig. 3
figure 3

Comparison of the components, \(\chi _1\) and \(\chi _2\), of the GAO, defined as the difference in the geodetic excitation function and sum of the AAM and OAM excitation function of polar motion with the HAM excitation functions computed from the flux variables of GLDAS models: GLDAS CLM, GLDAS MOSAIC, GLDAS NOAH and GLDAS VIC. Each plot also contains components, \(\chi _1\) and \(\chi _2\), of the gravimetric excitation functions computed from the GRACE CSR RL05 data

Fig. 4
figure 4

Comparison of the components, \(\chi _1\) and \(\chi _2\), of the GAO, defined as the difference in the geodetic excitation function (GAM) and sum of the AAM and OAM excitation functions of polar motion with the SM HAM excitation functions computed from the different layers of soil moisture and accumulated snow in GLDAS models: GLDAS CLM, GLDAS MOSAIC, GLDAS NOAH and GLDAS VIC. Each plot also contains components, \(\chi _1\) and \(\chi _2\), of the gravimetric excitation functions computed from the GRACE CSR RL05 data

These results indicate that the four solutions for the different GLDAS models (NOAH, CLM, VIC, MOSAIC) do not strongly agree in both the \(\chi _1\) and \(\chi _2\) oscillations (Fig. 3). Similarly, the hydrological solutions from sum of all layers of soil moisture also do not strongly agree with the GAO and the GRACE data (Fig. 4). These HAM and SM HAM excitations are comparable to each other for both the \(\chi _1\) and \(\chi _2\) components, with the exception of the HAM and SM HAM functions computed from the GLDAS CLM model. The amplitudes of HAM functions are higher than the amplitudes for SM HAM functions and are also comparable to GRACE excitation and the GAO.

The time series determined for the four selected GLDAS models are dominated by the annual oscillations (Figs. 3, 4). The four solutions from the different GLDAS models (NOAH, CLM, VIC, MOSAIC) do not strongly agree with each other for both the \(\chi _1\) and the \(\chi _2\) oscillations. Likewise, the time series determined for the sum of soil moisture in all layers are also dominated by the annual oscillation (Fig. 4).

Table 2 Zero-lag correlation coefficients between the GAO, global HAM excitation functions and gravimetric excitation functions of polar motion, \(\chi _1\) and \(\chi _2\), the GAO were calculated as differences between GAM (C04 series) and the sum of AAM (the NCEP/NCAR model) and OAM (the ECCO model)
Table 3 Zero-lag correlation coefficients between the GAO, global SM HAM excitation functions and gravimetric excitation functions of polar motion, \(\chi _1\) and \(\chi _2\), calculated from different soil moisture layers of GLDAS models. The GAO were calculated as differences between GAM (C04 series) and the sum of AAM (the NCEP/NCAR model) and OAM (the ECCO model)

The zero-lag correlation coefficient between the GAO and the different HAM and SM HAM time series and gravitational excitation functions was also computed and is shown in Tables 2 and 3. Focusing only on the GLDAS models, it is clear that the best agreement exists between the GAO and the HAM GLDAS NOAH excitation functions in both the \(\chi _1\) and \(\chi _2\) components of polar motion and the correlation coefficient values are 0.50 and 0.35, respectively. An examination of the other time series indicates that the best agreement occurs between the GAO and the GRACE time series in the \(\chi _2\) component of polar motion (correlation coefficient \(=0.70\)). In the case of the \(\chi _1\) component, agreement between the GAO and the GRACE data is higher than the agreement between the GAO and the HAM excitation functions.

The correlations between the GAO and SM HAM functions are lower for the \(\chi _1\) components than for the HAM functions; however, the correlations are higher for the \(\chi _2\) components (Tables 2, 3). Here, a very good agreement is reached between the GAO and the GLDAS MOSAIC soil moisture excitation functions. The gravimetric excitation functions agree the best with the GLDAS NOAH functions and their agreement is better than for SM HAM functions.

Fig. 5
figure 5

Comparison of the non-seasonal components, \(\chi _1\) and \(\chi _2\), of the GAO, defined as the difference in the geodetic excitation function (GAM) and sum of the AAM and OAM excitation function of polar motion with the global HAM excitation functions computed from flux variables of GLDAS models: GLDAS CLM, GLDAS MOSAIC, GLDAS NOAH and GLDAS VIC. Each plot also contains components, \(\chi _1\) and \(\chi _2\), of the gravimetric excitation functions computed from the GRACE CSR RL05 data

Fig. 6
figure 6

Comparison of the non-seasonal components, \(\chi _1\) and \(\chi _2\), of the GAO, defined as the difference in the geodetic excitation function (GAM) and the sum of the AAM and OAM excitation function of polar motion with the SM HAM excitation functions computed from different layers of soil moisture and accumulated snow in GLDAS models: GLDAS CLM, GLDAS MOSAIC, GLDAS NOAH and GLDAS VIC. Each plot also contains components, \(\chi _1\) and \(\chi _2\), of the gravimetric excitation functions computed from the GRACE CSR RL05 data

Focusing only on the agreement between GLDAS HAM functions indicates that general compatibility between each of the GLDAS functions is high (correlation coefficients\(>0.7\), excepting GLDAM CLM model) and is better for \(\chi _2\) than for \(\chi _1\) components. The GLDAS CLM HAM and SM HAM functions do not agree well with another GLDAS excitations for the \(\chi _1\) components.

After removing the 365.25, 180.0 and 120.0 day oscillations from these time series and from the GAO, we also compared non-seasonal hydrological and gravimetric excitation functions with each other as well as with the GAO and the results are shown in Figs. 5 and 6.

The GLDAS NOAH, and MOSAIC HAM and SM HAM excitation functions agree well with the non-seasonal GRACE and the GAO excitation functions for the \(\chi _1\) component. The GLDAS NOAH HAM and GLDAS SM HAM estimates also agree remarkably well with the GRACE excitation functions as shown in Fig. 6.

The spectra of the complex-valued components of \(\chi _1+i\chi _2\) for the GAO and for the HAM and SM HAM functions determined from different realizations of GLDAS models and from the GRACE RL05 data are shown in Fig. 7a, b. These spectra were computed using Fourier Transform Band Pass Filter [FTBPF; (Kosek 1995)] in the broadband with a 50 and 450 day cutoff. The spectra of all considered functions show annual, semiannual and 120 day period oscillations, both in the prograde and retrograde band.

In the spectra of both global HAM and SM HAM excitation functions, GLDAS NOAH, GLDAS CLM, GLDAS VIC, GLDAS MOSAIC and the GAO, the most energetic oscillation is the annual oscillation signal (Fig. 7a, b). The prograde and retrograde annual oscillation of the GAO have a greater amplitude than that of the considered HAM functions. This suggests that these hydrological excitation functions based on the different GLDAS realizations do not have sufficient energy to excite the observed polar motion. The comparison of the different GLDAS HAM realizations shows that the annual power of GLDAS NOAH is consistent with the GLDAS MOSAIC hydrological function while the power of GLDAS CLM HAM functions has the lowest annual power, especially in retrograde oscillations. The semiannual oscillations both in the prograde and retrograde band are smaller for all considered HAM excitations. The power of the GLDAS NOAH and MOSAIC excitation functions are closest to that of the GAO. The GRACE excitation does not have a comparable amplitude to the other hydrological excitations in the case of prograde oscillations (exception GLDAS CLM); however, in the case of the retrograde oscillations, the power of the GRACE function is almost the same as the power of the GAO.

Similarly, in the spectra of the global SM HAM excitation function computed from the sum of soil moisture in all layers, the most energetic oscillation is also the annual oscillation (Fig. 7b). The power of the SM HAM excitation functions computed from the soil moisture content is less energetic than in the HAM excitations. The SM HAM excitation functions computed from GLDAS MOSAIC, GLDAS VIC and GLDAS NOAH have comparable prograde annual amplitudes, which is similar to the GRACE power.

Fig. 7
figure 7

Fourier Transform Band Pass Filter (FTBPF) amplitude spectra of the different complex HAM (a) and SM HAM (b) excitation functions computed from different GLDAS models: GLDAS NOAH, GLDAS CLM, GLDAS VIC and GLDAS MOSAIC and the GAO and gravimetric excitation functions

Finally, we analyze the prograde and retrograde annual oscillations using phasor diagrams. We consider the GAO and the different HAM and SM HAM excitation functions determined from the full hydrological models as well as from the sum of the soil moisture in all layers (see Fig. 8a, b). We computed these prograde and retrograde components of polar motion excitation functions, \(\chi _1+i\chi _2\) for each of the considered hydrological and gravimetric excitation functions and for the GAO using the least squares method (Brzeziński et al. 2009). The fitted model comprises the second order polynomial and the sum of complex sinusoids with periods of 365.25, 180.0 and 120.0 days. These results suggest that the hydrological GLDAS excitation functions are characterized by similar amplitudes and phases. For the retrograde annual oscillations, a good agreement is present between the GAO and the HAM GLDAS excitation functions. In particular, the best agreement occurs between the GLDAS MOSAIC model and the GLDAS NOAH model. Similarly, the hydrological excitation functions computed for the total sum of soil moisture in all layers agree better with GAO in the retrograde annual oscillations (see Fig. 8b). The excitation functions computed from soil moisture data are in good agreement with the gravimetric excitation function as well as with the GAO and the GLDAS NOAH model. It is also worth noting that annual amplitudes of HAM and SM HAM excitation functions computed from the GLDAS CLM are very small.

Fig. 8
figure 8

Phasor diagrams of the prograde and retrograde annual oscillations, determined using the least squares method, of the GAO, gravimetric excitation functions GRACE CSR and of the a different global HAM excitation functions (GLDAS CLM, GLDAS MOSAIC, GLDAS VIC, GLDAS NOAH) as well as of the b different SM HAM excitation functions determined for different quantity of a soil moisture given in each GLDAS model separately. The data model fitted and removed from the time series comprises the order polynomial and a sum of complex sinusoids with periods 365.25, 180.0 and 120.0 days. Analysis was conducted in the period 2002.3 to 2014.9. Phase \(\phi \) is defined here by the annual term as \(sin(2\pi (t-t_0)+\phi )\), where \(t_0\) is a reference epoch for January 1 2002

5 Discussion

The results of this study indicate that regional HAM and SM HAM and gravimetric excitation functions of polar motion of the complex-valued components \(\chi _1+i\chi _2\) have strong variability in similar land areas, excluding the results from the GLDAS CLM model. There are some exceptions in the Greenland pattern, where hydrological and gravimetric excitation functions of polar motion did not occur in the GLDAS CLM model and in the GRACE data. However, this fact cannot be a determining factor for the proper determination HAM excitation functions. The power of the GLDAS CLM model is lower than that of other GLDAS models (Figs. 1, 2). The time series determined for the four selected GLDAS models are dominated by the annual oscillations (Figs. 3, 4). There are only small differences between the HAM and SM HAM excitation functions determined from the GLDAS NOAH and MOSAIC and VIC models (Figs. 4, 5). In addition, the correlation coefficient determined between these HAM functions is higher than 0.70 (Tables 2, 3).

Comparison between the GAO and global hydrological and gravimetric excitation functions of polar motion indicates that the amplitudes of HAM variations are smaller than the variations of the GAO at seasonal and non-seasonal time scales (Figs. 3, 4, 5, 6), except the HAM GLDAS NOAH and HAM GLDAS MOSAIC excitation functions in the case of the \(\chi _2\) components. Therefore, the HAM hydrological and gravimetric excitation functions do not have enough energy to substantially improve the agreement between geodetic and total geophysical fluid (AAM + OAM + HAM) excitation functions. Although the retrograde annual power of the GLDAS MOSAIC and GLDAS NOAH HAM excitation functions is closest to the retrograde annual retrograde power of the GAO, none of these models considered HAM and SM HAM and gravimetric excitation functions to have enough power to achieve full agreement with the GAO in both prograde and retrograde annual oscillations (Fig. 7a, b).

The phasor plots for annual signal of different global hydrological and gravimetric excitations shown in Fig. 8 indicate that the phases and amplitudes of annual oscillations are compatible with the GAO in retrograde oscillations for HAM and SM HAM excitation functions determined from the GLDAS NOAH and GLDAS MOSAIC models. The largest agreement occurs for amplitudes of the prograde annual oscillation for the GLDAS NOAH and MOSAIC models; however, these functions do not agree well in phase. Improvement of the angular momentum time series is still needed to improve our understanding of the hydrological contributions of the subsystems to angular momentum and hence to the Earth rotation variations.

Further studies of geophysical excitation of the observed polar motion are needed and should focus on determining which of the HAM estimates appears to be most reliable.