The Global S \(_1\) Tide in Earth’s Nutation

Surveys in Geophysics, Feb 2016

Diurnal S\(_1\) tidal oscillations in the coupled atmosphere–ocean system induce small perturbations of Earth’s prograde annual nutation, but matching geophysical model estimates of this Sun-synchronous rotation signal with the observed effect in geodetic Very Long Baseline Interferometry (VLBI) data has thus far been elusive. The present study assesses the problem from a geophysical model perspective, using four modern-day atmospheric assimilation systems and a consistently forced barotropic ocean model that dissipates its energy excess in the global abyssal ocean through a parameterized tidal conversion scheme. The use of contemporary meteorological data does, however, not guarantee accurate nutation estimates per se; two of the probed datasets produce atmosphere–ocean-driven S\(_1\) terms that deviate by more than 30 \(\upmu \)as (microarcseconds) from the VLBI-observed harmonic of \(-16.2+i113.4\) \(\upmu \)as. Partial deficiencies of these models in the diurnal band are also borne out by a validation of the air pressure tide against barometric in situ estimates as well as comparisons of simulated sea surface elevations with a global network of S\(_1\) tide gauge determinations. Credence is lent to the global S\(_1\) tide derived from the Modern-Era Retrospective Analysis for Research and Applications (MERRA) and the operational model of the European Centre for Medium-Range Weather Forecasts (ECMWF). When averaged over a temporal range of 2004 to 2013, their nutation contributions are estimated to be \(-8.0+i106.0\) \(\upmu \)as (MERRA) and \(-9.4+i121.8\) \(\upmu \)as (ECMWF operational), thus being virtually equivalent with the VLBI estimate. This remarkably close agreement will likely aid forthcoming nutation theories in their unambiguous a priori account of Earth’s prograde annual celestial motion.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs10712-016-9365-3.pdf

The Global S \(_1\) Tide in Earth’s Nutation

