Are Ocean Reanalyses Useful for Earth Rotation Research?

Abstract

Oceanic circulation and mass-field variability play important roles in exciting Earth's wobbles and length-of-day changes (ΔΛ), on time scales from days to several years. Modern descriptions of these effects employ oceanic angular momentum (OAM) series from numerical forward models or ocean state estimates, but nothing is known about how ocean reanalyses with sequential data assimilation (DA) would fare in that context. Here, we compute daily OAM series from three 1/4° global ocean reanalyses that are based on the same hydrodynamic core and input data (e.g., altimetry, Argo) but different DA schemes. Comparisons are carried out (a) among the reanalyses, (b) with an established ocean state estimate, and (c) with Earth rotation data, all focusing on the period 2006–2015. The reanalyses generally provide credible OAM estimates across a range of frequencies, although differences in amplitude spectra indicate a sensitivity to the adopted DA scheme. For periods less than 120 days, the reanalysis-based OAM series explain ∼40%–50% and ∼30%–40% of the atmosphere-corrected equatorial and axial geodetic excitation, similar to what is achieved with the state estimate. We find mixed performance of the reanalyses in seasonal excitation budgets, with some questionable mean ocean mass changes affecting the annual cycle in ΔΛ. Modeled excitations at interannual frequencies are more uncertain compared to OAM series from the state estimate and show hints of DA artifacts in one case. If users are to choose any of the tested reanalyses for rotation research, our study points to the Ocean Reanalysis System 5 as the most sensible choice.

Introduction

Mass motion and redistribution in geophysical fluids excite Earth's wobbles and length-of-day changes with varying efficacy on different time scales. While atmospheric contributions to these rotational fluctuations are relatively well determined through reanalysis data sets (Bizouard & Seoane, 2010; Gross et al., 2003, 2004; Neef & Matthes, 2012; Schindelegger et al., 2013), uncertainties increase as one invokes—apart from the somewhat elusive core processes (Kuang et al., 2019; Mound, 2005; Pais & Hulot, 2000)—oceanic (Harker et al., 2021; Marcus et al., 1998; Ponte et al., 1998; Quinn et al., 2019; Zhou et al., 2005), hydrological (Adhikari & Ivins, 2016; Meyrath & van Dam, 2016; Nastula et al., 2019), and cryospheric (J. L. Chen, Wilson, et al., 2013; Göttl et al., 2021) effects in the planet's angular momentum budget. Here we are primarily concerned with modeling the non-tidal oceanic component in polar motion and length-of-day excitations, on time scales from a few days out to several years. Different to all previous studies on the subject (see Harker et al., 2021; Quinn et al., 2019, for recent works), we draw estimates of oceanic angular momentum (OAM) from novel ocean reanalyses, rather than free-running numerical forward models or ocean state estimates. These ocean reanalyses are conceptually equivalent to atmospheric reanalyses, which emerged more than 20 years ago, but they have remained untested for Earth rotation applications.

Ocean reanalyses are based on an ocean general circulation model, which is fitted to in situ and satellite observations by means of data assimilation (DA) (Storto et al., 2019; Wunsch & Heimbach, 2013). During the applied DA, reanalyses often use filter approaches, which vary depending on the exact purposes of the reanalysis. In oceanography, common approaches are for example, some variant of Kalman Filter or three dimensional variational assimilation (3D-Var), where the state of the ocean is estimated sequentially at discrete times. Here, available observations are connected with model states, which can be both forecasts and backgrounds and where, due to several assimilation cycles, information from past observations are included. However, ocean reanalyses can be impacted by unphysical or abrupt changes in state variables, when the forward-propagating model state is allowed to jump toward observations to enforce consistency between model and data within given uncertainties (Pilo et al., 2018; Storto et al., 2019; Wunsch & Heimbach, 2013). This may be potentially critical for Earth rotation research, since violation of conservation laws and kinematic inconsistencies are likely to project on global budgets (Wunsch & Heimbach, 2013) and globally integrated quantities such as OAM. On the other hand, if the physical inconsistencies are small, ocean reanalyses will usefully complement OAM estimates from other sources, which have their own limitations (e.g., high latency or coarse model resolution).

To test ocean reanalyses for their OAM signals and shed light on the impact of sequential DA on these global kinematic quantities, we use members of an ensemble of eddy-permitting ocean reanalyses, each produced by a different weather prediction and ocean monitoring service. The ensemble members are based on the same numerical ocean model and the same volumes of oceanographic data, but they use different DA schemes. This special product was created with the intent to quantify uncertainties of the deduced ocean state, whereas here we are interested in the credibility and uncertainties of the associated OAM changes. A secondary objective in our work is to combine the reanalysis-based OAM values with excitation estimates for other geophysical fluids and assess how well the total modeled excitation agrees with observed rotation fluctuations on sub-seasonal, seasonal, and interannual time scales. Furthermore, we use a statistical combination of the OAM series from single reanalyses to infer a higher-quality excitation series with reduced levels of noise and systematic error. In the following, we introduce the excitation formalism and mathematical description to combine the OAM functions (Section 2), describe the ocean reanalyses and ancillary data sets (Section 3), discuss the results (Section 4), before drawing conclusions and making suggestions for future improvements.

Mathematical Background Excitation Formalism

The equation of motion, representing the rotational behavior of a non-rigid Earth with respect to a body-fixed reference frame, is the Liouville equation (Moritz & Mueller, 1987), a statement of the conservation of the angular momentum within the Earth system. If external torques are neglected, the Liouville equation implies that changes in the angular momentum of a geophysical fluid (e.g., atmosphere or ocean) are compensated by a corresponding perturbation of the angular velocity vector of the solid Earth. The mechanisms involved in these geophysical “excitations” of Earth rotation variations are (a) mass movements affecting the tensor of inertia, and (b) particle motion that generate additional angular momentum relative to the solid-body reference frame. Since the induced deviations from a state of uniform rotation are small, the Liouville equation is readily linearized and the fluid angular momentum changes are cast as dimensionless excitation functions $$\hat{\chi }$$ and χ₃ (e.g., Barnes et al., 1983; Gross, 2007). [Image Omitted. See PDF] [Image Omitted. See PDF]where A′ is the mean equatorial moment of inertia, C_m is the axial principal moment of inertia of the mantle and the crust, Ω denotes Earth's mean angular velocity, $${\Delta }\hat{I}={\Delta }{I}_{13}+\mathrm{i}{\Delta }{I}_{23}$$ $$(\mathrm{i}\equiv \sqrt{-1})$$ and ΔI₃₃ are time variable increments to the third-column elements in the tensor of inertia, and $$\hat{h}={h}_{1}+\mathrm{i}{h}_{2}$$ and h₃ are the equatorial and axial components of relative angular momentum. Scaling factors and moments of inertia in Equations 1 and 2 are consistent with the common assumption of complete core-mantle decoupling from daily to decadal periods (Gross, 2007; W. Chen, Ray, et al., 2013).

To quantify geophysical fluid effects in observed (“geodetic”) excitation of polar motion and changes in length-of-day ΔΛ, we insert the above excitation functions into system functions. In keeping with the deconvolution approach (Chao, 1985), the polar motion system function is [Image Omitted. See PDF]where $$\hat{p}=\hat{p}(t)={p}_{1}-\mathrm{i}{p}_{2}$$ is the reported position of the conventional reference pole in the terrestrial reference frame and $$\dot{\hat{p}}$$ represents the derivative with respect to time t, $${\hat{\sigma }}_{\mathrm{c}}=2\pi (1+\mathrm{i}/2Q)/{T}_{\mathrm{c}}$$ describes the complex-valued Chandler frequency with values of T_c = 433.0 days and Q_c = 179 for period and quality factor (Gross, 1992). Geophysical contributions to changes in length-of-day are given by [Image Omitted. See PDF]where ΔΛ is also called excess length-of-day.

To evaluate excitation functions from ocean reanalysis output, we use volume integrals over density (ρ) and horizontal velocities ($$u,\,v$$) in zonal and meridional directions. [Image Omitted. See PDF] [Image Omitted. See PDF]

where superscripts m and v indicate mass (i.e., inertia) and motion (i.e., relative angular momentum) terms, a is the mean geocentric radius of the Earth, and (ϕ, λ) represent geographical latitudes and longitudes. The radial integral is implemented by assuming a mean Earth radius of 6,370 km in all integration steps, but the layer thickness varies, meaning that we integrate from the ocean bottom all the way to the free surface η. Throughout the paper, we interchangeably use the terms “excitation function” and “angular momentum function” with the prefix indicating the respective subsystem (e.g., OAM function).

Combination of OAM Functions

In this article, we probe excitation functions from three individual ocean reanalyses, along with a combined series aimed at reducing the noise level (i.e., shortcomings in single-model excitation functions). In generic notation, let x(t) be a combination of different, uncorrelated time series xⁱ(t) [Image Omitted. See PDF]

where i = 1, …, N and the quality of each series is represented by weights [Image Omitted. See PDF]

which are normalized such that ∑_iω_i = 1. Var(ϵ) denotes the variance of the noise, which may be estimated using the three-cornered hat method (Koot et al., 2006). The three-cornered hat method assumes a stochastic process $${\mathbf{x}}^{i}$$ with i = 1, …, N and M samples to consist of signal S and noise ϵⁱ components [Image Omitted. See PDF]

The signal is the same for every time series, but the noise differs. One could take differences between time series to approximate ϵⁱ under the assumption that the noises are uncorrelated. In our case, such an approach is problematic, as the reanalyses are identically configured and rely on the same oceanographic observations (see Section 3.1). Tavella and Premoli (1994) proposed a generalization of the three-cornered hat method for correlated noise components. In this method, the time series xⁱ are stored in a matrix X with dimensions M × N, where each column contains one time series. The expected values, which are determined by [Image Omitted. See PDF]are similarly stored in a matrix $$\overline{\mathbf{X}}$$, where the columns again represent time series. In addition, a matrix R is then defined by [Image Omitted. See PDF]representing the covariance of the individual noises, that is, R^ij = Cov(ϵⁱ, ϵ^j). Different to Tavella and Premoli (1994) and Koot et al. (2006), who pursued numerical approaches, we compute the covariance matrix R by applying Equation 11 directly. Working backwards to Equation 7, we deduce combined excitation series for each coordinate direction, and for mass and motion terms separately. The underlying weights, listed in Table 1, are time-invariant and assume different numerical values across the three reanalyses. In particular, the weights of the motion terms $$\left({\chi }_{1,2,3}^{\mathrm{v}}\right)$$ and the axial mass term $${\chi }_{3}^{\mathrm{m}}$$ deviate considerably from 1/3, meaning that the combined series is not a simple average of the three individual series. As an example, the combined equatorial and axial motion terms are comparatively more influenced by one reanalysis (GLORYS, cf. Section 3.1, ω_i ≈ 0.4) than by the other two reanalyses (ω_i ≈ 0.3). Supplemental checks of the behavior of weights in different spectral bands showed that the ω_i cover a wider range of values on interannual time scales (∼0.2–0.5 for both mass and motion terms) than in the sub-seasonal band (∼0.3–0.4). However, for simplicity, we consider all frequency bands at once when combining the reanalyses. The resulting time series are referred to as “Combination” hereinafter.

Table 1 Weights Per OAM Function Component for the Determination of the Combined Series

	GLORYS	ORAS	FOAM
$${\chi }_{1}^{\mathrm{m}}$$	0.340	0.335	0.325
$${\chi }_{2}^{\mathrm{m}}$$	0.334	0.320	0.345
$${\chi }_{3}^{\mathrm{m}}$$	0.294	0.519	0.187
$${\chi }_{1}^{\mathrm{v}}$$	0.406	0.360	0.234
$${\chi }_{2}^{\mathrm{v}}$$	0.415	0.310	0.275
$${\chi }_{3}^{\mathrm{v}}$$	0.447	0.270	0.283

Data Sets Ocean Reanalyses

We derive OAM functions from three out of four members of the eddy-permitting ocean reanalysis ensemble provided by CMEMS (Copernicus Marine Environment Monitoring Service, Desportes et al., 2019). All reanalyses cover the period 1993–2019 and are based on the ocean model NEMO3 (Nucleus for European Models of the Ocean version 3). Sea surface temperature observations, daily sea level anomalies, sea ice concentration as well as temperature and salinity profiles are used to constrain the three-dimensional ocean state as NEMO3 is integrated forward in time. The model itself is configured on a 1/4° horizontal tri-polar grid with each 75 vertical layers. The thickness of these layers increases from 1 m at the surface, to 10 m at 100 m depth and to 200 m at the bottom. Six-hourly buoyancy and momentum fluxes from ERA-Interim are used as common forcing data (Dee et al., 2011).