Surv Geophys The Global S1 Tide in Earth's Nutation Michael Schindelegger 0 1 2 3 4 David Einsˇpigel 0 1 2 3 4 David Salstein 0 1 2 3 4 Johannes Bo¨ hm 0 1 2 3 4 0 School of Cosmic Physics, Dublin Institute for Advanced Studies , Dublin , Ireland 1 Department of Geodesy and Geoinformation, TU Wien , Gußhausstraße 27-29, 1040 Vienna , Austria 2 Michael Schindelegger 3 Atmospheric and Environmental Research, Inc. , Lexington MA , USA 4 Department of Geophysics, Charles University in Prague , Prague , Czech Republic Diurnal S1 tidal oscillations in the coupled atmosphere-ocean system induce small perturbations of Earth's prograde annual nutation, but matching geophysical model estimates of this Sun-synchronous rotation signal with the observed effect in geodetic Very Long Baseline Interferometry (VLBI) data has thus far been elusive. The present study assesses the problem from a geophysical model perspective, using four modern-day atmospheric assimilation systems and a consistently forced barotropic ocean model that dissipates its energy excess in the global abyssal ocean through a parameterized tidal conversion scheme. The use of contemporary meteorological data does, however, not guarantee accurate nutation estimates per se; two of the probed datasets produce atmosphere-ocean-driven S1 terms that deviate by more than 30 las (microarcseconds) from the VLBI-observed harmonic of 16:2 þ i113:4 las. Partial deficiencies of these models in the diurnal band are also borne out by a validation of the air pressure tide against barometric in situ estimates as well as comparisons of simulated sea surface elevations with a global network of S1 tide gauge determinations. Credence is lent to the global S1 tide derived from the Modern-Era Retrospective Analysis for Research and Applications (MERRA) and the operational model of the European Centre for Medium-Range Weather Forecasts (ECMWF). When averaged over a temporal range of 2004 to 2013, their nutation contributions are estimated to be 8:0 þ i106:0 las (MERRA) and 9:4 þ i121:8 las (ECMWF operational), thus being virtually equivalent with the VLBI estimate. This remarkably close agreement will likely aid forthcoming nutation theories in their unambiguous a priori account of Earth's prograde annual celestial motion. 1 Introduction Nutation Geophysical excitation Atmospheric Describing variations of our planet’s orientation in space is a multidisciplinary subject matter that has occupied the attention of mathematicians, astronomers, and geophysicists alike. Nutations, that is, periodic motions of a predefined physical or conventional reference axis with respect to an inertial system, have been classically modeled in a Lagrangian or Hamiltonian framework (Woolard 1953; Kinoshita 1977) as the rigid Earth response to gravitational lunisolar torques. The estimates of these standard treatments are accurate to a few tenths of mas (milliarcseconds), but the advent of precise observational data as well as the pursuit of insights into Earth’s internal constitution has stimulated the development of non-rigid nutation theories for a realistic Earth (Jeffreys and Vicente 1957; Molodensky 1961; Sasao et al. 1980) . In 1980, the International Astronomical Union (IAU) adopted theoretical values of the lunisolar forced nutations for an elliptical, oceanless, elastic Earth with a fluid outer core and solid inner core (Wahr 1981) , though the geophysical approximations and omissions within that model were soon to become larger than the requirements posed by space geodetic techniques such as VLBI (Very Long Baseline Interferometry). A timely formation of a new nutation series, which has been the reference model since its endorsement by the IAU in 2000, is documented in Mathews et al. (2002) (MHB for short) and explicitly allows for mantle anelasticity, inner core dynamics, and non-hydrostatic equilibrium effects. Basic Earth parameters that govern the nutation response to lunisolar and planetary torques are constrained to their ‘‘best estimates’’ from a least-squares fit of the theoretical nutation expressions to VLBI results. This semi-analytical approach to modeling Earth’s nutation leaves residuals with observational data below 0.1 mas. Some effort has been devoted by MHB to properly account for the nutation perturbations associated with the Earth’s fluid layers. These contributions range from a few tens of las (microarcseconds) to 1 mas in amplitude and can be understood as the manifestations of (quasi-)diurnal atmosphere–ocean dynamics in the terrestrial frame. Their underlying excitation mechanisms are twofold, comprising (1) the daily cycle of solar heating and (2) the differential gravitational forces that directly act upon the atmosphere and ocean and produce global-scale waves known as tides. At the major diurnal tidal frequencies, the oceanic variability is almost exclusively driven by the gravitational influence with minute modulations related to the hydrodynamic response to atmospheric forcing. Satellite altimetry provides an accurate global record of these signals in the modern ocean and is also typically used to infer the associated oceanic angular momentum (OAM) variations of the largest diurnal tides, K1, P1, O1, and Q1. Early OAM determinations for these four waves (Chao et al. 1996) were adopted by MHB to predict the full OAM spectrum across the diurnal band and subsequently correct the equations of motion and anelasticity in the nutation theory. Tides in the atmosphere around the central diurnal period of S1 are capable of exciting small additional, seasonal nutation waves that exceed the statistical uncertainties of VLBIbased parameters on the order of 10 las (Dehant et al. 2003). The gravitational components of these oscillations are—to the extent they have not been implicitly accounted for in the MHB model—negligibly small (Bizouard and Lambert 2002) and in fact overshadowed by thermal tides due to periodic radiation and absorption processes (Chapman and Lindzen 1970) . Nutation contributions of some minor radiational constituents, such as the P1 tide, have been thoroughly addressed by MHB, yet the main S1 wave proved to pose some sort of conundrum to the authors. Whereas the VLBI data (1979.8–1999.11) distinctly testified to the existence of an S1 influence in the form of a prograde annual nutation residual, available geophysical model estimates (Bizouard et al. 1998; Yseboodt et al. 2002) were deemed unreliable, and thus no ‘‘theoretical’’ account of the effect was incorporated into the dynamical equations. To compensate for this mismatch, MHB subtracted from the VLBI spectrum an a priori S1 harmonic of somewhat more than 100 las and superimposed the very same signal as a post-fit correction term to the final nutation series. From a practical point of view, this approach is legitimate and in fact aided by the pronounced harmonic character of the unexplained nutation variability in the prograde annual band. Corrections to the MHB model, published as daily celestial pole offsets by the International Earth Rotation and Reference Systems Service (IERS), are generally below 20 las in the S1 band, also for more recent years not included in the original MHB analysis; cf. Fig. 1. However, in the spirit of a theory that should be ultimately free of empirical adjustments (Fedorov et al. 1980) and also unambiguous in its account of the various effects at the prograde annual frequency, it is still worthwhile to strive for an independent S1 estimate from geophysical fluid models. Studies of this subject matter are required to accommodate not only the atmospheric portion of the tide but also the substantial 24-h oceanic mass redistributions driven by S1ðpÞ, the diurnal pressure variations at the sea surface. Hence, both atmospheric and oceanic oscillations are closely interrelated aspects of the same ‘‘global S1 tide,’’ labeled as such by Ray and Egbert (2004) as well as in the context of the present work. Owing to its radiational origin, the oceanic S1 variability might be also perceived as an indirect influence of the atmosphere on the rotation of the solid Earth, and it has therefore been occasionally classified as a ‘‘non-tidal’’ phenomenon (Brzezi n´ski et al. 2004; de Viron et al. 2004) —a terminology that shall, however, not be used in the following. Investigations of the S1 effect in nutation on the basis of dynamically coupled atmosphere–ocean models have been pursued primarily by A. Brzezin´ ski and collaborators; cf. Fig. 1 Prograde annual signal amplitude in the IAU2000 celestial pole offsets (CPO) w.r.t. the MHB model. Estimates are 3year sliding window fits, and error bars indicate standard deviation (SD) in amplitude that have been propagated rigorously from the CPO errors. The MHB VLBI analysis, documented in Herring et al. (2002), involves data up to November 1999 (red line) and features a threefold SD of 21 las in the prograde annual band (dashed black line) 60 50 )s40 a µ ( de30 u lit p Am20 10 0 Brzezin´ ski et al. (2004), Brzezin´ ski (2011) and references therein. These studies employed different atmospheric analyses and ocean models both in a full 3D baroclinic formulation as well as 2D (constant density) barotropic versions that efficiently capture short-period hydrodynamic processes. The inferred atmosphere–ocean excitation terms of the prograde annual nutation vary substantially, though, both among each other and from the geodetic VLBI value, with deviations usually being larger than 50 las. While the VLBI estimate itself is possibly perturbed by other, imperfectly modeled seasonal effects, we must proceed on the assumption that the probed general circulation models were not appropriately designed for S1-related investigations. Another though brief assessment of the radiational tidal influence on nutation is given in de Viron et al. (2004) and their seemingly close agreement (\15 las) with the VLBI value has been acknowledged by Dehant and Mathews (2009) (Sect. 10.11, ibid.). We think, however, that the results of de Viron et al. (2004) are questionable and actually affected by an incorrect conversion of excitation values to periodic nutation terms. This deficiency is particularly evident for their tabulated atmospheric contributions (P1, S1, and w1), which are inconsistent with what has been documented for the very same atmospheric dataset by Koot and de Viron (2011) and Bizouard et al. (1998) , even if one makes allowance for differences in the utilized transfer functions and the analyzed time spans. Note, e.g., that over the period 1991–2002, Fig. 2 of Koot and de Viron (2011) suggests 75 las for the S1 out-of-phase term, while de Viron et al. (2004) specify a value of 38 las. Without further insight into the actual (corrected) S1 nutation predictions of these authors, we will employ Brzezin´ ski (2011) as a reference study by which our results can be measured. Building on the elucidations of these pilot investigations, the key objective of the present work is to provide an up-to-date treatment of the global S1 tidal effect in nutation and ultimately arrive at an explanation of MHB’s empirical prograde annual nutation term. The modern-day aspect of our effort resides in the use of four of the currently most advanced atmospheric assimilation systems, comprising three constant-model, retrospective analyses (so-called reanalyses) for a principal time span from 1994 to 2013, and a shorter, operational dataset (2004–2013) that stems from the near real-time weather analysis with a steadily improving model. This atmospheric portion of our study can be rightly understood as a continuation of similar earlier assessments (Bizouard et al. 1998; Yseboodt et al. 2002) , and as such, it is partly motivated by the good agreement of atmospheric nutation estimates from reanalyses that are essentially the precursors of the presently tested models (Koot and de Viron 2011). To ensure conformity with these investigations, nutation values for the minor solar constituents (w1, P1, p1, and /1) are tabulated, even though our prime focus is on the S1 tide throughout. A second key theme of this work is to numerically model the dynamic ocean response to diurnal atmospheric pressure forcing, which has been a fruitful geophysical industry over the last decade; cf. Ray and Egbert (2004) , Dobslaw and Thomas (2005), Ponte and Vinogradov (2007) , or Carre`re et al. (2012) . Mere superpositions of these modeling results to our atmospheric excitation terms are invalid, though, and a rigorous treatment of the geophysically driven prograde annual nutation requires deducing hydrodynamic S1 solutions and respective OAM values that are consistent with the utilized atmospheric datasets. We adopt the recent barotropic time-stepping model of Einsˇpigel and Martinec (2015), designated as DEBOT (Barotropic Ocean Tide model developed by D. Einsˇpigel), and implement the necessary modifications for the problem in hand. Specifically, to obtain S1 tidal solutions that are on par with Ray and Egbert (2004) , ocean self-attraction and loading effects (SAL, Ray 1998b) are accounted for in an iterative fashion and the overestimation of sea surface elevations in deep water is mitigated by a parameterized expression for the barotropic-to-baroclinic energy conversion over abyssal hills (Bell 1975; Jayne and St Laurent 2001). Moreover, forcing the same hydrodynamic model with different pressure tide solutions should go some way to reveal the dependence of the global/regional character of S1 (and its OAM values) on variations in the barometric input data. This is a subtle issue that has yet not been addressed by the oceanographic community. Our approach is a climatological one inasmuch as the present formulation of DEBOT only allows for a strictly harmonic pressure loading by the S1 air tide, even though the temporal variability of this forcing can be large (Ray 1998a) . Seasonal modulations of S1ðpÞ by 1 cpy (cycle per year) correspond to the P1 and K1 constituents, inducing small radiational ocean tides that are automatically included in altimetric solutions of P1 and K1 and are thus of no practical significance. Variations on inter-annual timescales (e.g., Vial et al. 1994) pose, however, a more delicate challenge, which can be partly resolved by working with decadal-scale S1 averages for the coupled atmosphere–ocean system. To determine a favorable, i.e., inter-annually ‘‘quiet’’ averaging period, we assess the S1 variability both in the integrated atmospheric nutation values and in surface pressure. Specifically, the analysis of S1ðpÞ is conceived as a validation of model pressure tides against ‘‘ground truth’’ estimates from 50 island- and buoy-based barometers. This comparison is in fact a vital (though limited) measure in deciphering the fine margins in quality among the different models regarding the diurnal cycle. Further observational constraints on our model-based investigations are supplied by S1 determinations at 56 coastal tide gauges, which, to some extent, echo the varying degree of reliability of the simulated tidal heights from each atmospheric dataset. The paper is organized as follows. Section 2 places the present and previous nutation studies in the context of an evolving collection of meteorological assimilation systems and describes the main characteristics of the four models utilized herein. Atmospheric excitation time series are computed and mapped to nutation amplitudes in Sect. 3, complemented by the validation of model pressure tides against in situ estimates from pelagic barometers. The ocean model and its hydrodynamic configuration are thoroughly discussed in Sect. 4, and we assess the quality of our forward simulations both from an angular momentum perspective as well as in a comparison to coastal tide gauges. Section 5 finally synthesizes atmospheric, oceanic, and VLBI nutation results to address the mismatch between theory and observation in the S1 band. 2 Meteorological Data for Nutation Studies Solar tides in the atmosphere are not of immediate relevance to operational or retrospective analyses, yet a largely realistic model account of these oscillations is guaranteed by the use of insolational forcing physics in combination with in situ and remotely sensed meteorological data. Reanalyses, created by various weather agencies on the basis of an unchanging assimilation scheme over decades, are usually credited with a realistic longterm variability that also modulates the tides and, by implication, nutation amplitudes. Their products, issued with an invariable spatial and temporal resolution, have thus become the preferred means to investigate atmospheric effects in nutation. Table 1, taken from Schindelegger et al. (2015) , summarizes some basic information of presently available reanalysis datasets, sorted by a rough generation index (Dee et al. 2015) that is thought to reflect the varying degree of sophistication in terms of model physics, resolution, and Nameb a The generation index (GI) and information about horizontal resolution and assimilation technique are taken from Dee et al. (2015), while the model vintage (i.e., the fixation date of the agency’s operational model) has been extracted from the reanalysis-specific reference articles and differs from Dee et al. (2015) in individual cases. Yet inaccessible models (e.g., MERRA-2) and twentieth-century reanalyses that assimilate surface observations only (e.g., Compo et al. 2011) are not tabulated. We have also omitted citations of reanalyses that are not examined in the frame of the present work b Abbreviations NCEP National Centers for Environmental Prediction, ECMWF European Centre for Medium-Range Weather Forecasts, ERA ECMWF Reanalysis, JRA-25/JRA-55 Japanese 25-year/55-year Reanalysis, JMA Japan Meteorological Agency, MERRA Modern-Era Retrospective Analysis for Research and Applications, GMAO Global Modeling and Assimilation Office, CFSR NCEP Climate Forecast System Reanalysis c Abbreviations 3DVar/4DVar 3D/4D Variational Assimilation, IAU Incremental Analysis Update assimilation technique. Refer to the caption of Table 1 for any model abbreviations used in the following. Previous assessments of atmosphere-driven nutations have relied heavily on NCEP’s first-generation reanalysis R1, whose physical formulation and relatively coarse resolution (2 –2.5 for surface and vertical parameters) date back to 1995. Bizouard et al. (1998) , Yseboodt et al. (2002) , and Brzezin´ ski et al. (2004) derived R1-related estimates at tidal frequencies for particular reanalysis periods, while Koot and de Viron (2011) additionally analyzed NCEP R2 and the second-generation ERA-40 model over a common time span from 1979 to 2002. As a sole modern-day reanalysis, ERA-Interim (henceforth ERA) has been subject to an evaluation of nutation signals by Brzezin´ ski (2011). Ignoring differences due to varying analysis periods, the S1 estimates of these studies exhibit a fair agreement, roughly at 20–30 las, though individual outliers exist. In an attempt to document the same level of agreement or even further convergence for third-generation reanalyses, the present work derives nutation values for ERA and its contemporaries MERRA (Rienecker et al. 2011) and CFSR (Saha et al. 2010) . Operational models, designed for weather prediction on a daily basis, may be thought to be less suitable for long-term nutation studies. Indeed, S1 estimates of Yseboodt et al. (2002) from early operational systems of ECMWF, NCEP, and JMA differ by as much as 70 las, suggesting that the regular changes to the model and assimilation technique within each agency are capable of introducing artificial tidal variability (Koot and de Viron 2011). To some extent, though, the results of Yseboodt et al. (2002) are affected by time series limitations, e.g., data gaps of up to 30 % or occasionally short record lengths (3 years). If provided continuously, operational datasets might, in fact, still be proper vehicles for tidal studies, as the underlying analysis systems are optimized to represent the atmospheric variability on short timescales and as they are also readily adaptable to the introduction of new data types. By contrast, reanalyses assimilate an evolving observation record through a predefined framework and are thus prone to spurious variabilities (Dee et al. 2015). Moreover, and in the context of nutation, most model updates involve resolution and orography changes that possibly cause local discontinuities in continental surface pressure but are of no immediate consequence for all other components of the global S1 tide (i.e., the wind signal at higher altitudes, the pressure tide over the ocean, and, thus, the full oceanic S1 variability). Table 2 summarizes the main specifications of the globally gridded datasets from MERRA, CFSR, ERA, and the ECMWF operational model, denoted as EC-OP in the following. The analysis was initially conceived for the period of 2004–2013 (Schindelegger et al. 2015) but extended in retrospect to a 20-year window (1994–2013) for the three reanalyses. CFSR constitutes an exception, having been produced as a genuine, constant-model reanalysis until the end of 2010 with subsequent operational extensions that were, however, disregarded in the frame of the present work. We extracted standard 6-h analysis fields from the respective data archives for both ECMWF models as well as CFSR, whereas 3-h analysis/forecast combinations, designated as assimilated states, were utilized for MERRA. This high-resolution dataset has been a source of continued Earth rotation research at TU Wien, comprising also investigations of the semidiurnal atmospheric tide S2, which is not properly resolved by four-times-daily analysis products. To prepare for further subsequent work on S2, we also interlaced 3-h forecast data to the CFSR analysis fields after verifying that their inclusion had no evident effect on the S1 signature in surface pressure and nutation estimates. Atmospheric excitation is classically inferred from the two components of AAM (atmospheric angular momentum), comprising effects both due to particle movement (wind or motion term) and redistribution of matter (pressure or mass term). The evaluation of the wind term involves vertical integration over an appropriate number ([15) of isobaric levels and is thus computationally intensive. Yet, the higher-altitude horizontal convection associated with S1 constitutes a robust, large-scale signal that is integrated with sufficient accuracy from comparatively coarse mesh sizes of about 2 ; cf. Table 2. Given the pronounced small-scale (and possibly subgrid-scale) characteristics of the diurnal surface pressure variability over landmasses (e.g., Li et al. 2009), the mass component of AAM is preferably deduced from better resolved surface grids, although limitations are imposed by the intrinsic model resolution and our hardware resources. MERRA’s assimilated states of surface pressure are distributed at 1.25 in latitude and longitude, whereas 0.5 grids could be utilized for CFSR and ERA. Note that the native 80-km spacing of ERA (Table 1) is somewhat coarser than 0.5 , suggesting that these data were interpolated during the assignment process. By contrast, the chosen 1 mesh for the operational data is a largely downsampled version of the model’s fine intrinsic discretization, having improved in resolution from 40 km in 2004 to 16 km at the end of 2013. A preliminary comparison of two reanalyses with regard to their diurnal cycle is shown in Fig. 2 in the form of cotidal S1ðpÞ charts for MERRA and CFSR. Mean amplitude and phase lag values were obtained from a standard least-squares tidal analysis of 10-year pressure time series at each grid point location. These climatologies echo the well-known spatial characteristics of the diurnal barometric tide (see Ray and Ponte 2003, and references therein) but also exemplify that its representation in global analysis models can diverge. CFSR suggests higher S1 amplitudes almost throughout the world, particularly over landmasses in latitudes lower than 30 and in valleys for which the local diurnal oscillation is not resolved by MERRA (e.g., Sierra Nevada). Large regional-scale differences over flatter terrain (e.g., Sahara, Central Africa, India) portend to difficulties in the 240 210180150 120 models to represent the significant diurnal boundary-layer effects driven by sensible and latent heating from the ground (Dai and Wang 1999) . These non-migrating components can be excluded from the tidal spectrum by performing a Fourier decomposition of S1ðpÞ by wavenumber s (Chapman and Lindzen 1970) and retaining only the main Sun-synchronous, migrating S11 tide (s ¼ 1). Forced by tropospheric absorption processes, S11 corresponds to a longitudinally uniform wave that is clearly evident in Fig. 2 over oceanic areas. Latitudinal profiles of this migrating tide from MERRA and CFSR exhibit a fair agreement (not shown), although equatorial peak amplitudes differ slightly (64.1 Pa for MERRA, 67.2 Pa for CFSR) and an overestimation of about 10 Pa at latitudes close to 60 S can be observed for the CFSR solution. These discrepancies in pressure are likely to have a bearing on the simulation of the oceanic S1 tide. 3 Atmospheric Contributions to Nutation 3.1 Implementation Computation of the vector equatorial AAM mass term H~p ¼ Hxp þ iHyp (complex notation) involves weighted-area double integrals of surface pressure, whereas for the motion term H~w ¼ Hxw þ iHyw a full 3D summation of geometrically weighted horizontal winds is required. We applied the respective standard formulas, given, e.g., in Sect. 2.5 of Schindelegger et al. (2013) on the gridded datasets of Table 2, with lower boundaries in the vertical integration taken from the model-specific topographies. A priori corrections to the mass terms for an isostatic (inverted barometer, IB) ocean response to air pressure variations were categorically avoided, as the oceanic S1 tide is a dynamic phenomenon that will be rigorously estimated in Sect. 4. H~p;w are the basic excitation quantities that can be related to nutation through a proper dynamical theory. Sasao and Wahr (1981) devised corresponding expressions for the geophysically driven nutation from the angular momentum balance equations of a coupled two-layer Earth, comprising mantle and a fluid core which are allowed to deform elastically under the action of body tides and atmospheric (oceanic) loads at the Earth’s surface. Brzezin´ ski (1994) reformulated this pilot equation to a practicable broad-band excitation scheme for both nutation and polar motion. Yet, for reasons of consistency, the comparison of geophysical model estimates with the S1 post-fit correction terms of MHB and Koot et al. (2010) is better accomplished through the excitation scheme of Koot and de Viron (2011). Their formalism conforms to MHB’s nutation theory for an up-to-date Earth model with inner core dynamics and anelastic properties. Perturbations n~ r ð Þ ¼ dX þ idY of the celestial pole offsets in X and Y in response to changes of AAM at some Earth-referred, retrograde diurnal frequency r (in cycles per sidereal day, cpsd) are modeled as n~ r ð Þ ¼ ~p r H~0pðrÞ T ð Þ XðC AÞ H~0wðrÞ T~wðrÞ XðC AÞ where X is the nominal sidereal angular velocity, A and C denote the equatorial and polar principal moments of inertia of an axisymmetric solid Earth (i.e., mantle and crust), and T~p;wðrÞ are transfer functions describing Earth’s nutation response to atmospheric forcing as conveyed by the periodic terms H~0p;wðrÞ, which are defined below. The transfer functions read ð1Þ T~p;w r ð Þ ¼ X4 N~ip;w i¼1 r r~i comprising resonances at the frequencies r~i (cpsd) of the four rotational normal modes of a three-layer Earth: the Chandler wobble (CW), the free core nutation (FCN), the free inner core nutation, and the inner core wobble (Koot and de Viron 2011). The strengths of these resonances upon mass and motion excitation are characterized by the coefficients N~ip;w, specified in Table 3 with a truncation at i ¼ 2 that retains CW and FCN and excludes the inner core modes without loss of accuracy. If viewed from the surface of the rotating Earth, the FCN occurs as retrograde nearly diurnal oscillation, thus providing significant enhancement to excitation effects associated with atmosphere–ocean dynamics at S1 and adjacent tidal lines. Note, however, that nutations are much more efficiently driven by the mass term than by relative particle motion (Sasao and Wahr 1981; Brzezin´ ski 1994) ; cf. also the excess of N~p at the FCN frequency relative to N~w by a factor of 200 (Table 3). The forcing terms H~0p;wðrÞ in Eq. (1) are complex coefficients of time dependence / eirt as seen from the rotating reference frame (Koot and de Viron 2011). Yet, the amplitudes of these sinusoids are typically estimated in inertial space after translating the terrestrial AAM time series H~ðtÞ to their celestial counterparts H~0ðtÞ via the demodulation (Brzezin´ ski 1994) that is applicable to both pressure and wind effects (respective superscripts have been omitted for brevity). The exponent represents a sufficiently accurate linear approximation for the Greenwich sidereal time (Bizouard et al. 1998) , employing a phase offset of U0 referred to t0 at 12 h UT1, 1 January 2000 (J2000.0). The conventional expression for the Earth rotation angle of the IERS (Petit and Luzum 2010) is compatible with Eq. (3) and implies U0 ¼ 0:7790572732640 rad as well as X ¼ ð2prÞ rad per solar day, where r ¼ 1:00273781191135448 scales solar to sidereal time intervals; see also Koot and de Viron (2011). The demodulation procedure preserves amplitudes but maps the retrograde (r\0) quasi-diurnal spectral components to low frequencies, with the center frequency X (i.e., the K1 band) shifted to zero and S1 appearing as prograde annual line in the celestial frame. Dominant (intra-)seasonal signals in the original AAM series are mapped to high frequencies in space and are efficiently removed through filtering (Bizouard et al. 1998) . To this end, we applied an idealized rectangle filter with cutoff at 20 cpy (cycles per year) on a Theoretical, complex-valued frequencies of CW and FCN are given in the Earth-fixed frame and correspond to a terrestrial period of 396.06 solar days with quality factor Q ¼ 223 for the CW and a celestial period of 429:05 solar days with a (terrestrial) quality factor of Q ¼ 19736 for the FCN resonance the frequency transform of H~0ðtÞ and resampled the proper inverse transform in the time domain at daily intervals. Experiments with more customary time domain filters and a range of reasonable cutoff frequencies testified to the insensitivity of our nutation results to details in the filtering strategy. To convert the low-frequency, quasi-harmonic celestial AAM time variability associated with S1 and its seasonal modulations to periodic circular components H~0ðrÞ for use in Eq. (1) we imposed a Fourier decomposition ð4Þ on the complex-valued filter output. The non-dimensional frequencies fmjgj5¼1 of the demodulated tidal constituents fS1; w1; P1; /1; p1g are f1; 1; 2; 2; 3g=366:26 (Koot and 5 de Viron 2011), and respective phase values fujgj¼1 agree with those from the corresponding lunisolar nutation terms. The uj are readily computed from the fundamental arguments of nutation theory (Petit and Luzum 2010) , using the integer multipliers specified in Table 4. A standard least-squares fit of H~0ðtÞ onto these basis functions provides the unknown parameters a~j composed of real (in-phase, ip) and imaginary (out-ofphase, op) parts that form the complex-valued forcing terms in Eq. (1) at discrete terrestrial frequencies rj ¼ mj 1 (cpsd). The constant pole offset contribution from the K1 tide, conveyed by the estimated c~ term, was excluded from further consideration, as were the minute secular (precession) contributions associated with the time derivatives of nutation arguments. The harmonic decomposition in Eq. (4) exactly follows the model of Koot and de Viron (2011) except for our inclusion of the /1 component at 2 cpy that suggests some interannual modulation of the thermal S1 tide. While the existence of such a modulation is debatable (Bizouard et al. 1998) , its impact on nutation is at the level of 10 las and thus comparable to the usually modeled p1 term. Moreover, for numerical reasons, we performed the least-squares fit on the basis of prescaled AAM time series H~0ðtÞ=ðXðC AÞÞ (Eq. 1) in units of las, similar to the classical ‘‘celestial effective angular momentum functions’’ of Brzezin´ ski (1994). The excitation scheme of this author was tested briefly and found to yield nutation results well within 5 % of the estimates from the above formalism. a Periods and phases (referred to J2000.0) result from the linear combinations of series expansions as given by Petit and Luzum (2010) for each argument. The table is identical to Table 1 of Koot and de Viron (2011) except for the /1 term Period (solar days) 3.2 Results Figure 3 displays the superimposed pressure and wind nutation estimates n~ðrÞ in the prograde annual band, both as mean contributions over the model-specific time spans as well as yearly values obtained from repetitions of the analysis in Sect. 3.1 with a 3-year sliding window. Consistent with illustrations in Yseboodt et al. (2002) and Koot and de Viron (2011), little agreement is seen between all S1 curves in terms of their inter-annual variability, even though the post-2004 estimates are stable within 20 las for each model. Roughly 75 % of the observed fluctuations are driven by the pressure term and likely relate to random perturbations of the second-order tesseral harmonic in surface pressure, i.e., the only component of the S1ð pÞ wave that efficiently excites Earth’s nutational motion. Amplitudes of this mode do not exceed 10 Pa, so its representation in atmospheric assimilation systems is prone to noise interferences. By contrast, a signal with a possibly physical origin is evident for MERRA during 1997–2001, coinciding in time with a peak El Nin˜ o event in 1997/1998 and subsequent cold La Nin˜ a conditions up to 2001; cf., e.g., the Oceanic Nin˜ o Index tabulated at http:// www.cpc.noaa.gov/products/analysis_monitoring/ensostuff/ensoyears.shtml (accessed 29 September 2015). The warm phase of ENSO (El Nin˜ o–Southern Oscillation) has been previously suggested to alter the radiative forcing of the solar tide in the troposphere (Lieberman et al. 2007), thereby providing significant enhancement to diurnal pressure oscillations across the Pacific (Vial et al. 1994) . The response of the climate system to 80 ) s aµ 60 ( e s a h p f− 40 o − t u o 20 0 S 1 −80 2009 ENSO events is, however, not restricted to the Tropics but can entail atmospheric circulation changes in higher latitudes that may ultimately couple to nutation; see similar conjectures in Yseboodt et al. (2002) . Assessing whether the irregular nutation changes from MERRA in Fig. 3 are linked to ENSO or merely represent spurious variabilities in the wake of observing system changes (Robertson et al. 2011) is beyond the scope of this study, though. We will in fact avoid these signals in our selection of the mean analysis window below. In terms of multi-year nutation averages, the ECMWF-based solutions agree particularly well with each other and with Koot and de Viron (2011)’s results for the firstgeneration NCEP models from 1979 to July 2002. The ERA-40 estimate of these authors is anomalous in the ip component ( 58:2 las) due to erratic S1 variations up to the mid1990s (cf. Fig. 2 of Koot and de Viron 2011) . Disregarding the impact of the ‘‘ENSO swerve’’ on the MERRA solution, the only nutation anomaly in the present work is a large CFSR estimate, with op pressure term values (50 las as from 1998) exceeding the predictions from other reanalyses by 20–30 las. This overestimation traces back to the dubious CFSR pressure oscillations of more than 40 Pa in the Southern Ocean (Fig. 2), a region that is void of conventional in situ observations and sensitive to the details of radiance data assimilation. Note also that the quality with which atmospheric tides can be represented in analysis systems is tied to the time step of radiative processes (Poli et al. 2013) . Trading off computational costs and a fine spatial resolution, CFSR integrates its longwave radiation parameterization every 3 h (Saha et al. 2010) , significantly coarser than the hourly time step recommended by Poli et al. (2013) and employed within MERRA, ERA, and EC-OP. Numerical results for our nutation analysis are reported in Table 5, using—with some exceptions—an averaging period from 2004 to 2013 that has been specified after consulting station tide determinations in next section. If the somewhat deficient CFSR results are discarded, S1 estimates from third-generation reanalyses and EC-OP deviate from each other by less than 22 las, which slightly betters the agreement noted by Koot and de Viron (2011) and conforms with the threefold VLBI SD in the prograde annual band (Fig. 1). We have also assessed the stability of the S1 peak in terms of its frequency through a Morlet wavelet analysis of demodulated filter residuals H~0. Deviations from the nominal S1 ridge at 365.26 days (solar days) are well within 5 days for all datasets except for MERRA, which exhibits a transition from 355 days in 1998 to 385 days in 2002 before leveling off exactly at the annual period (not shown). These minor fluctuations contrast with the 30 day range deduced by Dehant et al. (2003) on the basis of NCEP R1 data during 1958–1999. We surmise, however, that the estimate of Dehant et al. (2003) is less reliable due to the inclusion of reanalysis products prior to 1979, i.e., a period that lacks both satellite retrievals and a broad network of in situ pressure observations in the southern hemisphere. Our nutation results for the minor solar constituents (w1, P1, /1, p1) can be compared with estimates tabulated in Bizouard et al. (1998) , Yseboodt et al. (2002) , Brzezin´ ski et al. (2004), or Koot and de Viron (2011). Here, we only point out that the harmonics fitted to the four atmospheric datasets agree well for P1 and p1 but differ substantially in the w1 band, which is of interest for studies of the FCN and Earth’s internal properties (Dehant and Defraigne 1997; Koot and de Viron 2011). Temporal variations of this tide can be large ( [ 100 las for individual models), and periods of its wavelet ridge vary within 30 days, probably as a reflection of a strong stochastic atmospheric influence. Nonetheless, a comparatively good inter-model agreement is found for the ip component of w1 ( 50 las). f o d o i r e p s i s y l a n a s trem tlao T ip d m in o rf W ip n o it e a r t u u s n s e o rP ip t s n o i t u b i r t n o c l e c d i r o e h M p s o m t a c i d o i r ea P 3 1 5 0 2 e – lab 004 rem T 2 T 5 9 7 7 7 7 2 2 1 4 2 2 2 6 4 9 9 0 3 8 p 4 4 8 5 6 4 5 2 2 4 4 4 4 4 1 1 o 0 5 3 3 3 3 0 0 2 1 5 5 1 3 3 1 4 5 op 27 26 3 2 2 4 4 3 4 4 9: 8: 3: 8: 9: 1: 1: 2: 3: 8: 1: 3: 9: 3: 0: 1: 8: 3: 2: 8: 5 0 2 1 1 0 5 8 9 0 0 0 6 7 7 2 1 7 2 1 2 3 3 3 2 3 4 5 7 5 1 1 1 1 1 1 b b b b A A A A A A A A R R b P R R b P R R b P R R b P þ ( n o r i V t o o K f o 2: 3: 2: 2: 2: n 0 0 0 0 0 io t n e 4: 6: 7: 3: 1: v p 1 1 0 0 0 n o o c d 2: 3: 2: 2: 2: n 0 0 0 0 0 a 9: 3: 1: 1: 1: e l 1 2 2 1 1 b a 2: 2: 2: 2: 2: f 0 0 0 0 0 o s t n 1: 0: 3: 8: 7: e 4 4 3 2 2 m u g r a m a d 2: 2: 2: 2: 2: n 0 0 0 0 0 u f e 6: 3: 1: 5: 2: th 2 2 2 2 2 rP ip d e r 0 r e 1 f 0 e 2 r – e 4 r 0 a 0 b s 2 A A t l R R b P n d e e o d R R SR A -O n ir o E E F R C o e p p M M M C E E ) r a e y 3 = 1 þ ( 1 p n A I 3.3 Validation of Pressure Tides against In Situ Data 30°N 0° 30°S (b) 30°N 0° 30°S 60°N 60°S 60°N 60°S Stations Buoys 0 60 120 180 240 300 360 semidiurnal L2 tide. Again, care was exercised in selecting stations at sufficiently small islands and atolls. The final subset of estimates, also derived by Schindelegger and Dobslaw (2016) , comprises 16 buoy locations that are part of the Tropical Moored Buoy System (McPhaden et al. 2010) . Further densifications were attempted but led to clusters and subsequent biases in the statistics given below. The median time series length is 6 years, the maximum is 26 years (Bermuda), and short time spans of only 2 years occurred for eight sites, most of them being buoys. Atmospheric model pressure values were evaluated by bilinear interpolation at the locations of the 50 ground truth stations and tidally analyzed in a moving 3-year window. RMS statistics and globally averaged amplitude differences from the comparison of these windowed S1 solutions to the in situ estimates (non-windowed) are shown in Fig. 5, both for the full global network as well as for a subset of 20 stations excluding latitudes lower than 18 . This restriction avoids an over-emphasis on the large migrating pressure tide near the equator and tests smaller signals in mid-latitudes, i.e., regions of increased importance for nutation. The resulting network (Fig. 4) is, however, very sparse and dominated by 13 stations in the Pacific. 10 (a) (c) ) a P (h 9 t u r td 8 n u rgo 7 o t S 6 M R 5 5 )a 4 P ( s 3 e c en 2 r e iff 1 d de 0 u it lp−1 m A−2 2005 Year 1995 2000 2010 1995 2005 Year 2010 −13995 2000 2005 Year 2010 −13995 2000 2005 Year 2010 10 (b) ) a P (h 9 t u r td 8 n u rgo 7 o t S 6 M R 5 5 )a 4 P (s 3 e c en 2 r iffe 1 d de 0 u t lip−1 m A−2 (d) MERRA CFSR ERA EC−OP 2000 In Fig. 5a (full compilation), ERA features the largest RMS misfits at about 9 Pa, accumulated through a significant overestimation of S1ð pÞ in the Tropics and an underestimation of the tidal amplitude elsewhere (Fig. 5d). These deficiencies can be expected to lead to a less reliable oceanic S1 tide. For equatorial Pacific stations, the second ECMWF dataset EC-OP also produces an excess in amplitude (by about 10 Pa) that maps into Fig. 5c but has no effect on the reduced network (Fig. 5d). Yet, median and normalized RMS differences are generally better for EC-OP than for the probed reanalyses and highlight the accuracy of the operational solution in a global domain. The CFSR statistics are comparable to those of MERRA and EC-OP, evidently unaffected by the southern hemispheric pressure anomalies (Sect. 3.2) as these regions are not sampled by our ground truth network. The 1998 El Nin˜ o and its subsequent reversal during 1999–2001/2002 introduce irregular tidal behavior and larger RMS values in all three reanalyses when validated against the climatological in situ solution. Following our prescription of minimal interannual S1 variability, we excluded model data up to 2002 from further consideration. A sufficiently long averaging period, required to somewhat conform with the 20-year mean S1 fit of MHB, might thus be realized by the overlapping 7-year window (2004–2010) common to all models. Disregarding the partially deficient CFSR analysis, we finally adopted the time span from 2004 to 2013 as the main analysis window for MERRA, ERA, and EC-OP, while CFSR results were averaged through 2004–2010, and MERRA excitation data for 2004–2010 were maintained as well, but only as a secondary option. This duality for MERRA is motivated by the deterioration of RMS values (Fig. 5a) and the moderate drifts in nutation estimates (Fig. 3a) as from the year 2010. A brief numerical comparison of the resulting model pressure tide climatologies against the 50-station set of S1 estimates is given in Table 6. Median absolute differences (MAD) are included as a more robust supplement to the averaged RMS misfits, underlining the reliability of all analysis models other than ERA. The general level of consistency with ground truth data (5–6 Pa) closely resembles values obtained by Ray and Ponte (2003) , but note that these authors have additionally applied a small phase shift Du to their model estimates. Median phase lag differences for MERRA and CFSR in Table 6 generally confirm Ray and Ponte (2003) ’s assumption of Du ¼ 5 , and imposing the corresponding correction on the MERRA tide does indeed reduce the RMS misfit with station data to 4.9 Pa. We have, however, abstained from revising the tidal phases to avoid inconsistencies with the diurnal cycle in AAM. a Model tides are 2004–2013 averages except for CFSR (2004–2010) 4 Numerical Modeling of the Oceanic S1 Tide 4.1 Ocean Model Configuration DEBOT (Einsˇpigel and Martinec 2015) is a recently developed ephemeris-forced, barotropic time-stepping model conceived to study the effect of the ocean flow on Earth’s magnetic field. In the present work, we create a spin-off of the model’s hydrodynamic core for individual partial tides, with a rigorous treatment of SAL and a parameterized drag term to account for the conversion from barotropic waves to baroclinic internal tides (IT) over rough bottom topography. The one-layer shallow water momentum and mass conservation equations define the horizontal velocity vector u and the tidal surface displacement f ou ot þ f z^ u ¼ gr f CDkuku H þ f fEQ fSAL fMEM P=gq CIT u H þ f þ AH r r of ot ¼ r ½ðH þ fÞu ð5Þ ð6Þ where f is the Coriolis parameter oriented along the local vertical unit vector z^, g is the nominal gravitational acceleration, r signifies the spherical del operator, H is the resting water depth, q is the average density of seawater, CD ¼ 0:003 denotes a dimensionless drag coefficient in the standard expression for quadratic bottom friction, and CIT is a location-dependent scalar (in units of m s 1) to represent the drag due to tidal conversion. The forcing terms in the gradient operator of Eq. (5) comprise the gravitational equilibrium tide fEQ, a combination of self-attraction/loading and ‘‘memory’’ elevations fSAL and fMEM to realize the SAL scheme of Arbic et al. (2004) , as well as the atmospheric pressure tide P ¼ S1ð pÞ. AH r r is a comparatively rigorous implementation of the horizontal turbulent eddy viscosity with a second-order tensor r related to the Reynolds stress tensor; see Einsˇpigel and Martinec (2015) for details. Here, we keep this term to eschew possible numerical instabilities in our medium-resolution runs (Egbert et al. 2004) , with horizontal viscosity AH set to the widely cited value of 103 m2 s 1. Our forward tidal solutions and OAM results are insensitive to the exact value of AH , unless it is used as a tuning parameter of inordinately large magnitude ( 105 m2 s 1); cf. Arbic et al. (2004) . Equations (5) and (6) were solved by finite difference time-stepping on a 1=3 C-grid, covering the latitude range from 78 S to 78 N with rigid walls assumed at the top and the bottom of the domain. This setting does not allow for accurate tidal modeling in the Weddell and Ross Sea, or in the Arctic Ocean. Yet, we readily accept such high-latitude limitations given our interest in the equatorial component of Earth’s rotation that has a peak sensitivity at 45 from the equator. The bottom topography was derived from the bedrock version of the fully global 10 10 ETOPO1 database (Amante and Eakins 2009) by choosing average values over each 1=3 model grid cell and setting depths between 10-m and the 0-m land–sea boundary to 10 m. Coastlines in the Antarctic come from a recent data-assimilative ocean model (Taguchi et al. 2014) and are similar to those of Padman et al. (2002) , with the cavities under the floating ice shelves considered as part of the ocean domain; cf. also Arbic et al. (2004) or Carre`re et al. (2012) . Blocking ice shelf areas as dry cells would lead to a noticeable increase of the tidal variability in southern hemisphere waters, also amplifying the oceanic contribution to the prograde annual nutation by roughly 10 las in both ip and op components. However, in these simulations, also the RMS misfit of the gravitationally forced constituents (M2 and O1; see below) to altimetry-based reference solutions, (FES2012, Carre`re et al. 2012) increases, consistent with similar control runs by Wilmes and Green (2014) . We thus proceed on the assumption that vertically displaceable ice shelves allow for a more realistic account of the tides, including S1. Other aspects of the DEBOT configuration closely follow Ray and Egbert (2004) . We prescribe equilibrium tidal forcing (fEQ) for M2 and O1 with amplitudes and solid Earth tide corrections taken from Table 1 of Arbic et al. (2004) . The resulting largermagnitude background variability appears to aid the fidelity with which S1 can be simulated in various basins and bays, but an extension to more than one diurnal and semidiurnal gravitational constituent is expendable as it alters our S1 OAM estimates by less than 3 %. Experiments with different equilibration periods (up to 90 days) showed that the spin-up time of the model could be reduced to 12 days with little effect (cf. Arbic et al. 2004) , although in very shallow waters (Gulf of Thailand, Java Sea) the convergence of S1 takes considerably more time than that of any gravitational tide, presumably due to the vagaries of the pressure forcing near landmasses (Fig. 2). With 12 days reserved for equilibration, we integrated the model in each of our runs for 40 days at a time step of 24 s, harmonically analyzing the last 28 days to deduce the tidal constants of S1 (as well as M2 and O1) in terms of sea level elevation and barotropic volume transports uH. 4.2 Effects of Self-Attraction and Loading (SAL) Gravitational self-attraction and yielding of the solid Earth to the weight of the water column (Hendershott 1972) are feedback effects to the tidal dynamics and included in Eq. (5) as an additional equilibrium-like tide fSAL. This term can be related to the (unknown) local tidal elevation f through convolution with the global SAL Green’s function G (Ray 1998b) fSALð/; kÞ ¼ qa2 ZZ fð/0; k0ÞGðwÞsin/0d/0dk0 ð7Þ where a is the Earth’s radius and w measures the angular separation of ð/; kÞ from the load with spherical coordinates ð/0; k0Þ. For our 1=3 model, values of GðwÞ were interpolated from the SAL kernel function tabulated in Stepanov and Hughes (2004) . Explicit usage of Eq. (7) in the momentum equations is computationally unfeasible (Egbert et al. 2004) , so alternative implementation schemes are required. To first order, the full convolution with G is approximated by a simple scalar multiplication fSAL bf (Accad and Pekeris 1978) , with b usually taken to be in the range of about 0.08 to 0.12. This widely used approximation is inappropriate for all locations in the ocean (Ray 1998b) and accurate tidal modeling necessitates a more rigorous handling of the effect. In tide models forced by a suite of individual constituents, the unparameterized formalism of Eq. (7) can be applied in a comparatively simple manner via iteration, that is, repeated model runs where each simulation employs a better approximated SAL term to gradually achieve convergence between the tidal elevations and fSAL. We applied the iteration method of Arbic et al. (2004), initialized by the scalar SAL estimate using a nominal value of b ¼ 0:12 that is an appropriate choice for diurnal tides; cf. Parke (1982) and Fig. 11 of Einsˇpigel and Martinec (2015). Once this initial run is completed and harmonically analyzed, the tidal components (sine and cosine terms) of M2, O1, and S1 are inserted into Eq. (7) in an intermediate offline a Tidal solutions have been computed from our hydrodynamic model using atmospheric pressure forcing from Ray and Egbert (2004) . Amplitudes are in units of 1023 kg m2 s 1 and cotidal phases are given relative to Greenwich noon, consistent with the Doodson convention for the S1 phase as given in Ray and Egbert (2004) . For the respective OAM formulas, refer to Chao et al. (1996) computation to derive a first solution of fSAL for each tidal constituent. The following simulation then time steps the sum of all partial SAL tides as well as an additional memory term (Arbic et al. 2004) fMEM ¼ bðf fPREV Þ ð8Þ that measures the departure of the tidal height f in the current (second) run from the cumulative M2-O1-S1 elevation fPREV in the previous (first) run. Subsequent iterations are performed in the same manner, drawing on continuously updated maps of fSAL and fMEM . The correction term in Eq. (8) guarantees rapid convergence of the SAL scheme, as exemplified by diminishing RMS discrepancies of the gravitational constituents against FES2012 tides in successive simulations. Specifically, with the choice of b optimized for the diurnal band, our forward solutions of O1 remain effectively unchanged after the first iteration, whereas sufficient accuracy for semidiurnal tides is reached after three iterations. More to the point, rapid equilibration of tidal dynamics is also observed for the radiational S1 tide. Table 7 presents successively updated OAM mass values of S1 as obtained from a three-times iterative DEBOT run with the pressure forcing S1ð pÞ taken from Ray and Egbert (2004) and IT drag (cf. next section) switched off. For equatorial components in particular, the scalar approximation appears to provide reasonably accurate initial OAM estimates, deviating by no more than 5 in phase and less than 10 % in amplitude from the (arguably) self-consistent third iteration. Yet, the scalar SAL relation is inadequate for both the axial OAM component and the comparison of simulated S1 surface elevations to coastal tide gauges (Sect. 4.4). Accordingly, results from all of our forward runs presented below have been inferred after completing the second model iteration. 4.3 Internal Tide (IT) Drag Scheme Consistent with previous studies of forward-modeled barotropic tides (Jayne and St Laurent 2001; Arbic et al. 2004; Egbert et al. 2004) , surface elevations and tidal energies are poorly represented in DEBOT unless allowance is made for the substantial amount of drag generated by internal tides over major bathymetric features. With this proper dissipation mechanism omitted, area-weighted RMS differences Df1 to the FES2012 reference tide f~R 1 Computed as Df ¼ sffiRffiRffiffiffijffif~ffiffiffiffif~ffiffiRffijffiffi2ffidffiffiAffi, equivalent to the time-averaged expression of Arbic et al. (2004), with 2RR dA grid points poleward of 66 and waters shallower than 1000 m excluded. (complex sinusoid) are as large as 14.2 cm and 3.0 cm for M2 and O1, respectively. These values translate to a mere 72 and 79 % of sea surface height variance explained. Moreover, S1 charts deduced from IT-free simulations display a number of apparent regional artifacts, such as persistently high amplitudes of the tide in the northern Atlantic ( 1 cm) or the South China Sea ( 2 cm) that have no correspondence in both the altimetric and hydrodynamic S1 solutions of Ray and Egbert (2004) . Parts of the OAM discrepancy of our initial control run (Table 7) to Ray and Egbert (2004) ’s benchmark values can be understood in this light. To increase the fidelity of our model tides and in particular S1, we implemented the linear tidal conversion formulation of Zaron and Egbert (2006) as described and slightly modified by Green and Nycander (2013). In this parameterization, the local drag coefficient is explicitly proportional to the slope of the scattering topography where C ¼ 50 is a non-dimensional constant, x denotes the frequency of the tidal motion, and theoretical buoyancy frequencies N follow from the prescription of a horizontally uniform abyssal stratification. Values of N at the ocean bottom (Nb) as well as vertical averages (N) over the entire water column are calculated from Green and Nycander (2013) CIT ¼ CHðrHÞ28Npb2Nx Nb ¼ N0e H=1300 with N0 ¼ 5:24 10 3 s 1, and H is the resting water depth (in m). Equation (9) is similar in form to the drag coefficient of Jayne and St Laurent (2001) and likewise ignores the influence of critical turning latitudes (where x ¼ f , the Coriolis parameter) on the internal wave propagation characteristics. Yet, through scaling by x, the scheme is still frequency-dependent and thus applicable to only one specified constituent or, less strictly, to a particular tidal species. Considering our emphasis on the diurnal band, we fixed x to X, though that choice was found to improve the elevation accuracy of semidiurnal fringes (M2) as well. The IT drag formulation of Zaron and Egbert (2006) rather relies on scaling arguments than on a solid theoretical description of the topographically induced energy flux and thus contains a free parameter (C) to optimize the performance of the scheme. For practical reasons, we set C ¼ 50 (Green and Nycander 2013) and applied a secondary independent multiplier c at the order of Oð1Þ. With S1ð pÞ taken from Ray and Egbert (2004) , we could choose c in such a way that our model emulates the OAM values of this reference study. Alternatively, given the resemblance of S1 to the global character of K1 and O1, the RMS misfit of forward-modeled diurnal gravitational constituents to altimetry-constrained solutions can be optimized. Both criteria do not lead to fully rigorous tuning experiments, as ocean dynamics vary from one tide to the other and allowance must be made for subtle differences of our time-stepping model with respect to Ray and Egbert (2004) . In two separate suites of 40-day simulations with forcing specified for either fM2; O1; S1g or the purely gravitational combination of fM2; O1; K1g, c was varied in steps of 0.5 within a range of 0.5–4. Tuning by RMS differences of K1 and O1 to the observed tide favored c ¼ 1:5, whereas the best match with Ray and Egbert (2004) in terms of OAM was achieved by c values in the vicinity of 3, although the eventual prograde annual nutation results appeared to be only weakly dependent on the exact value of c (within about 10 las). ð9Þ ð10Þ ð11Þ 1.23 (163 ) 2.88 (304 ) 2.23 (218 ) a Air pressure tides S1ðpÞ for MERRA, ERA, and EC-OP are averages over 2004–2013, while the CFSR run draws on a 2004–2010 average. For brevity, the additional MERRA solution computed for the reduced 2004– 2010 window is not tabulated but given in terms of nutation in Table 10. Amplitudes are in units of 1023 kg m2 s 1 and cotidal phases are given relative to Greenwich noon b S1ðpÞ from Ray and Egbert (2004) was deployed for the control run to validate our hydrodynamic model configuration including the tidal conversion scheme As a trade-off, we adopted c ¼ 2 as a ‘‘best estimate’’ for all of our S1 runs below. RMS discrepancies to FES2012 produced by this setting are 5.6 cm (M2), 1.7 cm (O1), and 2.3 cm (K1), implying more than 93 % of sea surface height variance explained for each tide; cf. similar statistics obtained by Arbic et al. (2004) with their 1=2 barotropic model. S1 amplitude and phase charts in our updated control runs with IT drag included are nearly indistinguishable from Fig. 3 of Ray and Egbert (2004) , with previously noted regional anomalies (South China Sea, North Atlantic) eliminated (not shown). Table 8 underlines this sound agreement on the level of OAM values; cf. also the marked improvement with respect to our original, drag-free solution in Table 7. For both mass and motion components, phases differ by less than 15 throughout and deviations in amplitude are within 0:2 1023 kg m2 s 1. We calculated the corresponding contributions to the prograde annual nutation by aid of a standard protocol noted below, obtaining 21:9 þ i46:4 las as a credible reproduction of the S1 excitation value 20:7 þ i54:8 las implied by the OAM terms of Ray and Egbert (2004) ; see Table 10. A moderate underestimation of the op component in DEBOT ( 8 las) likely relates to differences in bathymetry or the treatment of ice shelves. Note also that the IT drag has no correspondence in the shallow water dynamics formulated by Ray and Egbert (2004) , although a pertinent parameterization of tidal conversion, incorporated to the very same numerical model by Egbert et al. (2004), might have gone unmentioned. Whether our prescription of topographically generated drag at the S1 frequency is physically justified or not is a potentially interesting issue but not of immediate importance for the topic in hand. One of the vexing problems related to this question is that our forward model operates in the time domain, while internal waves are preferably studied in the frequency domain. Here, we have adopted a diagnostic approach, inferring the need for additional mid-ocean dissipation by comparing our initial results for S1 and other diurnal tides to established reference charts. On a side note, also the discrepancies to coastal tide gauge estimates of S1 (next section) are markedly lower when internal wave drag is parameterized. 4.4 Hydrodynamic Solutions and Validation with Tide Gauge Data Tidal elevation charts obtained from numerical modeling are shown in Figs. 6, 7 and 8 for MERRA, CFSR, and EC-OP, while ERA has been left aside as the solution with probably the least accurate forcing data; cf. Sect. 3.3. Similarities with published S1 charts (e.g., Dobslaw and Thomas 2005; Ponte and Vinogradov 2007) are readily apparent and particularly striking for our EC-OP model as compared to the S1 tide of Ray and Egbert (2004) , who also employed ECMWF operational analysis data. Measured against Fig. 4 of these authors, DEBOT appears to underestimate the tide in Baffin Bay and the Sea of Okhotsk, probably due to differences in seafloor topography or the specification of dissipative processes. Such small-scale deficiencies in high latitudes are, however, of little relevance for the global OAM integrals. All computed S1 realizations agree in terms of the global character of the tide, but basinwide features can vary substantially in response to different pressure forcing data. Specifically, the secondary peaks of S1ð pÞ around 60 S in the CFSR climatology (Fig. 2) induce sea level signals in the Southern Ocean that exceed the corresponding tidal 20 15 10 5 0 360 300 240 180 120 60 0 60°N 30°N 0° 30°S 60°S 60°N 30°N 0° 30°S 60°S 60°E 120°E 180° 120°W 60°W 0° 60°E 120°E 180° 120°W 60°W 0° Fig. 6 Amplitudes (top, in mm) and Greenwich phase lags (bottom, in deg) for the sea level signal due to forcing by the S1 atmospheric pressure tide from MERRA (2004–2013 average). Cotidal phases are relative to Greenwich noon variability from MERRA and EC-OP by about 1–5 mm. Amplitudes in the North Atlantic are also comparatively high in the CFSR solution, whereas EC-OP displays the largest S1 tide in the Tropical Pacific, consistent with the overestimation of equatorial pressure gradients as exposed by Fig. 5c. Empirical knowledge of S1 to validate our forward simulations comes from globally distributed coastal tide gauges. Such a point-wise verification may imply little with regard to Earth rotation, yet it is instructive to determine whether individual models systematically perform better than others. Harmonic estimates of S1 at some 200 places are available in the online datasets of Ponchaut et al. (2001) 2, who tidally analyzed multi-year time series of hourly sea level records assembled both by BODC (British Oceanographic Data Centre) and UHSLC (University of Hawaii Sea Level Center). We extracted a subset of 51 estimates from Ponchaut’s compilation, excluding sites where the tide is effectively zero (Hawaii, Japan, Maldives, Central Atlantic) or all model predictions are equivocally different from the observations, e.g., due to unresolved coastal geometries (Prudhoe Bay, Bluff Harbour). Moreover, in order to avoid biases toward densely sampled regions, only a 60°E 120°E 180° 120°W 60°W 0° Fig. 7 Same as Fig. 6 but for S1ðpÞ from CFSR (2004–2010 average) 2 http://www.bodc.ac.uk/projects/international/woce/tidal_constants/. Accessed 9 October 2015. 60°N 30°N 0° 30°S 60°S 60°N 30°N 0° 30°S 60°S 60°E 120°E 180° 120°W 60°W 0° 20 15 10 5 0 360 300 240 180 120 60 0 few locations were retained in the equatorial Pacific and further thinning was applied to higher-amplitude S1 estimates in close proximity to each other (Arabian Sea, Gulf of Alaska). Our final 56-station set also includes five tide gauges from Ray and Egbert (2004) (Karachi, Benoa, Broome, Bermuda, Gibraltar) and is presented in Fig. 9. Valuable additions to this network, e.g., in the seas of Southeast Asia or along the coast of Brazil were pinpointed in the PSMSL (Permanent Service for Mean Sea Level) holdings, yet we refrained from a thorough tidal analysis of these hourly data in the frame of the present work. The collected harmonics are aggregate measures of both the radiational S1 ocean tide and the much smaller gravitationally driven component, S1g. The latter must be removed from the in situ data to rigorously compare with our numerical solutions of S1 that are 60°E 120°E 180° 120°W 60°W 0° Fig. 8 Same as Fig. 6 but for S1ðpÞ from EC-OP (2004–2013 average) 60°N 30°N 0° 30°S 60°S 60°N 30°N 0° 30°S 60°S 60°E 120°E 180° 120°W 60°W 0° 20 15 10 5 0 360 300 240 180 120 60 0 n H io .t s n S e c s A t s e llo i W n ye zan K a M d n a l s I s o ra a a l N a T a r lt a B d n a l s I e k a W g n u K n e h C i h c a r a K h a l a l a S a c i r A a v u S e t e e p a P a u c s a P a d a l s I z e d n a n r e F n a u J e c n a r e p s E l u a P . t S h t e b a z li n E a t b r r o u P D a v i H u k u N n y h r n e P n i w r a D o r B s a m t s i r h s C co o C a w a r a n o t n a K T a r . a i am no llie g H v in s p n a w K o T l a k a l a a o g n e M n u e m it B o B e l a M a i r o t c i V t r o P a s a b m o M ° 0 m i s i S k i v a j k y e R w o n r o t S 0 6 W 0 2 1 W 0 8 1 0 2 1 a d a g l e D a t n o P r a lt a r b i G e p a C ity D C n tr t a e n S p e u c s tta eR re u c C k n a ir Y P m m S t r o P o t r e u P O a z − n o ra C d e E a p e s s e E D P R S F C A R R E M e g u a G e d i n e l e u g r e w K a n o s M ° ° ° ° N N MERRA CFSR ERA EC-OP Amplitude \8 mma RMSb a Of 56 stations, 31 feature in situ amplitudes less than 8 mm b Numbers in parentheses are amplitude-normalized RMS differences solely forced by atmospheric pressure. To that end, we evaluated the Sg1 chart given in Appendix ‘‘The Gravitational S1 Ocean Tide’’ at the locations of our 56 gauges and changed the phase reference from the tide-generating potential to the simple radiational S1 argument; see Ray and Egbert (2004) for details. Tidal components after subtraction of the gravitational signal (usually 1–3 mm) are displayed as phasors in Fig. 9 and cover an amplitude range from 56 mm at Darwin (10-year mean estimate) down to 1.3 mm at St. Helena (4-year mean). Confidence intervals for all small-magnitude S1 determinations from only a few years of data appear to be sufficiently tight in the analysis of Ponchaut et al. (2001) to warrant the inclusion of these stations in our network. The median time series length over all 56 tide gauges is 8 years. Simulated S1 signals at each gauge location were taken from the nearest pelagic point in our 1=3 model and are illustrated for MERRA, CFSR, and EC-OP in Fig. 9. In general, all hydrodynamic solutions agree reasonably well with the observations, even though disparities on the order of a few mm must be accepted at most sites. The unusually large tide at Darwin (56.0 mm after reduction of the gravitational signal) has been addressed by Ray and Egbert (2004) and appears to be a very local modulation of S1 in a shallow (5-m) inlet that is approximated by a coarse gridpoint of 10 m depth in our bathymetry. MERRA and CFSR amplitudes at Darwin are 35 mm and 37 mm and thus somewhat closer to the observation than the model estimate of Ray and Egbert (2004) . In a broader context, the collection of tide gauges across the Atlantic testify to the shortcomings of the CFSR solution, evident, e.g., from the amplitude excess at Puerto Deseado, Port Stanley, and Esperanza. Moreover, most of the EC-OP estimates in the equatorial Pacific are too high, implying that this model must be treated with caution for studies of axial changes in Earth’s rotation. We have also attempted to express the varying accuracies of our hydrodynamic solutions by global statistical measures in Table 9. Median absolute differences as well as RMS misfits, given both as absolute and amplitude-normalized values, show little variations among the four models if all 56 tide gauge locations are considered. This result conforms with Fig. 9 inasmuch as the simulated tide at larger-amplitude sites tends to differ from the observation in the same way for all models; see, e.g., Karachi, Port Victoria, Yakutat, or all Australian stations. Somewhat more instructive statistics are obtained if the network is limited to stations below a certain amplitude threshold. Table 9 specifies results for an 8-mm threshold that preserves 31 tide gauges, most of them being located at mid-ocean islands. In this variant, MERRA outperforms all other models both in terms of MAD and RMS values, of which the reduction from 1.9 mm to 1.6 mm is significant at the 0.15 level. Results for the second MERRA run associated with the 2004–2010 pressure tide average (not shown) are slightly inferior, comparable with the MAD statistics of the two ECMWF models. Note also that the ERA tide, for which we have made reservations with regard to the forcing data (Sect. 3.3), is among the best models in Table 9 owing to a particularly good match with tide gauge estimates in the North Pacific. 4.5 Contribution of the Oceanic S1 Tide to Nutation Global OAM integrals derived from our numerical modeling efforts are compiled in Table 8 and exhibit a considerable scatter in accordance with the large-scale inter-model differences noted in the previous section. There is, however, a broad consensus that the x and y mass terms, i.e., the two single most important components with regard to nutation, are in the order of 1.0 1023 kg m2 s 1 (160 phase lag) and 3.0 1023 kg m2 s 1 (0 phase lag), respectively. The tabulated harmonics were translated to nutation values in essentially the same manner as AAM through multiplication of mass and motion terms with the proper transfer function coefficients T~p;wðrÞ; cf. Eqs. (1) and (2). As this scheme is initialized by a demodulation of angular momentum series in the time domain, we first discretized x and y OAM components over a 3-year window by a cosine function using the Doodson argument for the radiational S1 tide (see Appendix A of Ray and Egbert 2004) T þ 180 uH ð12Þ where T is Universal Time and uH are respective phase lag values as given in Table 8. Demodulated equatorial OAM series were then cleansed from non-seasonal signals (that is, the S1 contribution to prograde polar motion), fitted to the periodic forcing model (Eq. 4), and expressed as nutation harmonics through Eq. (1) with phases referred to the fundamental arguments of gravitational diurnal tides. Table 10 summarizes the various estimates of the oceanic S1 effect in nutation. Formal errors have been omitted as they are effectively zero given our usage of perfect sinusoids Table 10 Periodic oceanic contributions to the prograde annual nutation (las) as inferred from the modelspecific OAM values of Table 8a Total ip op Ray/Egbert Control MERRA MERRAb CFSRb ERA EC-OP Mass ip a Results are split up into the contributions from tidal heights (mass term) and currents (motion term). Inand out-of-phase components are referred to the fundamental arguments of nutation (Table 4) and the sign convention is that of Koot and de Viron (2011) b Forced by the respective pressure tide average from 2004–2010 for the angular momentum time series. Overall, the nutation results from all model runs are reasonably consistent, ranging from 0 to 20 las in the ip terms and roughly 40 to 60 las for the op component, with 90 % of the signal coming from the mass component. MERRA (2004–2013 average) and EC-OP produce a particularly close match within 6 las, and the moderate increase in magnitude for the reduced MERRA time span (2004–2010) is in fact expected on grounds of the time-variable amplitudes of S1ð pÞ (Fig. 5). For CFSR, the large-scale enhancement of tidal heights in the Indian Ocean and the North Atlantic (Fig. 7) combine to yield an excessive op estimate of 84.7 las. Nonetheless, the spread of nutation values is significantly smaller than that of previous inter-model comparisons, conducted, e.g., by Brzezin´ ski (2008 ) based on IB-corrected OAM values from much coarser ([1 ) barotropic and baroclinic models. We have also mapped the fine-resolution S1 tide of FES2012 to the prograde annual band, finding a nutation estimate of 2:1 þ i49:1 las that roughly matches our ERA harmonic. This agreement likely relates to similarities in the barometric forcing data, as the hydrodynamic core of FES2012 includes pressure loading from the 3-h ECMWF delayed cutoff stream (Carre`re et al. 2012) . Prograde annual nutation estimates for FES2012 as well as the model of Ray and Egbert (2004) are also tabulated in Schindelegger et al. (2015) , albeit with an internal conversion error at the order of 10 las which has been corrected in the frame of the present study. 5 Comparison with Geodetic Observations At an amplitude of 25.6 mas, the prograde annual nutation is among the principal signal components in Earth’s celestial motion and driven almost exclusively by the action of the solar gravitational torque on the equatorial bulge. In the MHB theory, the term is modulated to a minor degree by anelasticity ( 10 i4 las), electromagnetic torques ( 14 þ i6 las for both core mantle and inner core boundaries), geodesic nutation ( 30 þ i0 las), and the angular momentum exchange of the solid Earth with the gravitational ocean tide, S1g ( 21 þ i22 las); see also Table 2 of Brzezin´ ski et al. (2004). With these contributions accounted for, theory and observation of the prograde annual nutation produce a mismatch of 10:4 þ i108:2 las that has been attributed by MHB to the thermal atmospheric S1 tide and, implicitly, to the radiational S1 tide in the ocean. Realizations of the very same residual have been also derived by Koot et al. (2010) in the frame of a time domain Bayesian inversion of nutation observations including nonlinearities and additional terms in the functional model. Koot et al. (2010) used 10 years of additional VLBI data compared to MHB but employed identical corrections for geodesic nutation and the gravitational ocean tide. It is thus not surprising that the empirical S1 estimate of these authors is numerically very similar to the MHB residual; from a joint inversion of three nutation series from different analysis centers Koot et al. (2010) deduced a harmonic of 0 þ i107 las. Corresponding SD in both ip and op components are 4 las but probably underestimated and arguably better represented by the single-solution error of 7 las; cf. Table 1 of Herring et al. (2002). Residual VLBI-based nutations obtained after reduction of known effects do not necessarily provide a clean account of the rotational signal associated with the global S1 tide. Both unconsidered Sun-synchronous effects as well as inaccuracies in the incorporated relativistic or geophysical corrections at the prograde annual frequency might perturb empirical S1 estimates. However, theoretical values of geodesic nutation are known to great precision (Fukushima 1991), whereas anelastic and electromagnetic coupling contributions to the S1 band are too small (\20 las) to leave room for significant changes even if the MHB treatment of these effects is revised. The contribution from the gravitational ocean tide is somewhat larger (see above) and in fact subject to uncertainties owing to the manner in which it has been included in the nutation formalism. In detail, MHB inferred a harmonic of 21 þ i22 las from OAM estimates of K1, P1, O1, and Q1 (Chao et al. 1996) via scaling relationships that were optimized for the diurnal band on a broad scale instead of particular tidal lines. We therefore recomputed the effect based on S1g OAM integrals deduced in Appendix ‘‘The Gravitational S1 Ocean Tide’’ (Table 12), applying essentially the same time domain discretization as in Sect. 4.5 but with phases referred to the present-day argument of the gravitational S1 tide, that is, T þ 295:66 þ 90 (Ray and Egbert 2004) . Multiplication of adjusted mass and motion term coefficients (Eq. 4) with the respective transfer ratios (Eq. 1) yielded a harmonic of 15:2 þ i16:8 las. This value is about 5 las smaller than the intrinsic MHB estimate in both ip and op components, and a similar decrement is assumed for the analysis of Koot et al. (2010), who also utilized OAM data of Chao et al. (1996) . The corresponding correction was imposed on the prograde annual nutation residuals of both studies, resulting in the empirical S1 terms given in Table 11. Additional regard must be paid to the distortion of observed nutations through Sunsynchronous thermal deformations of some or all VLBI telescopes (Herring et al. 1991). This effect is now rigorously accounted for in VLBI analyses by means of a conventional procedure using on-site values of temperature, but the matter of discussion is whether a proper deformation correction was employed in the computation of nutation series that underlie the studies of MHB as well as Koot et al. (2010). Here, we draw on different evidences to argue that the effect was sufficiently well modeled, e.g., by early reduction schemes similar to Sovers et al. (1998) (Sect. G, ibid.). MERRA MERRAb CFSRb ERA EC-OP Brzezin´ski et al. (2004) Brzezin´ski (2011) VLBI (MHB) VLBI Koot and de Viron (2011) 1r error In-phase 8:0 6:1 32:3 38:7 9:4 113.1 60:6 16:2 5:8 7 a MERRA, CFSR, ERA, and EC-OP results are superpositions of the harmonics from Tables 5 and 10. For comparison, earlier estimates from Brzezin´ski et al. (2004) and Brzezin´ski (2011) are also shown and have been multiplied by 1 to account for differences in the definition of nutation amplitudes. VLBI values, with formal errors taken from Herring et al. (2002), have been cleared of the gravitational S1 tide influence by using the results of Appendix ‘‘The Gravitational S1 Ocean Tide’’; see the text for further details b Forced by the respective pressure tide average from 2004–2010 Single-session runs with our in-house VLBI software (Bo¨ hm et al. 2012) showed that a spurious prograde annual nutation variability of about 20 las in both ip and op components is incurred by analyses that explicitly omit corrections for solar heating of VLBI antennas. Hence, a persistent bias of 30 las should be evident in the comparison of nutation series from present-day VLBI analyses with the MHB model, assuming that for the latter diurnal deformation signals were neglected. This comparison is actually realized in the form of the IERS CPO data given w.r.t. the MHB series, of which a windowed Fourier analysis in the prograde annual band has already been presented in Fig. 1. The absence of any systematic distortion at the order of 30 las is readily apparent, in particular in the post-2000 period that features little uncertainty in the CPO estimates. Moreover, Table 3 of Koot et al. (2010) itself implies a proper modeling of solar heating in previous VLBI solutions. Among the three nutation series inverted by these authors, a thermal deformation correction is unambiguously identified for the IAA (Institute of Applied Astronomy, Moscow) data; see the corresponding documentation available at ftp://ivsopar.obspm.fr/vlbi/ ivsproducts/eops/ (accessed 14 October 2015). The treatment of the heating effect is undisclosed in the description of the other two series but, reassuringly, the associated prograde annual nutation residuals are not systematically offset from the IAA solution. Note also that the MHB estimate blends in well with Koot et al. (2010)’s results for various analysis centers. Artificial nutation signals related to antenna structure changes can be therefore deemed insignificant and the collected VLBI-based empirical S1 terms should be accurate enough to serve as reference values for our geophysical model estimates. This excitation balance is elaborated in Table 11 as well as in Fig. 10 and represents the core result of our study. Both MERRA (either solution) and EC-OP estimates agree with geodetic observations of the prograde annual nutation at the 10 las level, well below the threefold SD of the VLBI solutions. A discrepancy of only 3 las is found between MERRA (2004–2010) and the joint inversion residual of Koot et al. (2010), even though such a close fit might be fortuitous considering the time variability of S1 excitation terms 3σ VLBI (Koot) 3σ VLBI (MHB) MERRA CFSR ERA EC−OP 0 20 40 60 80 100 120 out−of−phase (µas) 140 160 180 and the uncertainties of the involved numerical models. By and large, the consistency of MERRA with VLBI data is in keeping with its good performance in comparison with atmospheric and oceanic ground truth data. Further correlations between model-specific in situ statistics and the nutation results in Fig. 10 are less obvious, but a closer inspection of tide gauge estimates in Durban, Port Elizabeth, Esperance, and Mawson (Fig. 9) revealed that ERA systematically underestimates the ocean tide in the South Indian Ocean, by about 2 mm compared to the observed sea level signal and to other simulations (MERRA, EC-OP). We implemented a split-up of the retrograde OAM mass term into contributions from different basins, showing that the broadscale S1 features in the South Indian Ocean have indeed a significant bearing on the op component of the ocean-driven prograde annual nutation. The 30-las deviation of the ERA model from, e.g., MHB’s nutation estimate is thus attributed in large part to a regional signal loss in terms of tidal elevation in the southern hemisphere. Table 11 also places our final results in the context of previous geophysical modeling efforts by Brzezin´ ski et al. (2004) and Brzezin´ ski (2011). In these studies, a balance with VLBI observations has been mostly impeded by deficiencies in the ip component of the modeled S1 excitation. To some extent, the choice of meteorological data (NCEP R1, ERA) is critical in either investigation, and further errors relate to insufficiencies of the utilized ocean models regarding the simulation of the S1 tide. Specifically, the model deployed by Brzezin´ ski (2011) has been optimized for a range of timescales and full baroclinic variability, for which coarse horizontal resolutions (1.875 ) greatly reduce computational costs. We tested the impact of a 1 discretization in DEBOT, obtaining somewhat anomalous S1 charts and increasingly negative ip components of the oceandriven prograde annual nutation ( 15 las decrement). Brzezin´ ski et al. (2004) analyzed the output of a barotropic model with a comparably coarse domain representation (1.125 spacing) but also with SAL dynamics neglected. This omission alters nutation amplitudes by roughly 30 las, and perturbations of similar size occur if dissipative processes are imperfectly accounted for. Such shortcomings have been redressed in the present work, by drawing on modern insights into the forward modeling of global ocean tides. 6 Concluding Discussion S1 tidal excitations of nutation in the order of 3.5 mm at the Earth’s surface ( 120 las) have constituted an anomaly to non-rigid nutation theories for decades. We have put forth an explanation of geodetic observations of the effect based on reanalysis data from MERRA and operational ECMWF analysis fields, complemented by numerical hydrodynamic solutions for the radiational S1 tide in the ocean. Atmospheric contributions averaged over 2004–2013 are 21:9 þ i45:8 las (MERRA) and 22:4 þ i67:5 las (EC-OP) and combine well with the respective oceanic estimates (13:9 þ i60:2 las, 13:0 þ i54:3 las) to match the VLBI-observed S1 terms within 10 las. No attempt was made to rigorously quantify the uncertainty of these geophysical model estimates, but we suppose that the errors are comparable to the threefold VLBI SD in the prograde annual band (21 las). In particular, the atmospheric mass term is among the least robust components of the global S1 excitation given its dependence on the small second-order tesseral surface pressure wave. Table 5 documents an inter-model spread of 25 las (excepting CFSR) for the pressure-driven nutation, comprising also uncertainties due to inter-annual S1 variations that have not been completely removed by the chosen 10-year averaging window. In contrast, our forward simulations of the radiational ocean tide should be fairly reliable on condition that the barometric forcing data themselves are accurate. Only weak (\10 las) and possibly counterbalancing influences of bathymetry and drag parameterization have been noted in Sect. 4. We also conducted DEBOT runs on C-grids finer than 1=3 , obtaining nutation harmonics of only a few las deviation with respect to the estimates given in Table 10. Differences in the diurnal cycle of modern atmospheric assimilation systems have played one of the recurring themes throughout this paper and are not necessarily smaller than the high-frequency disparities among earlier generation reanalyses; recall, e.g., our assessment of the CFSR pressure data. A more coherent representation of air tides is evidently tied to a near-global observing system with continuous sub-daily sampling (Schindelegger and Dobslaw 2016) but also depends on other aspects of the (re)analysis framework. Poli et al. (2013) emphasized the importance of at least hourly radiation time steps—a condition that is met neither by CFSR nor by JRA-55 (Kobayashi et al. 2015), which was also examined in a preliminary stage of our study but led to deficient AAM/ OAM phasors. Moreover, the formulation of the assimilation technique can have implications for tides, considering in particular that the variational analysis (3DVar or 4DVar; see Table 1) is usually performed in sequential 12-h windows without accounting for continuity of state variables at the transition epochs. The resulting perturbations occur at integer fractions of a solar day and potentially fold to an artificial S1/S2 variability. Such spurious signals are, however, minimized in the special case of MERRA through its Incremental Analysis Update method (Rienecker et al. 2008) , which might ultimately figure into the good performance of MERRA throughout our study. Finally, the accuracy of S1 in global analysis models is closely linked to the fidelity with which moist convection and latent heat flux can be simulated. Deficiencies in these quantities relate to imperfect physical paramaterizations or uncorrected biases in observations (Meynadier et al. 2010) and are, e.g., well documented for ERA (Dee et al. 2015). Displaying little long-term variability both in the celestial pole offsets (Fig. 1) and in the atmospheric S1 excitation, the 2004–2013 period has provided the ideal setting to study the mean harmonic atmosphere–ocean contribution to the prograde annual nutation. A reliable estimation of the temporal evolution of nutation amplitudes is still challenging, though. These signals differ substantially among the probed models and are masked by noise interferences as well as spurious variabilities when the frozen assimilation routines of reanalyses are confronted with new types and volumes of observations. Judging from Figs. 1, 3, and similar analyses in Bizouard et al. (1998) , an upper bound of 30 las appears to be a plausible estimate for the irregular departures from a simple sinusoidal S1 term in nutation. These vacillations dictate the likely accuracy of upcoming nutation models but also serve as an incentive for future foundational research, relating climate signals in geodetic observations with the time-variable excitation quantities from geophysical fluid models. Acknowledgments Open access funding provided by Austrian Science Fund (FWF). We are indebted to R. Ray for supplying various datasets and patiently answering questions about tides and ocean modeling. Comments from A. Brzezin´ski, H. Dobslaw, and M. Green are acknowledged, and we also thank C. Mayerhofer for coding some of the basic routines used in this work. The analyzed meteorological data were provided by the ECMWF, NASA’s GMAO, and the Research Data Archive (RDA) of NCAR. FES2012 was produced by Noveltis, Legos and CLS Space Oceanography Division and distributed by Aviso, with support from Cnes (http://www.aviso.altimetry.fr/). We greatly appreciate the Austrian Science Fund for financial support within project I1479-N29. D. Salstein is sponsored in part by Grant ATM-0913780 from the US National Science Foundation (NSF), and D. Einsˇpigel thanks the Charles University in Prague for supporting him under Grant SVV-2015-260218. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 2 1 360 300 240 180 120 60 Appendix 1: The Gravitational S1 Ocean Tide Following Sect. 2c of Ray and Egbert (2004) , a modern-day chart of S1g in the global ocean is readily computed from observations of the gravitational P1 and K1 tides, which are separated from the S1 band by only 1 cpy. Tidal heights of K1 were extracted from the FES2012 atlas on a 1=16 mesh, moderately downsampled, and scaled to local Sg1 amplitudes using a factor of 1.98/368.74, that is, the ratio of gravitational potentials at S1 and K1. Greenwich phase lags were calculated as averages from both the K1 and P1 charts and are illustrated together with the amplitudes of Sg1 in Fig. 11. Associated barotropic 60°E 120°E 180° 120°W 60°W 0° Fig. 11 Cotidal elevation charts for the gravitational S1g tide in the global ocean as computed from FES2012 solutions of K1 and P1. Amplitudes (in mm) are shown in the upper panel, Greenwich phase lags ( ) in the lower panel 60°N 30°N 0° 30°S 60°S 60°N 30°N 0° 30°S 60°S 60°E 120°E 180° 120°W 60°W 0° currents follow from the velocity (u) grids of K1 and P1 in the FES2012 model, based on the same interpolation procedures as for the local elevation. To serve Sect. 5, we have mapped heights and currents of S1g to global OAM values as documented in Table 12. Note that these harmonics might be also derived from a direct application of admittance relationships to the angular momentum values of the K1 and P1 tides. Dee D, Fasullo J, Shea D, Walsh J (2015) The climate data guide: atmospheric reanalysis: overview & comparison tables. https://climatedataguide.ucar.edu/climate-data/atmospheric-reanalysis-overviewcomparison-tables. Accessed 13 July 2015 Dehant V, Defraigne P (1997) New transfer functions for nutations of a nonrigid Earth. J Geophys Res 102(B12):27659–27687. doi:10.1029/97JB02347 Dehant V, Arias F, Bizouard C et al (1999) Considerations concerning the non-rigid earth nutation theory. Celest Mech Dyn Astron 72:245–310 Dehant V, Feissel-Vernier M, de Viron O, Ma C, Yseboodt M, Bizouard C (2003) Remaining error sources in the nutation at the submilliarc second level. J Geophys Res 108(B5):2275. doi:10.1029/ 2002JB001763 Dehant V, Mathews PM (2009) Earth rotation variations. In: Herring T (ed) Treatise on Geophysics, vol 3. Geodesy, pp 295–350 de Viron O, Boy JP, Goosse H (2004) Geodetic effects of the ocean response to atmospheric forcing in an ocean general circulation model. J Geophys Res 109 (B03411). doi:10.1029/2003JB002837 Dobslaw H, Thomas M (2005) Atmospheric induced oceanic tides from ECMWF forecasts. Geophys Res Lett 32 (L10615). doi:10.1029/2005GL022990 Egbert GD, Ray RD, Bills BG (2004) Numerical modeling of the global semidiurnal tide in the present day and in the last glacial maximum. J Geophys Res 109:C03003. doi:10.1029/2003JC001973 Einsˇpigel D, Martinec Z (2015) A new derivation of the shallow water equations in geographical coordinates and their application to the global barotropic ocean model (the DEBOT model). Ocean Model 92:84–100. doi:10.1016/j.ocemod.2015.05.006 Fedorov EP, Smith ML, Bender PL (1980) Nutation and the Earth’s rotation. In: Proceedings of the of IAU symposium no. 78, D Reidel. Dordrecht, Netherlands Fukushima T (1991) Geodesic nutation. Astron Astrophys 244:L11–L12 Green JAM, Nycander J (2013) A comparison of tidal conversion parameterizations for tidal models. J Phys Oceanogr 43:104–119. doi:10.1175/JPO-D-12023.1 Hendershott M (1972) The effects of solid earth deformation on global ocean tides. Geophys J R Astron Soc 29:389–402 Herring TA, Buffett BA, Mathews PM, Shapiro II (1991) Forced nutations of the Earth: influence of inner core dynamics: 3. Very long interferometry data analysis. J Geophys Res 96(B5):8259–8273. doi:10. 1029/90JB02177 Herring TA, Mathews PM, Buffett BA (2002) Modeling of nutation-precession: very long baseline interferometry results. J Geophys Res 107(B4):2069. doi:10.1029/2001JB000165 Jayne SR, St Laurent LC (2001) Parameterizing tidal dissipation over rough topography. Geophys Res Lett 28:811–814. doi:10.1029/2000GL012044 Jeffreys H, Vicente RO (1957) The theory of nutation and the variation of latitude. Mon Not R Astron Soc 117:142–161 Kinoshita H (1977) Theory of the rotation of the rigid Earth. Celestial Mech 15:277–326 Kobayashi S, Ota J, Harada J et al (2015) The JRA-55 reanalysis: general specifications and basic characteristics. J Meteorol Soc Jpn 93:5–48. doi:10.2151/jmsj.2015-001 Koot L, Dumberry M, Rivoldini A, de Viron O, Dehant V (2010) Constraints on the coupling at the coremantle and inner core boundaries inferred from nutation observations. Geophys J Int 182:1279–1294. doi:10.1111/j.1365-246X.2010.04711.x Koot L, de Viron O (2011) Atmospheric contributions to nutations and implications for the estimation of deep Earth’s properties from nutation observations. Geophys J Int 185:1255–1265 Li Y, Smith RB, Grubisˇic´ V (2009) Using surface pressure variations to categorize diurnal valley circulations: experiments in Owens Valley. Mon Weather Rev 137:1753–1769. doi:10.1175/2008MWR2495.1 Lieberman RS, Riggin DM, Ortland DA, Nesbitt SW, Vincent RA (2007) Variability of mesospheric diurnal tides and tropospheric diurnal heating during 1997–1998. J Geophys Res 112 (D20110). doi:10.1029/ 2007JD008578 Mathews PM, Herring TA, Buffett BA (2002) Modeling of nutation and precession: new nutation series for nonrigid Earth and insights into the Earth’s interior. J Geophys Res 107(B4):2068. doi:10.1029/ 2001JB000390 McPhaden MJ, Ando K, Bourle`s B et al. (2010) The Global Tropical Moored Buoy Array. In: Hall J, Harrison DE, Stammer D (eds) Proceedings of the Oceanobs’09: Sustained Ocean Observations and Information for Society Conference, vol 2. ESA Publ WPP-306 Meynadier R, Bock O, Gervois S, Guichard F, Redelsperger JL, Agusti-Panareda A, Beljaars A (2010) West African monsoon water cycle. 2: assessment of numerical weather prediction water budgets. J Geophys Res 115 (D19107). doi:10.1029/2010JD013919 Molodensky MS (1961) The theory of nutation and diurnal earth tides. Comm Obs Roy Belgique 188:25–56 Accad Y , Pekeris CL ( 1978 ) Solution of the tidal equations for the M2 and S2 tides in the world oceans from a knowledge of the tidal potential alone . Philos Trans R Soc London 290 : 235 - 266 Amante C , Eakins BW ( 2009 ) ETOPO1 1 arc-minute global relief model: procedures, data sources and analysis . NOAA Technical Memorandum NESDIS NGDC-24 , National Geophysical Data Center, NOAA. doi: 10 .7289/V5C8276M Arbic BK , Garner ST , Hallberg RW , Simmons HL ( 2004 ) The accuracy of surface elevations in forward global barotropic and baroclinic tide models . Deep Sea Res II 51 : 3069 - 3101 Bell TH ( 1975 ) Topographically generated internal waves in the open ocean . J Geophys Res Oceans 80 : 320 - 327 . doi: 10 .1029/JC080i003p00320 Bizouard C , Brzezin´ski A , Petrov S ( 1998 ) Diurnal atmospheric forcing and temporal variations of the nutation amplitudes . J Geod 72 : 561 - 577 Bizouard C , Lambert S ( 2002 ) Lunisolar torque on the atmosphere and Earth's rotation . Planet Space Sci 50 : 323 - 333 . doi: 10 .1016/S0032- 0633 ( 01 ) 00120 - 9 Bo¨hm J , Bo¨hm S , Nilsson T , Pany L , Plank L , Spicakova H , Teke K , Schuh H ( 2012 ) The new Vienna VLBI Software VieVS . In: Kenyon S, Pacino MC , Marti U (eds) Proceedings of the IAG scientific assembly 2009 , International Association of Geodesy Symposia Series, vol 136 , pp 1007 - 1011 . doi: 10 .1007/ 978-3- 642 -20338-1_ 126 Brzezin´ski A ( 1994 ) Polar motion excitation by variations of the effective angular momentum function, II: Extended model . Manuscr Geod 19 : 157 - 171 Brzezin´ski A , Ponte R , Ali A ( 2004 ) Nontidal oceanic excitation of nutation and diurnal/semidiurnal polar motion revisited . J Geophys Res 109 ( B11407 ). doi:10.1029/2004JB003054 Brzezin´ski A ( 2008 ) On the influence of diurnal atmospheric tides on Earth rotation . In: Capitaine N (ed) Proceedings of the Journe´es Syste`mes de Re´fe´rence Spatio-Temporels 2007 , pp 180 - 183 Brzezin´ski A ( 2011 ) Diurnal excitation of Earth rotation estimated from recent geophysical models . In: Capitaine N (ed) Proceedings of the Journe´es Syste`mes de Re´fe´rence Spatio-Temporels 2010 , pp 131 - 186 Carre`re L , Lyard F , Cancet M , Guiloot A , Roblou L ( 2012 ) FES 2012: a new global tidal model taking advantage of nearly 20 years of altimetry . Paper presented at the Symposium 20 Years of Progress in Radar Altimetry , ESA SP , Venice, Italy, 24 -29 Sept Chao BF , Ray RD , Gipson JM , Egbert GD , Ma C ( 1996 ) Diurnal/Semidiurnal polar motion excited by oceanic tidal angular momentum . J Geophys Res 101 ( B9 ): 20151 - 20163 . doi: 10 .1029/96JB01649 Chapman S , Lindzen RS ( 1970 ) Atmospheric tides . D Reidel , Dordrecht Compo GP , Whitaker JS , Sardeshmukh PD et al ( 2011 ) The Twentieth Century Reanalysis Project . Q J R Meteorol Soc 137 : 1 - 28 . doi: 10 .1002/qj.776 Dai A , Wang J ( 1999 ) Diurnal and semidiurnal tides in global surface pressure fields . J Atmos Sci 56 : 3874 - 3891 Dee DP , Uppala SM , Simmons AJ et al ( 2011 ) The ERA-Interim reanalysis: configuration and performance of the data assimilation system . Q J R Meteorol Soc 137 : 553 - 597 . doi: 10 .1002/qj.828 Padman L , Fricker HA , Coleman R , Howard S , Erofeeva L ( 2002 ) A new tide model for the Antarctic Ice shelves and seas . Ann Glaciol 34 : 247 - 254 Parke ME ( 1982 ) O1, P1, N2 models of the global ocean tide on an elastic earth plus surface potential and spherical harmonic decompositions for M2, S2, and K1 . Mar Geod 6 : 35 - 81 Petit G , Luzum B , ( 2010 ) IERS Conventions 2010 . IERS technical note no. 36. Verlag des Bundesamtes fu¨ r Kartographie und Geoda¨sie, Frankfurt am Main, Germany Poli P , Hersbach H , Tan T et al. ( 2013 ) The data assimilation system and initial performance evaluation of the ECMWF pilot reanalysis of the 20th-century assimilating surface observations only (ERA-20C) . ERA report series no. 14 , European Centre for Medium-Range Weather Forecasts . Reading, UK Ponchaut F , Lyard F , Le Provost C ( 2001 ) An analysis of the tidal signal in the WOCE sea level dataset . J Atmos Ocean Technol 18 : 77 - 91 Ponte RM , Vinogradov SV ( 2007 ) Effects of stratification on the large-scale ocean response to barometric pressure . J Phys Oceanogr 37 : 245 - 258 Ray RD ( 1998a ) Diurnal oscillations in atmospheric pressure at twenty-five small oceanic islands . Geophys Res Lett 25 : 3851 - 3854 . doi: 10 .1029/1998GL900039 Ray RD ( 1998b ) Ocean self-attraction and loading in numerical tidal models . Mar Geod 21 : 181 - 192 Ray RD , Ponte RM ( 2003 ) Barometric tides from ECMWF operational analyses . Ann Geophys 21 : 1897 - 1910 . doi: 10 .5194/angeo-21-1897-2003 Ray RD , Egbert GD ( 2004 ) The global S1 tide . J Phys Oceanogr 34 : 1922 - 1935 Rienecker MM , Suarez MJ , Todling R et al. ( 2008 ). The GEOS-5 data assimilation system- Documentation of versions 5.0.1, 5.1.0, and 5.2.0. NASA technical report series on global modeling and data assimilation , vol 27 . Greenbelt , Maryland, USA Rienecker MM , Suarez MJ , Gelaro R et al ( 2011 ) MERRA: NASA's Modern-Era Retrospective Analysis for Research and Applications . J Clim 24 : 3624 - 3648 . doi: 10 .1175/ JCLI-D- 11- 00015 . 1 Robertson FR , Bosilovich MG , Chen J , Miller TL ( 2011 ) The effect of satellite observing system changes on MERRA water and energy fluxes . J Clim 24 : 5197 - 5217 . doi: 10 .1175/2011JCLI4227.1 Saha S , Moorthi S , Pan HL et al ( 2010 ) The NCEP Climate Forecast System Reanalysis . Bull Am Meteorol Soc 91 : 1015 - 1057 Sasao T , Okubo S , Saito M ( 1980 ) A simple theory on the dynamical effects of a stratified fluid core upon nutational motion of the Earth . In: Fedorov EP , Smith ML , Bender PL (eds) Proceedings of the IAU symposium no. 78, Nutation and the Earth's Rotation , pp 165 - 183 Sasao T , Wahr JM ( 1981 ) An excitation mechanism for the free 'core nutation' . Geophys J R Astron Soc 61 : 729 - 746 Schindelegger M , Bo¨hm S , Bo¨hm J , Schuh H ( 2013 ) Atmospheric effects on Earth rotation . In: Bo¨hm J , Schuh H (eds) Atmospheric Effects in Space Geodesy . Springer, pp 181 - 231 Schindelegger M , Bo¨hm J , Salstein DA ( 2015 ) The global S1 tide and Earth's nutation . In: Malkin Z, Capitaine N (ed) Proceedings of the Journe´es Syste`mes de Re´fe´rence Spatio-Temporels 2014 , pp 145 - 150 Schindelegger M , Dobslaw H ( 2016 ) A global ground truth view of the lunar air pressure tide L2 . J Geophys Res Atmos 121 : 95 - 110 . doi: 10 .1002/2015JD024243 Smith A , Lott N , Vose R ( 2011 ) The Integrated Surface Database: recent developments and partnerships . B Am Meteorol Soc 92 : 704 - 708 . doi: 10 .1175/2011BAMS3015.1 Sovers OJ , Fanselow JL , Jacobs CS ( 1998 ) Astrometry and geodesy with radio interferometry: experiments, models, results . Rev Mod Phys 70 : 1393 - 1454 Stepanov VN , Hughes CW ( 2004 ) Parameterization of ocean self-attraction and loading in numerical models of the ocean circulation . J Geophys Res 109 ( C03037 ). doi:10.1029/2003JC002034 Taguchi E , Stammer D , Zahel W ( 2014 ) Inferring deep ocean tidal energy dissipation from the global highresolution data-assimilative HAMTIDE model . J Geophys Res Oceans 119 : 4573 - 4592 . doi: 10 .1002/ 2013JC009766 Vial F , Lott F , Teitelbaum H ( 1994 ) A possible signal of the El Nin˜o-Southern Oscillation in time series of the diurnal tide . Geophys Res Lett 21 : 1603 - 1606 . doi: 10 .1029/94GL01016 Wahr J ( 1981 ) The forced nutations of an elliptical, rotating, elastic and oceanless earth . Geophys J Int 64 : 705 - 727 Wilmes SB , Green JAM ( 2014 ) The evolution of tides and tidal dissipation over the past 21,000 years . J Geophys Res Oceans 119 : 4083 - 4100 . doi: 10 .1002/2013JC009605 Woolard EW ( 1953 ) Theory of the rotation of the earth around its center of mass . Astron J 58:2 Yseboodt M , de Viron O , Chin TM , Dehant V ( 2002 ) Atmospheric excitation of the Earth's nutation: comparison of different atmospheric models . J Geophys Res 107 ( B2 ): 2036 . doi: 10 .1029/2000JB000042 Zaron ED , Egbert GD ( 2006 ) Estimating open-ocean barotropic tidal dissipation: the Hawaiian Ridge . J Phys Oceanogr 36 : 1019 - 1035


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs10712-016-9365-3.pdf

Michael Schindelegger, David Einšpigel, David Salstein, Johannes Böhm. The Global S \(_1\) Tide in Earth’s Nutation, Surveys in Geophysics, 2016, 643-680, DOI: 10.1007/s10712-016-9365-3