Our analysis period is from 2006 to 2015, which allows us to examine oceanic excitations and Earth rotation fluctuations with sub-seasonal, seasonal and (to some extent) interannual frequencies. The choice of this 10-year period is dictated by computational constraints, but also by the availability of complementary AAM (atmospheric angular momentum) series and other auxiliary data sets (e.g., satellite gravimetry). From the CMEMS product, we use the Global Ocean Reanalysis and Simulation 2 version 4 (GLORYS2v4, for short: GLORYS), the Ocean Reanalysis System 5 (ORAS5, for short: ORAS) and the Forecast Ocean Assimilation Model - Global Seasonal forecast system version 5 (FOAM-GloSea5, for short: FOAM). These ensemble members were selected mainly because GLORYS and ORAS use different DA schemes, whereas ORAS and FOAM use a similar DA scheme but different assimilation windows, which may or may not impact the reconstructed state and thus OAM quantities. GLORYS employs the SAM2 (Système d’Assimilation Mercator version 2) method based on a singular evolutive extended Kalman filter (SEEK) formulation and a 7-day assimilation window. Future and past observations relative to the window mid-point are used to perform the analysis in 7-day intervals (Garric & Parent, 2017; Lellouche et al., 2013). The DA software for ORAS and FOAM is NEMOVAR, an incremental three-dimensional variational assimilation approach (Blockley et al., 2014; Zuo et al., 2017). One difference between these two reanalyses consists in the assimilation window, which is 5 days for ORAS and 1 day for FOAM (cf. Table 2). In addition, ORAS accounts for representation errors in observation and structure, as well as analysis errors in surface forcing. These uncertainty estimates were derived by perturbing initial conditions, observations, and forcing, and performing the ocean state reconstruction for a total of five times (Zuo et al., 2019).

Table 2 Selected Components of the Ocean Reanalyses Used in This Study^a

Reanalysis	GLORYS2V4	ORAS5	FOAM-GloSea5v13
Processing center	Mercator Ocean	ECMWF	UK Met Office
Ocean model	NEMO3.1	NEMO3.4.1	NEMO3.4
Surface nudging	No, but flux correction on precipitation	Yes	Yes
DA scheme	SAM2 (SEEK)	NEMOVAR (3D-Var)	NEMOVAR (3D-Var)
	7-day assimilation window	5-day assimilation window	1-day assimilation window
Uncertainties	Observation and background error variances from statistical methods (Lellouche et al., 2013)	Representation errors in observation and structure and analysis errors (Zuo et al., 2017)	Observation and background error variances from statistical methods (Blockley et al., 2014)

^aTable adapted from Desportes et al. (2019).

All three reanalyses use climatological, seasonally varying river discharge, and GLORYS additionally considers seasonal ice shelf discharge. Combined with evaporation minus precipitation over the ocean, the continental freshwater input leads to global ocean mass fluctuations, which are relevant for excitations of ΔΛ. In this context, the surface nudging scheme of the reanalyses becomes interesting (Table 2). Surface nudging is not applied in GLORYS in favor of a flux correction on precipitation. The correction corresponds to an addition or removal of a thin, spatially uniform layer of mass at each analysis step—a reasonably accurate approach given the strong tendency toward an equilibrium response to loading by variable freshwater fluxes (Ponte, 2006).

From the selected reanalyses, we use daily fields of potential temperature and salinity to compute density ρ, eastward and northward velocities and the sea surface height η for calculating the angular momentum functions (Equations 5 and 6). The data are available at CMEMS (Copernicus Marine Service, 2019). Consistent with the model bathymetry adopted in the three reanalyses, we use a 1/4°-averaged version of the 60-arcminute ETOPO1 data set (Amante & Eakins, 2009; NOAA National Geophysical Data Center, 2009) as lower bound in the vertical integration of dynamical fields (Equations 5 and 6). Because the reanalyses treat ice-shelf cavities as land, the water bodies underneath ice shelves do not contribute to the global OAM integrals.

ECCOv4

OAM mass and motion terms from the ECCOv4 (Estimating the Circulation and the Climate of the Ocean) Version 4 Release 3 state estimate (ECCOv4 for short) are used as a point of comparison in this study. The ECCOv4 state estimates are iterative fits of the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al., 1997) to most oceanic in-situ and satellite data, including in situ hydrographic profiles, Argo float observations, sea surface height estimates from satellite altimetry, and bottom pressure anomalies from the Gravity Recovery and Climate Experiment (GRACE) mission. The adjustment is accomplished via the adjoint method, which propagates model-data misfits—distributed in space and time—to variations of uncertain model inputs, such as initial conditions, forcing fields, and coefficients of subgrid scale parameterisations. A forward integration of the MITgcm under the adjusted inputs yields new model-data misfits, and the procedure is repeated until an acceptable fit is found. Being an exact solution to a general circulation model, the ECCOv4 state estimates preserve dynamic and kinematic consistency, such that there are no spurious sources or sinks of tracers, volume, momentum, and derived global quantities (e.g., OAM). The MITgcm setup underlying ECCOv4 is Boussinesq volume-conserving, has a nominal horizontal resolution of 1°, and contains 50 layers in the vertical, ranging from 10 m thickness at the surface to 460 m at the bottom (ECCO Consortium et al., 2017; Forget et al., 2015). Newer releases, which contain forcing by barometric pressure, are not considered in this study, given that pressure loading is also absent in the ocean reanalyses. From the ECCOv4 OAM series made available by the IERS SBO (International Earth Rotation and Reference Systems Service, Special Bureau for the Oceans, https://isdc.gfz-potsdam.de/ggfc-oceans/, Last access: 14.04.2022), we choose the version that includes net effects of freshwater flux from continental and atmospheric reservoirs into the ocean.

Ancillary Data Sets

A possible means of inferring fluctuations in the mass term associated with terrestrial hydrology and ice bodies is to use satellite gravimetry data of the GRACE mission. From 2002 to 2017, GRACE monitored surface mass changes in the Earth system, available as unconstrained global gravity field solutions (or derived quantities) with a nominal sampling period of 1 month (Tapley et al., 2019). Here, we use Release-06 GRACE mascon gravity fields provided by the Jet Propulsion Laboratory (JPL), which are solutions in terms of mass concentration blocks instead of spherical harmonics (Watkins et al., 2015). The surface mass changes are computed for each equal-area 3° × 3° spherical cap mascon and are then eventually sampled to a 0.5° latitude-longitude grid. A process-based Coastline Resolution Improvement filter is employed to separate mass signals near the land-ocean boundary by source region. As is standard, the JPL processing includes a replacement of the zonal degree-2 coefficient with a more accurate estimate from Satellite Laser Ranging (SLR). We subset the monthly mascon solutions to total water storage anomalies (TWSA) over land and ice sheets by clipping oceanic areas, using a land-ocean mask provided with the data set. Over the 10 years considered, the GRACE data contains 13 gaps of 30 or 60 days duration, but none between 2006 and 2010. We fill these gaps with seasonal sinusoids, estimated from available epochs at each location, plus a residual derived from an principal component analysis of globally gridded non-seasonal fluctuations. The resulting gridded TWSA are readily converted into angular momentum functions, which we abbreviate as “hydrology/ice” in figures, tables, and some of the text below.

For consistency with the atmospheric forcing in the ocean reanalyses, AAM series are based on ERA-Interim data (see Schindelegger et al., 2013, for details). We use mass term estimates corrected for the inverted barometer (IB) effect and time-average 6-hourly to daily values centered at midnight.

The Earth rotation data in our excitation budget considerations are the SPACE2018 series by Ratcliff and Gross (2019), deduced from a Kalman filter-based combination of various space geodetic measurements (Very Long Baseline Interferometry, Global Positioning System, SLR). Here, we take daily polar motion and ΔΛ estimates sampled at midnight. From the ΔΛ values, we remove tidal contributions at 80 spectral lines using the model by Ray and Erofeeva (2014). Similarly, long-period tidal effects in polar motion are accounted through the conventional model (Petit & Luzum, 2010), with the fortnightly component replaced by the Mf solution of Ray and Egbert (2012).

Results Signal Content

To illustrate characteristics of the OAM functions across a range of frequencies, we show the amplitude spectra for equatorial mass and motion terms in Figure 1. Since the axial terms are small in magnitude and discussed as time series to some extent below, they are omitted. Estimates of power P were first computed with a 512-point Fast Fourier Transform, using the method of Welch (1967), and then converted into amplitudes $$A=\sqrt{2P}$$.

Figure 1. Amplitude spectra (in mas) for the equatorial mass terms χ1,2m ${\chi }_{1,2}^{\mathrm{m}}$ (a, b) and motion terms χ1,2v ${\chi }_{1,2}^{\mathrm{v}}$ (c, d) for ECCOv4, the three reanalyses, and the combined series, which is dubbed “Combination.” The RMS of each time series (in mas) is specified on the right side of the panels.

Both mass terms (Figure 1, panels a and b) follow a red spectrum, but the $${\chi }_{2}^{\mathrm{m}}$$ component is more energetic than $${\chi }_{1}^{\mathrm{m}}$$ for periods longer than 20 days; cf. Gross et al. (2003) and Harker et al. (2021). The ECCOv4 and ORAS mass terms bear a close resemblance to each other in all bands, except near the annual frequency in $${\chi }_{1}^{\mathrm{m}}$$, where ORAS amplitudes exceed those of ECCOv4 by ∼2 mas (milliarcseconds). GLORYS exhibits 2-day to 5-day oscillations not seen in other reanalyses, while FOAM has enhanced variability on interannual time scales. Neither feature is carried through to the Combination, as it suppresses fluctuations in OAM when covariance among the reanalyses is lacking. In $${\chi }_{1}^{\mathrm{m}}$$ (Figure 1a), the spectra of all reanalyses are similar, except for the mentioned variability of GLORYS. Evident in all estimates of $${\chi }_{2}^{\mathrm{m}}$$ (Figure 1b) is a suite of oscillations with periods between 20 and 50 days, interrupted by a pronounced trough of energy around 40 days. These spectral patterns are likely related to strong barotropic variability in the Bellingshausen Basin (see Fu, 2002; Fukumori et al., 1998; Weijer, 2015) and Australian-Antarctic Basin (Harker et al., 2021; Weijer, 2010), both leading to a clear signal in the χ₂ component.

Turning to the motion term $${\chi }_{1}^{\mathrm{v}}$$ (Figure 1c), FOAM stands out with enhanced power across all frequencies but particularly at periods between 20 and 100 days. By contrast, the spectra of GLORYS and ORAS are less steep at intraseasonal periods and only moderately pick up power beyond the annual frequency. The Combination contains weighted information of all reanalyses with FOAM being less weighted than the other two reanalyses, resulting in spectral characteristics comparable to those of GLORYS and ORAS. Variability in the $${\chi }_{2}^{\mathrm{v}}$$ component (Figure 1d) is distributed over frequency in a similar fashion across all models (including ECCOv4), with noticeable cusps of energy near 50-day and 90-day periods. Other features, including a broad and comparatively energetic peak in ORAS around the annual frequency, remain specific to only one reanalysis.

In summary, the OAM series from the three reanalyses are far from identical, despite their origin in a common numerical model, constrained by the same oceanographic data. The differences seen in Figure 1 thus reflect the impact of the DA scheme (Table 2) or global parameter choices (e.g., drag coefficients, to which we have no access to). It is very unlikely that the freshwater flux schemes (Table 2) cause some of the $$\hat{\chi }$$ differences, since the contribution of these mass fluxes to the OAM variance in the equatorial direction does not exceed 1.5% (value estimated from comparing two ECCOv4 OAM solutions with and without freshwater loads, Quinn et al., 2019). Moreover, currents involved in the response to freshwater fluxes are restricted to time scales of a few days (Durand et al., 2019), implying that differences among reanalyses in the motion terms must be explained by other processes. In fact, the $${\chi }_{1}^{\mathrm{v}}$$ and $${\chi }_{2}^{\mathrm{v}}$$ spectra provide strong indications that each DA method introduces its own, somewhat arbitrary perturbations to the dynamical model state, as speculated in Section 1. Below we assess whether these perturbations to OAM quantities are relevant in comparisons to Earth rotation data.

Sub-Seasonal Band

For the sub-seasonal band, we consider periods below the small ter-annual wobble and thus high-pass filter time series at a cutoff frequency of 1/120 days⁻¹, as in Harker et al. (2021). The agreement between geodetic and geophysical excitation is quantified in terms of the root-mean-square (RMS) of their difference and the percentage of variance explained (PVE), as presented in Table 3. To compute the PVE of the oceanic excitation, IB-corrected atmospheric effects are removed from the geodetic excitation and the residual is then compared with the OAM functions. During 2006–2015, the atmosphere accounts for 56.2% of the deconvolved polar motion variance and 93.7% of sub-seasonal ΔΛ signals (cf. Gross et al., 2003; Gross et al., 2004; Harker et al., 2021), such that the RMS of the residuals is 25.7 mas for $$\hat{\chi }$$ and 34.7 μs in the χ₃ component.

Table 3 Excitation Budget for Sub-Seasonal Oscillations 2006–2015^a

	χ1	χ2	$$\hat{\chi }$$	χ3
RMS Observation	21.1	32.9	39.1	140.1
PVE by atmosphere
ERA-Interim (IB)	53.7 (14.3)	57.3 (21.3)	56.2 (25.7)	93.7 (34.7)
PVE by ocean in atmosphere-corrected geodetic excitation
ECCOv4	35.6 (11.5)	45.0 (15.9)	42.1 (19.6)	30.9 (29.0)
GLORYS	27.5 (12.2)	50.3 (15.2)	43.3 (19.5)	28.8 (32.7)
ORAS	40.4 (11.1)	51.9 (14.9)	48.4 (18.6)	29.9 (27.8)
FOAM	−1.1 (14.5)	33.6 (17.6)	22.9 (22.7)	−15.8 (37.8)
Combination	42.0 (11.0)	56.2 (14.3)	51.8 (18.0)	40.1 (27.2)
PVE by sum of atmosphere and ocean
ERA-Interim + ECCOv4	70.2	76.5	74.7	95.7
ERA-Interim + GLORYS	66.4	78.8	75.2	95.8
ERA-Interim + ORAS	72.4	79.5	77.4	95.5
ERA-Interim + FOAM	53.2	71.6	66.2	92.7
ERA-Interim + Combination	73.1	81.3	78.9	96.2
PVE by secondary terms in residual series
Hydrology/Ice	−0.5	−1.7	−1.2	0.4

^aValues are PVE, except for the first line, and the corresponding RMS of residuals is in parentheses. Units are [mas] for equatorial terms and [μs] for the axial term.

Looking in detail at the statistics of the OAM solutions in Table 3, we see that the variance reduction in the equatorial components is greater than in the axial component, as is well known (Gross et al., 2003, 2004). For most models, the RMS of residuals drops to ∼19 mas in $$\hat{\chi }$$ and ∼29 μs in χ₃. FOAM forms an exception and produces negative PVE values in the χ₁ and χ₃ components, meaning that the variance increases when removing FOAM from the atmosphere-corrected geodetic excitation. One interpretation of this result could be that the frequent (daily) adjustments of the FOAM forward model toward the data increase the amount of unphysical signals in the reanalysed state and thus the noise in the OAM series. As for the other OAM data sets, ECCOv4 has a slightly lower PVE in $$\hat{\chi }$$ compared to GLORYS and ORAS, while in the axial component the PVE for all three data sets hovers around 30%. Hence, some of the reanalyses, despite their possible shortcomings, match or improve upon a dynamically consistent state estimate in excitation budgets. This result is an indication for the benefits of high horizontal model resolution (1/4° in the reanalyses vs. 1° in ECCOv4) consistent with findings in Harker et al. (2021). More to the point, by conflating the three reanalysis series in a statistical meaningful way (Section 2.2), we achieve highest variance reduction with the Combination, amounting to 51.8% $$\left(\hat{\chi }\right)$$ and 40.1% (χ₃). Accordingly, when comparing the sum of atmosphere and ocean with the geodetic excitation, the Combination also has the highest PVE, cf. Table 3. In detail, we find variance reductions of up to 78.9% and 96.2% for $$\hat{\chi }$$ and χ₃, suggesting that atmospheric and oceanic excitations explain a large fraction, but not all of the observed variability on sub-seasonal time scales.

When adding hydrological and cryospheric contributions to these considerations, the PVE in residual series of geodetic, atmospheric and oceanic excitation are very small. Values are negative in the equatorial component and they do not exceed 0.4% in the axial direction. Thus, secondary excitation processes in other Earth system components—as represented by the coarse-resolution GRACE fields—fail to account for the residual geodetic excitation. The remaining gaps in Earth's sub-seasonal rotation budget may be rather due to missing effects of barometric pressure loading in the ocean models, over- or underestimation of tropospheric winds (Schindelegger et al., 2013), imperfect representation of (a) atmospheric planetary waves, (b) mass exchange across ocean basins (Afroosa et al., 2021) or (c) topographically constrained oceanic excitations (Harker et al., 2021), or simply high-frequency noise in the rotation data (Dill et al., 2020).

Seasonal Oscillations

Figure 2 shows phasor plots for annual and semi-annual oscillations in polar motion excitation, as deduced from a least squares harmonic analysis. Here, we compare the geodetic excitation, corrected for hydrology/ice effects, with the combined atmosphere-ocean excitation signals for the different OAM solutions. In the prograde annual component (Figure 2a), imposing the Combination or any of the reanalyses onto the atmospheric excitation gives a phasor sum that is very close to the observation, to within 0.3–1.5 mas. In comparison, the ECCOv4 estimate has the wrong phase and is too large in magnitude (∼4 mas), similar to previous ECCO state estimates (Gross et al., 2003). In the retrograde annual wobble excitation (Figure 2b), all OAM solutions, including ECCOv4, form a cluster with phases between −2° to 52°, meaning that they point away from the corrected geodetic excitation when added to the ERA-Interim phasor. Owing to enhanced power in the mass term, GLORYS and ORAS have the largest amplitudes in the retrograde annual band, more than twice as large as the amplitude of FOAM. Complementary to these illustrations, we list amplitude and phase estimates for annual harmonics of polar motion excitation from individual fluids, various sums of them, and rotation data in Table A1. Compared with earlier such decompositions (e.g., Gross et al. (2003) for the period 1980–2000, see also Dobslaw et al. (2010)), the annual wobble is more energetic over the time span analyzed here. Our results indicate that this is due to an enhanced atmospheric contribution in the prograde part. The ∼12 mas amplitude in the retrograde part probably arises from several excitation processes or simply a favorable phase constellation, but the imperfect budget closure in that component prohibits more solid conclusions.

Figure 2. Phasor plots for the annual (a, b) and semi-annual (c, d) polar motion excitation during 2006–2015, split up in prograde (left column) and retrograde oscillations (right column). Geodetic excitation minus hydrology/ice (black phasors) should be compared with the sum of atmospheric excitation (gray) and oceanic excitation from any of five sources (ECCOv4, GLORYS, ORAS, FOAM, Combination). Reference date for phase is 1 January 2006, 12 UTC.

Looking at the bottom panels in Figure 2, polar motion excitation signals at the semi-annual frequency are generally a factor of 3–4 smaller than at the annual frequency. However, oceanic effects remain at 2–3 mas in magnitude (see Table A1), turning them into the single most effective excitation process for the prograde semi-annual wobble over the 10-year period considered. All reanalyses and ECCOv4 yield very similar phasors in the prograde band, which—when added to the atmospheric contribution—agree with the corrected geodetic excitation to within 1.5 mas. By contrast, adding the oceanic estimates for retrograde semi-annual oscillations to that of the atmosphere (Figure 2d) draws the modeled excitation farther away from the observation. The magnitude of the disparity (2.1 mas) is about 40% of the observed excitation amplitude, suggesting that there are still considerable errors in presently available angular momentum data sets. Nonetheless, and of interest in the context of this study, the seasonal polar motion excitation signals deduced from ocean reanalyses are broadly consistent with the adopted reference solution (ECCOv4) except for the prograde annual term.

Results from least squares adjustment of seasonal sinusoids in the χ₃ component are depicted in Figure 3. In the annual component (Figure 3a), the length-of-day observation ΔΛ, corrected for hydrology/ice effects, is largely consistent with the sum of ERA-Interim and all OAM solutions, although discrepancies remain. The χ₃ signal is dominated by the mass term $$\left({\chi }_{3}^{\mathrm{m}}\right)$$, which we list in Table 4, along with the counteracting contribution from hydrological and (small) cryospheric effects (Dill & Dobslaw, 2019). All OAM data sets suggest a $${\chi }_{3}^{\mathrm{m}}$$ peak in October (phase ∼280°), but amplitudes vary from 69 μs in ECCOv4 to about 105 μs in GLORYS and FOAM. This difference is potentially revealing shortcomings in the two reanalyses, as ECCOv4 incorporates monthly GRACE solutions of ocean bottom pressure (OBP) and should therefore provide the most credible estimate of $${\chi }_{3}^{\mathrm{m}}$$. Before exploring the matter somewhat further, we take a brief look at the semi-annual χ₃ oscillation in Figure 3b. Atmospheric effects account for nearly all of the observed length-of-day change, leaving a residual of merely ∼15 μs, similar to what we see at the annual frequency. However, the hydrological excitation signal at the semi-annual period is weak (and also uncertain, 8.0 ± 2.7 μs), implying that only a very small oceanic contribution is required to close the excitation budget. In this light, the spread of oceanic χ₃ estimates in Figure 3b is not too surprising. ECCOv4, ORAS, and FOAM phasors generally point toward the hydrology-corrected ΔΛ estimate, although the magnitude appears to be too small and the separation into mass and motion term contributions remains unclear (cf. the diverging $${\chi }_{3}^{\mathrm{m}}$$ estimates in Table 4). Most anomalous in the semi-annual band is the GLORYS phasor, which is governed by a very large signal in the mass term (20 μs) that inevitably carries through to the Combination.

Figure 3. Phasor plots for the annual (a) and semi-annual (b) ΔΛ excitation during 2006–2015. Geodetic excitation minus hydrology/ice (black phasors) should be compared with the sum of atmospheric excitation (gray) and oceanic excitation from any of five sources (ECCOv4, GLORYS, ORAS, FOAM, Combination). Reference date for phase is 1 January 2006, 12 UTC.

Table 4 Phase and Amplitude for the Annual and Semi-Annual $${\chi }_{3}^{\mathrm{m}}$$ Component

	Annual		Semi-annual
	Amplitude [μs]	Phase [°]	Amplitude [μs]	Phase [°]
ECCOv4	68.88 ± 0.68	277.23 ± 0.56	5.42 ± 0.68	188.15 ± 7.14
GLORYS	102.49 ± 0.95	284.65 ± 0.53	21.98 ± 0.94	62.53 ± 2.46
ORAS	75.02 ± 0.82	280.53 ± 0.62	2.90 ± 0.82	330.16 ± 16.09
FOAM	105.42 ± 1.77	278.37 ± 0.96	1.10 ± 1.76	348.16 ± 92.09
Combination	88.69 ± 0.80	281.45 ± 0.52	6.68 ± 0.80	47.75 ± 6.89
Hydrology/Ice	60.40 ± 2.65	108.51 ± 2.53	7.96 ± 2.66	136.74 ± 19.14

Shifting the focus again to the annual oscillation, we illustrate time series of $${\chi }_{3}^{\mathrm{m}}$$ from the three reanalyses and ECCOv4 in Figure 4a, where a second-order peaking filter (Orfanidis, 1996) was applied to all excitation functions to extract the time-variable annual signal. In keeping with Table 4, the $${\chi }_{3}^{\mathrm{m}}$$ series from ECCOv4 and ORAS agree very well with each other (save a 1-year period in 2011–2012), while FOAM and GLORYS have excess amplitudes of a few tens of μs and more variable phases throughout. To examine the anomalous behavior in more detail, we include in Figure 4a a plot of the daily sampled FOAM excitation function before applying the peaking filter. The series has several spikes (both positive and negative) and abrupt transitions, for example, in early 2007, late 2009 and 2011, or throughout the year 2013, which clearly contribute to the large annual amplitude of $${\chi }_{3}^{\mathrm{m}}$$ from FOAM in Table 4. Erratic features of this kind also occur in the GLORYS OAM function, although the intraseasonal noise is smaller than in FOAM (not shown).

Figure 4. (a) Annual excitation signals in χ3m ${\chi }_{3}^{\mathrm{m}}$ for 2006–2015 derived from ECCOv4, GLORYS, ORAS and FOAM, and the daily OAM series of FOAM (arbitrary offset applied). In panel (b), we show mass term excitation functions computed from global ocean mass fluctuations in FOAM (black curve) and ECCOv4 (olive curve) with daily sampling. The underlying FOAM ocean mass series is also plotted, with scale and units (cm of equivalent water height, EWH) indicated by the y-axis on the right. Signal content with periods longer than 12 months has been removed from all time series.

Given that $${\chi }_{3}^{\mathrm{m}}$$ is sensitive to changes of the global ocean mass (see, e.g., Yan & Chao, 2012), one hypothesis to test is whether erroneous freshwater fluxes received from the atmosphere account for the large annual oscillation in FOAM and GLORYS. This is a physically plausible thread to follow, because the local mass gained or lost via atmospheric freshwater fluxes is spread evenly over the global ocean surface in a matter of days (Durand et al., 2019). Hence, we computed changes in the global ocean mass in FOAM and quantified—using a vertically integrated form of Equation 6—how much it contributes to $${\chi }_{3}^{\mathrm{m}}$$ in Figure 4a. The resulting time series, illustrated in Figure 4b (blue curve), is of appreciable variance and shares many features with the actual $${\chi }_{3}^{\mathrm{m}}$$ function from FOAM, but fails to echo a number of spikes (e.g., in spring 2007, 2011, and 2013). Thus, deficiencies in the axial OAM mass term in some of the reanalyses can only be partially associated with erroneous freshwater fluxes. Further suspects are a too energetic response to atmospheric wind stress torques (cf. Ponte et al., 2001) or spurious dynamics incurred during DA, which is not an unreasonable assumption given the results with FOAM in other parts of our analysis.

Interannual Variability

Figure 5 shows the interannual signals of oceanic excitation computed from ECCOv4, the three reanalyses and the Combination in comparison to geodetic observations corrected for atmospheric, hydrological and cryospheric effects. Both AAM and terrestrial water storage changes are important excitation processes on interannual time scales, which is why we use them for the computation of the residuals. To isolate the interannual signal, we have applied a high-pass filter with a cutoff frequency of 1/365 days⁻¹ after removing the mean, trend and the seasonal oscillations estimated in a least squares adjustment. The interannual signal was then obtained by subtracting the high-pass-filtered series from the seasonally corrected excitation function. For all solutions and all components, the mass term dominates over the motion term, typically by factors of 1.3 in χ₁, 4 in χ₂, and 10 in χ₃. With regard to the GRACE-based mass terms used to correct the geodetic excitation, it should be kept in mind that these contain both hydrological and cryospheric effects, but additional tests by us have shown that the contribution from ice has a much smaller impact on interannual time scales than land hydrology.

Figure 5. Interannual excitation signals χ1,2,3 as deduced from ECCOv4, GLORYS, ORAS, FOAM, the Combination and GRACE-based hydrology/ice in comparison with the geodetic excitation–AAM–HAM from 2006 through 2015. Four months were clipped both at the beginning and at the end of the time series to avoid filter artifacts.

From Figure 5 we see that both in χ₁ and χ₂, ECCOv4 agrees reasonably well with the reduced observations. Likewise, the Combination—providing a middle ground to the three reanalyses—captures most of the peaks and troughs in the reduced observations, especially before 2013. FOAM has the largest fluctuations in time, which is consistent with the amplitude spectra shown in Figure 1. However, the pronounced positive excitation signals in FOAM around 2010 in χ₁ (∼10 mas) and around 2012 in both χ₂ (∼15 mas) and χ₃ (∼100 μs) have no correspondence in the observations, suggesting that they are spurious and possibly DA artifacts.

From early to late 2009, a sharp decline of χ₂ is evident in all solutions. We interpret this feature as the manifestation of a positive OBP anomaly in the ocean, projecting onto a negative OAM anomaly due to the longitude of its likely source region. In fact, Boening et al. (2011) found a record increase in GRACE-based Southern Ocean OBP in the Bellingshausen Basin, caused by an anomalous anti-cyclonic flow potentially related to El Ni$$\tilde{\mathrm{n}}$$o in 2009/2010. Another clear and physically motivated excitation signal is a large (±200 μs) oscillation with a period of 5–6 years in the reduced length-of-day observations (χ₃ component). The causative mechanism is the exchange of angular momentum between the mantle and the core, mostly governed by the gravitational coupling between the two interior components (J. Chen et al., 2019).

Continuing the budget considerations of the preceding Sections, Table 5 shows the PVE and the correlation coefficients for excitations on interannual time scales (2006–2015). As before, the PVE is computed by comparing residual series with the oceanic excitation, while for the computation of the correlation coefficients, we use Pearson's R. The PVE results reflect what can be seen in Figure 5. ECCOv4 explains the atmosphere-corrected polar motion excitation by more than 72%, whereas the three reanalyses have smaller explained variances. Yet again, our combined OAM function mitigates noise and systematics of individual reanalysis series, such that in the χ₁ component, the PVE of the Combination (77%) is higher than that of the reanalyses (17–55%) and also ECCOv4. In χ₂, the skill of ORAS is only slightly worse than that of the Combination, and similar margins are seen in the complex-valued component $$\hat{\chi }$$. Residual ΔΛ signals are explained to some extent by ECCOv4, GLORYS and ORAS (10–23%), while the PVE of FOAM is negative. In fact, FOAM performs worst in all three components, with the discrepancies being caused by the large, and arguably unphysical long-period fluctuation apparent in Figure 1 (panels a, b) and in Figure 5.

Table 5 Modeled Versus Observed Excitations on Interannual Time Scales, 2006–2015^a

	χ1	χ2	$$\hat{\chi }$$	χ3
PVE by atmosphere
ERA-Interim (IB)	28.4 (0.53)	−43.9 (−0.19)	−19.5 (0.06)	21.1 (0.50)
PVE by ocean in atmosphere-corrected geodetic excitation
ECCOv4	72.7 (0.86)	72.4 (0.85)	72.5 (0.85)	10.0 (0.38)
GLORYS	54.3 (0.63)	57.1 (0.38)	56.2 (0.44)	22.5 (0.51)
ORAS	54.8 (0.57)	71.8 (0.72)	66.1 (0.66)	17.7 (0.51)
FOAM	16.6 (0.66)	0.5 (0.83)	5.9 (0.76)	−15.6 (0.09)
Combination	76.8 (0.70)	74.0 (0.81)	74.9 (0.75)	17.2 (0.45)
PVE by sum of atmosphere and ocean
ERA-Int. + ECCOv4	78.8 (0.89)	70.1 (0.85)	73.6 (0.86)	29.9 (0.57)
ERA-Int. + Combination	81.9 (0.79)	71.8 (0.59)	75.9 (0.66)	35.4 (0.59)
PVE by sum of atmosphere, ocean and hydrology/ice
ERA-Int. + ECCOv4 + Hyd./Ice	85.2 (0.92)	84.7 (0.92)	84.9 (0.92)	19.1 (0.49)
ERA-Int. + Combination + Hyd./Ice	87.4 (0.93)	85.6 (0.93)	86.2 (0.93)	25.5 (0.53)

^aValues are PVE and the corresponding correlation coefficient is in parentheses.

For completeness, Table 5, also includes statistics for the direct comparison of observed excitation with the sum of modeled excitations from different geophysical fluids. Joint consideration of atmospheric and oceanic effects yields PVE values of ∼74% and ∼76% in $$\hat{\chi }$$, and ∼30% and ∼35% in χ₃, dependent on whether ECCOv4 or the Combination is used as OAM source. Addition of hydrology and ice contributions reduces the total ΔΛ PVE to 19% and 26%, respectively, but these number should be interpreted with caution, given that we have not removed (or modeled) the ∼5-year oscillation due to core processes. By contrast, incorporating the GRACE-based mass term in the polar motion excitation budget improves the PVE to ∼85% (ECCOv4) and ∼86% (Combination), with the correlation coefficient being as high as 0.93 in the latter case. We conjecture that the residual in the comparison of geodetic and geophysical excitation of Earth's wobbles in Table 5 is due to errors in 3D wind fields and atmospheric forcing data over the ocean, representation errors in ocean models (e.g., bathymetry and the omission of circulation underneath ice shelves) or uncertainties of long wavelength features in the GRACE-based mass terms.

Summary and Outlook

We have evaluated whether or not angular momentum estimates from ocean reanalyses are useful for Earth rotation research. By and large, one can answer in the affirmative, although the quality of the OAM series varies with time scale and the specifics of sequential DA applied. At sub-seasonal periods from 2 to 120 days, the reanalyses offer similar skill in explaining atmosphere-corrected geodetic excitation as an established ocean state estimate (∼43%–52% in $$\hat{\chi }$$ and ∼29%–40% in χ₃), a result partially attributed to the benefits of high horizontal and vertical model resolution. As evident from the (low) PVE values for FOAM in Table 3, reanalyses are not free from error at short periods, but a statistical combination can successfully suppress noise and DA artifacts inherent to such single-reanalysis OAM series. While for annual frequencies, all tested OAM estimates blend in well with atmospheric and hydrological excitations to produce reasonably well-closed rotation budgets, larger discrepancies occur for semi-annual frequencies. Analysis of interannual variability is somewhat limited by the relatively short (10-year) time window, but from results in Section 4.4 it is clear that here, the reanalyses may not always compete with ECCOv4. The spread among reanalyses is appreciable and the large anomaly in FOAM-based OAM series after 2012 (particularly in χ₃) is a dubious feature as highlighted by our excitation budget considerations. If users are to choose any of the tested reanalyses, we recommend to use ORAS by ECMWF. Interestingly, DA products by ECMWF are also a well established source for computing AAM series, and they have been shown to fulfill global kinematic constraints across a range of time scales (Schindelegger et al., 2013).

The encouraging results notwithstanding, there is obvious room for improvement of ocean reanalyses in the context of Earth rotation research. We especially propose to consider the dynamic ocean response to atmospheric pressure loading, which plays an important role in forcing rapid (sub-weekly) rotation signals (Ponte & Ali, 2002). To some extent, it also acts on monthly to interannual time scales, at odds with a perfect IB behavior (Piecuch et al., 2022). Following the example of ECCOv4, it should be relatively straightforward for ocean reanalyses to assimilate monthly GRACE gravity field solutions as bottom pressure observations. Such development could help better constrain the equatorial and axial OAM mass terms, as well as global ocean mass (and thus ΔΛ) fluctuations caused by freshwater fluxes. In this context, a mass balance constraint could be incorporated in each assimilation step to ensure that there are no spurious fluctuations in the total ocean mass due to sequential DA. Only as much mass should be drawn from, or added to the ocean as the net effect of freshwater fluxes from the atmosphere, continental hydrology, and cryospheric components implies. Given that evaporation and precipitation fields from atmospheric models are still afflicted with errors, it would also be desirable to find ways for correcting these fluxes during the DA and bring them into consistency with the reconstructed ocean state (Quinn et al., 2019).

A AppendixSeasonal Prograde and Retrograde Oscillations in $$\hat{\chi }$$

Table A1 shows the annual and semi-annual amplitudes and phases for prograde and retrograde oscillations in $$\hat{\chi }$$ estimated in a least squares adjustment.

Table A1 Annual and Semi-Annual Amplitude and Phase^a Estimates for Prograde and Retrograde Oscillations in $$\hat{\chi }$$

	Prograde		Retrograde
	Amplitude [mas]	Phase [°]	Amplitude [mas]	Phase [°]
Annual oscillation
Observed	18.68 ± 0.23	−71.68 ± 0.70	11.80 ± 0.40	−123.35 ± 1.92
Atmosphere
ERA-Interim (IB)	17.91 ± 0.04	−85.65 ± 0.13	14.00 ± 0.06	−96.69 ± 0.25
Ocean
ECCOv4	3.97 ± 0.36	−25.03 ± 5.20	2.85 ± 0.40	−1.75 ± 7.98
GLORYS	2.50 ± 0.00	89.46 ± 0.08	4.31 ± 0.33	24.29 ± 4.43
ORAS	2.32 ± 0.37	17.01 ± 9.26	5.23 ± 0.24	52.04 ± 2.64
FOAM	2.53 ± 0.27	46.95 ± 6.03	1.99 ± 0.30	39.50 ± 8.68
Combination	2.12 ± 0.21	51.18 ± 5.66	3.81 ± 0.26	39.18 ± 3.88
Atmosphere + Ocean
ERA-Int. + ECCOv4	20.16 ± 0.16	−75.78 ± 0.47	14.04 ± 0.06	−85.03 ± 0.24
ERA-Int. + Combination	16.43 ± 0.10	−80.59 ± 0.36	11.57 ± 0.07	−83.43 ± 0.36
Atmosphere + Ocean + Hydrology/Ice
ERA-Int. + ECCOv4 + Hyd./Ice	23.22 ± 1.04	−55.38 ± 2.58	8.82 ± 1.11	−127.11 ± 7.20
ERA-Int. + Combination + Hyd./Ice	19.35 ± 0.95	−58.38 ± 2.81	6.49 ± 1.27	−134.63 ± 11.23
Semi-annual oscillation
Observed	4.37 ± 0.36	116.52 ± 4.69	4.75 ± 0.35	115.58 ± 4.17
Atmosphere
ERA-Interim (IB)	2.06 ± 0.63	9.98 ± 17.51	3.43 ± 0.36	124.75 ± 6.07
Ocean
ECCOv4	2.30 ± 0.30	137.57 ± 7.40	1.23 ± 0.02	−92.16 ± 0.70
GLORYS	2.50 ± 0.26	134.41 ± 5.97	0.35 ± 0.32	29.40 ± 52.41
ORAS	3.38 ± 0.23	125.28 ± 3.90	1.53 ± 0.34	−30.94 ± 12.79
FOAM	2.95 ± 0.34	149.63 ± 6.54	1.11 ± 0.08	−77.70 ± 4.31
Combination	2.86 ± 0.24	136.30 ± 4.89	0.86 ± 0.24	−44.27 ± 16.12
Atmosphere + Ocean
ERA-Int. + ECCOv4	1.94 ± 0.13	80.29 ± 3.89	2.56 ± 0.61	141.61 ± 13.71
ERA-Int. + Combination	2.33 ± 0.01	90.98 ± 0.30	2.60 ± 0.37	121.14 ± 8.09
Atmosphere + Ocean + Hydrology/Ice
ERA-Int. + ECCOv4 + Hyd./Ice	3.38 ± 1.75	127.06 ± 29.69	2.25 ± 0.22	85.65 ± 5.61
ERA-Int. + Combination + Hyd./Ice	3.87 ± 1.67	130.44 ± 24.73	2.63 ± 1.07	65.34 ± 23.36

^aIB is inverted barometer; reference date for phase is 1 January 2006, 12 UTC.

Acknowledgments

We are grateful for comments and suggestions provided by two anonymous reviewers. This work was supported by the German Research Foundation (DFG, Project no. 459392861). Open Access funding enabled and organized by Projekt DEAL.

Data Availability Statement

The data sets used in this study are available from the following links: Ocean reanalyses (https://resources.marine.copernicus.eu/product-detail/GLOBAL_REANALYSIS_PHY_001_031/INFORMATION), OAM from ECCOv4r3 (https://isdc.gfz-potsdam.de/ggfc-oceans/oam/), ETOPO1 (https://www.ncei.noaa.gov/access/metadata/landing-page/bin/iso?id=gov.noaa.ngdc.mgg.dem:316), GRACE JPL Release 06 mascons (https://podaac-tools.jpl.nasa.gov/drive/files/GeodeticsGravity/tellus/L3/mascon/RL06/JPL/v02/CRI/netcdf), and SPACE2018 rotation data (https://keof.jpl.nasa.gov/combinations/). All angular momentum estimates analysed in this study (except those from ECCOv4r3) are provided in Börger and Schindelegger (2022).