# Constraints on the coupling at the core-mantle and inner core boundaries inferred from nutation observations

L. Koot, M. Dumberry, A. Rivoldini, O. de Viron, V. Dehant

<figure><img src="image_1.png" /></figure>

## To cite this version:

L. Koot, M. Dumberry, A. Rivoldini, O. de Viron, V. Dehant. Constraints on the coupling at the core-mantle and inner core boundaries inferred from nutation observations. Geophysical Journal International, 2010, 182, pp.1279-1294. 10.1111/j.1365-246X.2010.04711.x. insu-03605291

HAL Id: insu-03605291

https://insu.hal.science/insu-03605291

Submitted on 11 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Constraints on the coupling at the core–mantle and inner core boundaries inferred from nutation observations

L. Koot,<sup>1,2</sup> M. Dumberry,<sup>1</sup> A. Rivoldini,<sup>2</sup> O. de Viron<sup>3</sup> and V. Dehant<sup>2</sup>

Accepted 2010 June 21. Received 2010 May 31; in original form 2010 March 9

## SUMMARY

We present an inversion of nutation observations in terms of parameters characterizing the Earth's interior properties. We use a Bayesian inversion in the time-domain, allowing us to take fully into account non-linearities in the nutation model and to reduce the loss of information occurring in frequency-domain inversions. Among the parameters we retrieve are two complex parameters, $K_{\text{CMB}}$ and $K_{\text{ICB}}$, referred to as 'coupling constants', characterizing the mechanical coupling at the core–mantle boundary (CMB) and the inner core boundary (ICB), respectively. Based on a joint inversion of nutation observations provided by different analysis centres, we find $\text{Im}(K_{\text{CMB}}) = (-1.78 \pm 0.02) 10^{-5}$, $\text{Re}(K_{\text{ICB}}) = (1.01 \pm 0.02) 10^{-3}$ and $\text{Im}(K_{\text{ICB}}) = (-1.09 \pm 0.03) 10^{-3}$ (where the errors correspond to 99.7 per cent confidence intervals). While our value of $\text{Im}(K_{\text{CMB}})$ is similar to previous estimates, our new values of $\text{Re}(K_{\text{ICB}})$ and $\text{Im}(K_{\text{ICB}})$ are significantly different. This is mainly because of the different inversion strategy that we use and also because of the lengthier record of observation available. We show that, based on existing coupling models, neither viscous nor electromagnetic coupling alone can explain our new values of $\text{Re}(K_{\text{ICB}})$ and $\text{Im}(K_{\text{ICB}})$. A combination of these two mechanisms is required and necessitates a radial magnetic field at the ICB of total rms strength between 6 and 7 mT and a kinematic viscosity of the fluid core at the ICB should be between 10 and $30 \text{ m}^2 \text{ s}^{-1}$, depending on the exact partition between viscous and electromagnetic coupling.

**Key words:** Inverse theory; Earth rotation variations; Core, outer core and inner core.

## 1 INTRODUCTION

The gravitational (tidal) forces from the Moon, the Sun and other planets apply a torque on the elliptical Earth. In response to this torque, the orientation of the Earth's rotation and figure axes varies in time when observed from a frame fixed in space, giving rise to precession/nutation motion. The precession is the large secular part of the motion while nutations are defined as the small periodic variations with periods larger than two days.

The amplitudes of the nutations depend on the Earth's ellipticity, internal density distribution, and rheological properties. For a three-layer Earth comprised of a mantle, fluid outer core and solid inner core, each of these individual regions would respond differently to the tidal torque on account of their different densities and ellipticities. As the outer core is fluid, differential rotations between the three layers can occur. The Earth's nutations are then also dependent on the interaction between the fluid core and mantle at the core–mantle boundary (CMB) and that between the fluid core and inner core at the inner core boundary (ICB).

Moreover, the amplitudes of the forced nutations are amplified by the presence of two normal modes, the 'free core nutation' (FCN)

and the 'free inner core nutation' (FICN), whose frequencies are close to those of the tidal forcing. These modes are characterized by a differential rotation between the mantle, the fluid outer core and the solid inner core and mechanical couplings at boundaries affect their natural frequencies. Similarly, the damping of these modes depends on the energy dissipated through this coupling. The observed amplitude and phase of the forced nutations depend thus directly on the frequencies of the normal modes, which in turn depend on the mechanical coupling at the CMB and ICB.

Many of the physical properties at these boundaries on which the couplings depend are not well known. For instance, viscous coupling requires a knowledge of the viscosity of the fluid core close to the boundaries. Measurements are performed in laboratories on liquid iron alloys (e.g. Rutter et al. 2002a; Rutter et al. 2002b) but those experiments are carried on at much lower pressures (~6 GPa) and temperatures (~2000 K) than those typical of the outer core. Similarly, electromagnetic (EM) coupling depends on the electrical conductivity of both the fluid core and of the solid side of the boundary (mantle or inner core) as well as on the magnetic field at the boundary. We only have partial knowledge of the magnetic field at the CMB, the short-wavelength being masked by the crustal

<footnote n="1"><sup>1</sup>Department of Physics, University of Alberta, Edmonton, Canada. E-mail: koot@ualberta.ca</footnote>

<footnote n="2"><sup>2</sup>Royal Observatory of Belgium, Brussels, Belgium</footnote>

<footnote n="3"><sup>3</sup>Université Paris Diderot and Institut de Physique du Globe de Paris, associés au CNRS, Paris, France</footnote>
magnetic field (e.g. Stacey & Davis 2008). The magnetic field at the ICB cannot be observed at all due to the screening effect of the conducting core.

The tidal torque is known very precisely from celestial mechanics (Bretagnon et al. 1998; Roosbeek & Dehant 1998; Souchay et al. 1999). Without a precise knowledge of the key parameters that participate in the coupling at boundaries, it is then difficult to predict accurately the nutational response of the Earth. However, conversely, precise observations of Earth’s nutations offer an opportunity to obtain information on these internal quantities for which we have little or no other constraints.

The nutational response of the Earth to the tidal torque has been the subject of many studies, including the model of Mathews et al. (2002) on which our study is based. This model predicts the nutation motion of an Earth model at hydrostatic equilibrium comprised of an elastic solid inner core, an inviscid fluid outer core and an anelastic mantle surrounded by an ocean layer. The gravitational interaction between the internal regions as well as the fluid pressure acting on the elliptical boundaries are included in the model. Other possible coupling mechanisms are not explicitly modelled but are taken into account by means of two complex parameters, referred to as ‘coupling constants’, $K_{CMB}$ and $K_{ICB}$, characterizing the coupling at the CMB and ICB, respectively.

The key internal parameters of the nutation model, among them $K_{CMB}$ and $K_{ICB}$, can thus be determined by an inversion procedure, where one infers knowledge on these parameters from a comparison between the nutation observations and model. A first estimation of $K_{CMB}$ and $K_{ICB}$ was obtained by Mathews et al. (2002). They used a linearized least-squares inversion in the frequency domain, only including the most dominant forced nutations that could be reliably estimated from the observations (Herring et al. 2002). In Koot et al. (2008), we developed a different inversion procedure directly in the time domain, thus avoiding any loss of information from retaining only a few of the forced nutations. Additionally, we used Bayesian inversion method, allowing us to take better into account the non-linearity of the nutation model.

In this study, we further improve on the inversion procedure developed by Koot et al. (2008). We now include in the model non-periodic variations in the tidal forcing, the so-called ‘Poisson terms’, and the freely excited FCN mode. In addition, we use nutation time-series computed by different Very Long Baseline Interferometry (VLBI) analysis centres. This allows us to determine the accuracy with which we can recover the parameters that have an influence on nutations.

Our analysis is focused on the coupling constants $K_{CMB}$ and $K_{ICB}$. We find that our retrieved value for $K_{CMB}$ is in agreement with that of Mathews et al. (2002) while $K_{ICB}$ is markedly different. Since $K_{ICB}$ characterizes the coupling at the ICB, this difference implies that previous estimates of physical parameters entering coupling models at the ICB must be revised. The coupling constant $K_{ICB}$ has been previously interpreted in terms of EM coupling at the ICB (Buffett et al. 2002). Mathews & Guo (2005) and Deleplace & Cardin (2006) generalized the model of Buffett et al. (2002) in order to take viscous coupling into account. In this paper, we concentrate our attention on these two coupling mechanisms. We investigate the fluid viscosity and the amplitude and structure of the magnetic field at the ICB that are compatible with our new value of $K_{ICB}$.

Our paper is organized as follow. We start by presenting a brief description of the nutation model in Section 2. In Sections 3 and 4, we describe the nutation data sets that we use in our study and our inversion strategy, respectively. Results are presented in Sections 5 and 6, the latter also including the interpretation of our

coupling constants in terms of physical parameters. We conclude by a summary and a discussion of our results in Section 7.

# 2 NUTATION MODEL

## 2.1 Forced nutations

From the perspective of an observer rotating with the Earth, the action of the gravitational (tidal) potential ($\tilde{\phi}$) on the elliptical Earth leads to small deviations from a uniform angular rotation $\Omega_0$ about an axis $\mathbf{i}_3$ aligned with the geometrical figure. This deviation is termed ‘wobble’. The perturbed rotation vector is then $\Omega = \Omega_0(\mathbf{i}_3 + \mathbf{m})$, where $\mathbf{m}$ is a non-dimensional vector characterizing the wobble. The perturbation in the rotation vectors of the fluid core ($\Omega_f$), and solid inner core ($\Omega_s$) can be similarly defined as $\Omega_f = \Omega_0(\mathbf{i}_3 + \mathbf{m} + \mathbf{m}_f)$ and $\Omega_s = \Omega_0(\mathbf{i}_3 + \mathbf{m} + \mathbf{m}_s)$, where $\mathbf{m}_f$ and $\mathbf{m}_s$ are the wobbles of the fluid core and the solid inner core, respectively.

A prediction of the response of the Earth (in terms of these wobbles) from a given tidal forcing depends on the Earth’s interior properties and on the coupling between the mantle, fluid core and inner core. For our study, we use the model of Mathews et al. (2002), an updated version of the model developed in Mathews et al. (1991a). We only present here a brief description of the model and refer the interested reader to the original articles. The model consists of a system of three dynamic equations relating the equatorial angular momentum of the whole Earth, outer core and inner core, to the equatorial torques acting on them. A fourth kinematic equation describes the orientation of the inner core. As the deviations from the steady rotation state are expected to be very small, the dynamic equations are developed to the first order in the wobbles.

The angular momentum of each internal region is given by the product of its moment of inertia tensor and its instantaneous rotation vector. The former are expressed in terms of mean principal moments of inertia ($A, A_f$ and $A_s$ in the equatorial plane, $C, C_f$ and $C_s$ in the polar direction, respectively for the whole Earth, the fluid core, and the solid inner core) and of small corrections $c_{ij}$, $c_{ij}^f$ and $c_{ij}^s$ (respectively for the whole Earth, the outer core and the inner core) due to deformations. These are from three different sources: (1) tidal forces (characterized by the tidal potential $\tilde{\phi}$), (2) departure of the centrifugal force from that of the steady axial rotation (characterized by the wobbles $\mathbf{m}, \mathbf{m}_f$ and $\mathbf{m}_s$) and (3) the time-dependent surface loading from the ocean (characterized by a loading potential $\tilde{\phi}_L$). To the first order, the $c_{ij}$'s can be expressed as (Sasao et al. 1980; Sasao & Wahr 1981; Mathews et al. 1991a):

$$ \tilde{c}_3 = c_{31} + i c_{32} = A[\kappa(\tilde{m} - \tilde{\phi}) + \xi \tilde{m}_f + \zeta \tilde{m}_s + \chi \tilde{\phi}_L] \quad (1a) $$

$$ \tilde{c}_3^f = c_{31}^f + i c_{32}^f = A_f[\gamma(\tilde{m} - \tilde{\phi}) + \beta \tilde{m}_f + \delta \tilde{m}_s + \eta \tilde{\phi}_L] \quad (1b) $$

$$ \tilde{c}_3^s = c_{31}^s + i c_{32}^s = A_s[\theta(\tilde{m} - \tilde{\phi}) + \alpha \tilde{m}_f + \nu \tilde{m}_s + \lambda \tilde{\phi}_L], \quad (1c) $$

where the tilde stands for the complex combination of the two equatorial components, for example, $\tilde{m} = m_1 + im_2$. Eqs (1) introduce 12 so-called ‘compliances’ or ‘generalized Love numbers’ ($\kappa, \xi, \zeta, \chi, \gamma, \beta, \delta, \eta, \theta, \alpha, \nu$ and $\lambda$), characterizing the deformation of each region in response to a given forcing. For an elastic Earth, these compliances are real numbers. Here, we follow Mathews et al. (2002) and take mantle anelasticity into account, in which case the compliances are complex numbers. Since nutations are located in a very narrow frequency band, the compliances are assumed to be frequency independent.
In addition to the external tidal torque, misalignments between the rotation vector of each region lead to internal torques between them. Inertial and gravitational torques have been computed by Sasao et al. (1980) and Mathews et al. (1991a) to the first order in the perturbations induced by the tidal potential. The torques from other coupling mechanisms, such as electromagnetic, viscous and topographic, are not modelled explicitly in Mathews et al. (2002) but are introduced in a parametrized way. The amplitude of these torques depends on the differential rotation at region boundaries. To first order, the equatorial torque applied on the outer core by the mantle (Γ<sub>CMB</sub>) and that applied on the inner core by the outer core (Γ<sub>ICB</sub>) can be written as (Mathews et al. 2002)

$$
\tilde{\Gamma}_{\text{CMB}} = -i \Omega_0^2 A_f K_{\text{CMB}} \tilde{m}_f \quad (2a)
$$

$$
\tilde{\Gamma}_{\text{ICB}} = -i \Omega_0^2 A_s K_{\text{ICB}} (\tilde{m}_s - \tilde{m}_f), \quad (2b)
$$

where $\tilde{m}_f$ and $(\tilde{m}_s - \tilde{m}_f)$ characterize the differential rotation be-
tween the mantle and outer core, and between the outer and in-
ner cores, respectively. The non-dimensional parameters $K_{\text{CMB}}$ and
$K_{\text{ICB}}$ are constants characterizing the strength of the coupling. They
are complex parameters: their norm characterizes the strength of the
coupling and the imaginary part the amount of dissipation.

The first-order dynamic equations form a system of coupled lin-
ear differential equations in the variables ˆm, ˆm<sub>f</sub>, ˆm<sub>s</sub>, and ˆn<sub>s</sub>, this
latter variable representing the deviation of the figure axis of the
inner core from that of the mantle. Solutions for the dynamic vari-
ables depend on the forcing from the tidal potential ˆφ. Because the
motion of the celestial bodies is nearly periodic, to a very good
degree of approximation, the tidal potential can be written as a sum
of periodic terms

$$
\hat{\phi}(t) = \sum_{l=1}^{N} \hat{\phi}_0(\sigma_l) e^{i\sigma_l \Omega_0 t}, \quad (3)
$$

where $\sigma_l$ is the angular frequency in cycles per sidereal day (cpsd),
as seen from an Earth-fixed reference frame.

As the dynamic differential equations are linear, they can be
solved independently for each term $\hat{\phi}_0(\sigma_l)$ of the tidal potential
(3). The four dynamic variables can also be expressed as a sum of
periodic terms, such as

$$
\hat{m}(t) = \sum_{l=1}^{N} \hat{m}(\sigma_l) e^{i\sigma_l \Omega_0 t} \quad (4)
$$

and similar expressions for ̂ṁ(t), ̂ṁ_s(t) and ̂ṅ_s(t). Time derivatives are replaced by iσᵢΩ₀ and solutions for each coefficient ̂ṁ(σᵢ) can be found independently in terms of the coefficient ̂φ₀(σᵢ) of the tidal potential.

The dynamic equations form a system of coupled linear algebraic
equations written concisely in matrix form as

$$
\mathbf{M}(\sigma_l) \cdot \mathbf{x}(\sigma_l) = \mathbf{y}(\sigma_l) \hat{\phi}_0(\sigma_l) + \mathbf{y}_L(\sigma_l) \hat{\phi}_L(\sigma_l) + \mathbf{y}_h(\sigma_l) \hat{h}(\sigma_l), \quad (5)
$$

where the vector **x** = [**m**, **m**<sub>f</sub>, **m**<sub>s</sub>, **n**<sub>s</sub>]<sup>T</sup> and **φ**<sub>0</sub>(**σ**<sub>l</sub>), **φ**<sub>L</sub>(**σ**<sub>l</sub>) and **h**(**σ**<sub>l</sub>) are the three sources of excitation of the wobble motion, all in the frequency domain: the tidal and ocean loading potentials, as well as the changes in the relative angular momentum of the ocean represented by **h**(**σ**<sub>l</sub>). Among these, the tidal potential is by far the most important. The 4 × 4 matrix **M** and the four-components vectors **y**, **y**<sub>L</sub> and **y**<sub>h</sub> describe the rotational response of the Earth to these three types of excitation. It includes the physical properties of the Earth characterized by the principal moments of inertia of each region, the 12 compliances, and the coupling constants **K**<sub>CMB</sub> and **K**<sub>ICB</sub>. The complete expression of **M** and **y** are

given in Mathews et al. (1991a) (and in Mathews et al. (2002) for the inclusion of the coupling constants). The vectors **y**<sub>L</sub> and **y**<sub>h</sub> are given by: **y**<sub>L</sub>(**σ**<sub>l</sub>) = [(1 + **σ**<sub>l</sub>)(**τ** − **χ**), −**σ**<sub>l</sub>**η**, −**σ**<sub>l</sub>**λ**, 0]<sup>T</sup>, and **y**<sub>h</sub>(**σ**<sub>l</sub>) = [−(1 + **σ**<sub>l</sub>), 0, 0, 0]<sup>T</sup>, where **τ** = **a**<sup>5</sup> **Ω**<sub>0</sub><sup>2</sup>/(3GA), **a** is the mean radius of the Earth, and **G** is the gravitational constant.

The solution of the wobble in the frequency domain is

$$
\begin{equation}
\begin{split}
\hat{m}(\sigma_l) = {}& [\mathbf{M}^{-1}(\sigma_l) \cdot \mathbf{y}(\sigma_l)]_1 \hat{\phi}_0(\sigma_l) + [\mathbf{M}^{-1}(\sigma_l) \cdot \mathbf{y}_L(\sigma_l)]_1 \hat{\phi}_L(\sigma_l) \\
& + [\mathbf{M}^{-1}(\sigma_l) \cdot \mathbf{y}_h(\sigma_l)]_1 \frac{\hat{h}(\sigma_l)}{A \Omega_0},
\end{split}
\tag{6}
\end{equation}
$$

where the notation [ . . . ]<sub>l</sub> indicates the first component of the vector.

Eqs (4) and (6) allow us to compute the wobble ̃m(t) for given excitation sources. This wobble is kinematically related to the nutation of the figure axis, that is, the variations in the direction of this axis in space. The latter is expressed in terms of two variables: the ‘nutation in obliquity’ (Δε), defined as the angle between the figure axis and the pole of the ecliptic, and the ‘nutation in longitude’ (Δψ) corresponding to the motion of the vernal point (the intersection between the equator and the ecliptic) along the ecliptic. This 2-D motion is represented by the complex nutation variable ̃η(t) defined as

$$
\tilde{\eta}(t) = \Delta\psi(t) \sin \epsilon_0 + i\Delta\epsilon(t), \quad (7)
$$

where $\epsilon_0 = 23^\circ 26'14.4''$ is the mean obliquity of the equator at J2000
(i.e. on 2000 January 1). The relation between the wobble $\hat{m}(t)$ and
the nutation $\tilde{\eta}(t)$, referred to as the Euler kinematic relation, is given
by (e.g. Moritz & Mueller 1987; Mathews & Bretagnon 2003)

$$
i \frac{d\tilde{\eta}(t)}{dt} = \Omega_0 \tilde{m}(t) e^{i\Omega_0 t}. \quad (8)
$$

This expression allows us to obtain a prediction of the nutation ˜η(t)
based on the wobble solution ˜m(t) of our model. This prediction
can be compared to the actual observations and the parameters that
enter the nutation model can then be estimated.

2.1.1 *Secular terms*

Let us first consider the particular case of a periodic wobble with
frequency σᵢ = -1, namely ̂ṁ(t) = ̂ṁ(-1)e⁻iσ₀t. In this case,
eq. (8) becomes: d̃η/dt = -iΩ₀̂ṁ(-1). As the right-hand side is
constant, the solution for ̃η(t) is linear in time. This is the secular part
of the motion which, for the component in longitude, corresponds
to the precession. This secular term can be written in the form

$$
\tilde{\eta}^{\text{sec}}(t) = \{P \sin \epsilon_0 + i \Delta \dot{\epsilon}\}(t - t_0) + c_{\psi} \sin \epsilon_0 + i c_{\epsilon}, \quad (9)
$$

where *P* and *Δ*Δ*ε* are constants called the precession and obliquity
rates, respectively, and *c*<sub>*ψ*</sub> and *c*<sub>*ε*</sub> are offset constants characterizing
*tilde{\eta}^{\text{sec}}* at the initial time *t*<sub>0</sub>. Here, we fix *t*<sub>0</sub> to J2000. The constants *P*
and *Δ*Δ*ε* depend directly on *H* ≡ (C − *A*)/*C*, and for this reason *H*
is referred to as the ‘precession constant’.

2.1.2 *Periodic terms*

For a wobble $\tilde{m}(t)$ given by eq. (4) with $\sigma_l \neq -1, \forall l$, eq. (8) gives the associated nutations

$$
\tilde{\eta}^{\text{per}}(t) = \sum_{l=1}^{N} \hat{\eta}_0(\sigma_l) e^{i(1+\sigma_l)\Omega_0 t} \quad (10)
$$

with

$$
\hat{\eta}_0(\sigma_l) = - \frac{\hat{m}(\sigma_l)}{(1 + \sigma_l)}. \qquad (11)
$$
Thus, the nutation associated with a periodic wobble of frequency $\sigma_i \neq -1$ and amplitude $\hat{m}(\sigma_i)$ is also periodic with a frequency $(1 + \sigma_i)$ and an amplitude $-\hat{m}(\sigma_i)/(1 + \sigma_i)$. The difference of one cycle per day in the angular frequencies of the wobble and the associated nutation results from the different reference frames in which they are defined: the wobble is defined in the terrestrial frame, itself rotating with an angular velocity of one cycle per day with respect to the inertial frame in which the nutation is defined. With eqs (6), (10) and (11), a prediction for the nutation variable $\tilde{\eta}(t)$ can be computed.

In order to emphasize the role played by deformations and internal couplings on the nutation motion, the solution is classically given in terms of the nutation $\tilde{\eta}_R$ that the Earth would have if it were completely rigid. The relation between the rigid-Earth nutation and the tidal potential is given by (Mathews et al. 1991a)

$$ \hat{\eta}_0^R(\sigma_i) = \frac{e_R}{(\sigma_i - e_R)(1 + \sigma_i)} \hat{\phi}_0(\sigma_i), \quad (12) $$

where $e_R$ is the dynamic ellipticity [i.e. $(C - A)/A$] of the rigid-Earth model. Rigid-Earth nutations $\hat{\eta}_0^R(\sigma_i)$ can be computed from the very precise ephemerides of the solar system bodies. Several rigid Earth nutation models are available; here, we use the model REN2000 computed by Souchay et al. (1999).

Written in terms of the rigid-Earth nutation, the prediction for the nutation of the non-rigid Earth becomes

$$ \hat{\eta}_0(\sigma_i) = \text{TF}^G(\sigma_i) \hat{\eta}_0^R(\sigma_i) + \text{TF}^L(\sigma_i) \hat{\phi}_L(\sigma_i) + \text{TF}^h(\sigma_i) \frac{\hat{h}(\sigma_i)}{A\Omega_0}, \quad (13) $$

where the coefficients $\text{TF}^G(\sigma_i)$, $\text{TF}^L(\sigma_i)$ and $\text{TF}^h(\sigma_i)$ are called ‘transfer functions’ and are given by

$$ \text{TF}^G(\sigma_i) = \frac{e_R - \sigma_i}{e_R} [\mathbf{M}^{-1}(\sigma_i) \cdot \mathbf{y}(\sigma_i)]_1 \quad (14a) $$

$$ \text{TF}^L(\sigma_i) = - \frac{[\mathbf{M}^{-1}(\sigma_i) \cdot \mathbf{y}_L(\sigma_i)]_1}{1 + \sigma_i} \quad (14b) $$

$$ \text{TF}^h(\sigma_i) = - \frac{[\mathbf{M}^{-1}(\sigma_i) \cdot \mathbf{y}_h(\sigma_i)]_1}{1 + \sigma_i} \quad (14c) $$

$\text{TF}^G(\sigma_i)$ accounts for all deformations and internal couplings characterizing the non-rigidity of the Earth. $\text{TF}^L(\sigma_i)$ and $\text{TF}^h(\sigma_i)$ describe the nutational response of the Earth to the loading potential and relative angular momentum of the ocean, respectively. Writing the nutation motion in the form (13) allows for a separation of two different problems: the modelling of the rigid-Earth nutation, which is a purely astronomical problem that can be solved from celestial mechanics, and the modelling of the transfer functions, a geophysical problem that depends on parameters characterizing the Earth’s interior.

The transfer function $\text{TF}^G(\sigma_i)$ can be rewritten in the form (Mathews et al. 2002)

$$ \text{TF}^G(\sigma) = \frac{e_R - \sigma}{1 + e_R} \frac{H}{H_R} \left[ 1 + (1 + \sigma) \sum_{i=1}^{4} \frac{N_i}{\sigma - \sigma_i} \right], \quad (15) $$

where $H$ and $H_R$ are the precession constants of the Earth and of the rigid-Earth model, respectively. The $\sigma_i$ are four resonance frequencies, corresponding to the frequencies of the free modes allowed by the system (5). These corresponds to the frequencies for which $\text{Det}[\mathbf{M}(\sigma)] = 0$. They are: the Chandler wobble (CW), the FCN, the FICN and the inner core wobble (ICW). Eq. (15) shows that the nutational response of the Earth to an external forcing

depends on the forcing frequency in a resonant way: the closer is the forcing frequency to that of one of the free modes, the largest is the amplitude of the corresponding nutation. To the first order, the frequencies of the free modes are given, in the terrestrial frame (Mathews et al. 2002)

$$ \sigma_{CW} \approx \frac{A}{A_m} (e - \kappa) \quad (16a) $$

$$ \sigma_{FCN} \approx -1 - \left(1 + \frac{A_f}{A_m}\right) \left(e_f - \beta + K_{CMB} + K_{ICB} \frac{A_s}{A_f}\right) \quad (16b) $$

$$ \sigma_{FICN} \approx -1 + \left(1 + \frac{A_s}{A_m}\right) (\alpha_2 e_s + \nu - K_{ICB}) \quad (16c) $$

$$ \sigma_{ICW} \approx e_s(1 - \alpha_2), \quad (16d) $$

where $e$, $e_f$ and $e_s$ are the dynamic ellipticities of the Earth, the outer core and the inner core, respectively, and $\alpha_2$ is a parameter introduced by Mathews et al. (1991a) in the computation of the inertial and gravitational torque on the inner core. If there was no dissipation inside the Earth, the free mode frequencies would be real. In our model, dissipation occurs through mantle anelasticity and couplings at the CMB and ICB and the free modes have complex frequencies. Two of the free modes have frequencies that are close to diurnal and therefore close to those of the tidal potential, the FCN and FICN. These modes are then of major importance for the nutation motion.

The oceanic effects on nutation are described in eq. (13) by $\hat{\phi}_L(\sigma_i)$ and $\hat{h}(\sigma_i)$. Following Mathews et al. (2002), we only consider the effects from ocean tides and we do not take into account other dynamic effects of the oceanic circulation. The oceanic tidal angular momentum (OTAM) can be separated into two terms, called respectively the ‘height’ ($\hat{H}^L$) and ‘current’ ($\hat{H}^c$) terms, which are directly related to $\hat{\phi}_L$ and $\hat{h}$ by $\hat{H}^L = -A\tau\Omega_0\hat{\phi}_L$ and $\hat{H}^c = \hat{h}$. The OTAM can be computed from models of the main ocean tides (see e.g. Chao et al. 1996). For the other frequencies, the OTAM is computed from an interpolation procedure introduced by Wahr & Sasao (1981). Further details on the computation of the OTAM can be found in Koot et al. (2008).

### 2.1.3 Poisson terms

Although we approximated $\tilde{\phi}$ in (3) as a sum of periodic terms, the current precision of nutation observations is such that the effect of small departures from a purely periodic tidal potential can be detected. Very small linear corrections, $\hat{\phi}_1(\sigma_i)$, called ‘Poisson terms’, have to be introduced

$$ \tilde{\phi}(t) = \sum_{i=1}^{N} \{\hat{\phi}_0(\sigma_i) + \hat{\phi}_1(\sigma_i)t\} e^{i\sigma_i\Omega_0 t} \quad (17) $$

These small $\hat{\phi}_1(\sigma_i)$ terms are assumed to give rise to a nutation motion of the form

$$ \hat{\eta}^P(t) = \sum_{i=1}^{N} \hat{\eta}_1(\sigma_i) t e^{i(\sigma_i+1)\Omega_0 t} \quad (18) $$

Folgueira et al. (2007) showed that, to the first order, $\hat{\eta}_1(\sigma_i)$ is given by

$$ \hat{\eta}_1(\sigma_i) = \text{TF}^G(\sigma_i) \hat{\eta}_1^R(\sigma_i), \quad (19) $$

where $\hat{\eta}_1^R(\sigma_i)$ are the Poisson terms of the rigid-earth model, which are directly related to the tidal potential terms $\hat{\phi}_1(\sigma_i)$.
## 2.2 Free nutations

In addition to their amplification effect on the forced nutations, the free nutations FCN and FICN, excited by other means, participate in the total nutation motion. Their amplitude is expected to be time-dependent, depending on the time-history of their excitation source. While the free nutation arising from the FICN mode has not been observed yet, the FCN mode can actually be observed in the nutation data (e.g. Herring *et al.* 2002), with an amplitude approximately five orders of magnitude smaller than that of the main forced nutation (less than 0.5 mas over the period 1979–2010). The FCN mode is presumed to be excited, at least partly, by diurnal variations in atmospheric pressure (Sasao & Wahr 1981; Lambert 2006).

Unfortunately, models of the atmospheric angular momentum (AAM), computed from general circulation models, are not precise enough in the diurnal band to allow for a modelling of the excited FCN mode (Lambert 2006). For this reason, we model the free nutation empirically, as a pseudo-periodic term with an amplitude varying in time

$$ \tilde{\eta}^{fcn}(t) = -i a_{FCN}(t) e^{iF_{FCN}(t-t_0)}, \quad (20) $$

where $F_{FCN}$ is the frequency of the FCN mode in the celestial frame, that is, $F_{FCN} = [1 + \text{Re}(\sigma_{FCN})] \Omega_0$. Following Herring *et al.* (2002), we choose to model the time-varying amplitude $a_{FCN}(t)$ as a piecewise linear function, that is, a function which is linear between given time nodes. The time nodes are chosen from a preliminary analysis to be at 1979, 1984, 1990, 1997, 2001 and 2010.

## 2.3 Corrections to the model

The nutation model is based on classical mechanics and does not take into account relativistic effects such as the geodesic nutation. This effect has been computed by Fukushima (1991) and contribute a change of $-0.0304$ and $-0.0004$ mas on the prograde annual and semi-annual nutations, respectively, and $0.0304$ and $0.0004$ mas on the retrograde annual and semiannual nutations, respectively. We have added these corrections to our model.

Additionally, as the equations are developed to the first order, the model does not take into account non-linear effects. Second order effects have been computed by Lambert & Mathews (2006), who showed that these effects almost compensate each other so that the only non-negligible contributions are on the 18.6 yr nutation and are $(0.0070 - i 0.0003)$ mas for the prograde component and $(-0.0070 - i 0.0013)$ mas for the retrograde component. These corrections are also added to our model.

## 2.4 Parameters to be estimated

Among the parameters of the precession/nutation model, some have a large influence on the motion and can thus be estimated from the observations. These parameters are listed in Table 1. Parameters that have a lower influence on nutation are assigned fixed numerical values and are listed in Table 2.

### 2.4.1 Geophysical parameters

The focus of our study is on the geophysical parameters. Among them, the precession constant $H$ [or equivalently the dynamic ellipticity $e$ which is related to $H$ by $H = e/(1 + e)$] determines the precession rate and has the largest influence on the precession/nutation motion. Additionally, it also directly influences the nutation transfer

<table><caption>Table 1. Parameters to be estimated from the nutation observations.</caption><thead><tr><th>Symbol</th><th>Definition</th></tr></thead><tbody><tr><td colspan="2"><strong>Geophysical parameters</strong></td></tr><tr><td><i>e</i></td><td>Dynamic ellipticity of the whole Earth</td></tr><tr><td><i>e</i><sub>f</sub> + Re(K<sub>CMB</sub>)</td><td>Dynamic ellipticity of the fluid core + CMB coupling constant (real part)</td></tr><tr><td><i>κ</i><sup>el</sup></td><td>Elastic compliance of the whole Earth</td></tr><tr><td><i>γ</i><sup>el</sup></td><td>Elastic compliance of the fluid core</td></tr><tr><td>Im(K<sub>CMB</sub>)</td><td>CMB coupling constant (imaginary part)</td></tr><tr><td>Re(K<sub>ICB</sub>)</td><td>ICB coupling constant (real part)</td></tr><tr><td>Im(K<sub>ICB</sub>)</td><td>ICB coupling constant (imaginary part)</td></tr><tr><td colspan="2"><strong>Secular terms parameters</strong></td></tr><tr><td>Δ<i>ε</i></td><td>Obliquity rate</td></tr><tr><td><i>c</i><sub>ψ</sub></td><td>Constant offset in longitude</td></tr><tr><td><i>c</i><sub>ε</sub></td><td>Constant offset in obliquity</td></tr><tr><td colspan="2"><strong>Amplitude of the FCN</strong></td></tr><tr><td>Re(a<sub>FCN</sub>)</td><td>Real parts of the complex amplitudes</td></tr><tr><td>Im(a<sub>FCN</sub>)</td><td>Imaginary parts of the complex amplitudes</td></tr><tr><td colspan="2"><strong>Atmospheric contribution to the amplitude of the prograde annual term</strong></td></tr><tr><td>Re(a<sub>ap</sub>)</td><td>Real part of the contribution</td></tr><tr><td>Im(a<sub>ap</sub>)</td><td>Imaginary part of the contribution</td></tr><tr><td colspan="2"><strong>Probabilistic modelling parameter</strong></td></tr><tr><td>σ<sub>M</sub></td><td>Standard deviation of the Gaussian modelling uncertainty</td></tr></tbody></table>

function (see eq. 15) so that it affects the amplitude of all of the forced nutations as well.

The most important parameters for the nutation motions are those that determine the frequencies and damping of the three main free modes ($\sigma_{CW}$, $\sigma_{FCN}$ and $\sigma_{FICN}$) and the strengths of the associated resonances ($N_1, N_2, N_3$ in eq. 15). The expressions for $N_1$ and $N_2$ are (Mathews *et al.* 2002)

$$ N_1 \approx -\frac{A}{A_m} \left( \frac{e - \kappa}{e} \right) \quad (21a) $$

$$ N_2 \approx \frac{A_f}{A_m} \left( \frac{e - \gamma}{e} \right) \quad (21b) $$

the expression of $N_3$ is not given as it does not introduce other geophysical parameters.

Eqs (16a), (21a) and (21b) show the importance of the compliances $\kappa$ and $\gamma$. These compliances are complex numbers and can be written as the sum of their (real) value for an elastic earth model and a small complex correction for mantle anelasticity effects: $\kappa = \kappa^{el} + \Delta\kappa^{AE}$ and $\gamma = \gamma^{el} + \Delta\gamma^{AE}$. Following Mathews *et al.* (2002), we fix the value of the anelastic contributions $\Delta\kappa^{AE}$ and $\Delta\gamma^{AE}$ to that computed from an Earth interior model and we estimate the elastic term, namely $\kappa^{el}$ and $\gamma^{el}$.

Eqs (16b) and (16c) show combinations between $K_{CMB}$, $K_{ICB}$, $\beta$, $\nu$, $e_f$ and $e_s$. Following Mathews *et al.* (2002), we choose to fix the values of the compliances (except $\kappa^{el}$ and $\gamma^{el}$) to those computed from an Earth interior model. The compliances for an elastic Earth have been computed by Mathews *et al.* (1991b) for the PREM model (Dziewonski & Anderson 1981) and are listed in Table 2. Small contributions to the compliances from mantle anelasticity can be computed, assuming no bulk dissipation, from the shear wave quality factor $Q_\mu$ of the PREM model. We computed these contributions following the procedure described in Mathews *et al.* (2002) and the numerical values are listed in Table 2. Also, we fix $e_s$ to its hydrostatic value as the whole of the inner core is not expected to depart much from hydrostatic equilibrium (e.g. Defraigne
<table><caption>**Table 2.** Numerical values of other parameters entering the
nutation model.</caption><thead><tr><th>Parameter</th><th>Numerical values</th></tr></thead><tbody><tr><td><strong>Principal moments of inertia</strong></td><td></td></tr><tr><td>A</td><td>8.0115 × 10<sup>37</sup> kg m<sup>2</sup></td></tr><tr><td>A<sub>f</sub></td><td>9.0583 × 10<sup>36</sup> kg m<sup>2</sup></td></tr><tr><td>A<sub>s</sub></td><td>5.8531 × 10<sup>34</sup> kg m<sup>2</sup></td></tr><tr><td>e<sub>s</sub></td><td>2.422 × 10<sup>-3</sup></td></tr><tr><td><strong>Elastic compliances</strong></td><td></td></tr><tr><td>θ<sup>el</sup></td><td>6.794 × 10<sup>-6</sup></td></tr><tr><td>ξ<sup>el</sup></td><td>2.222 × 10<sup>-4</sup></td></tr><tr><td>β<sup>el</sup></td><td>6.160 × 10<sup>-4</sup></td></tr><tr><td>α<sup>el</sup></td><td>-7.536 × 10<sup>-5</sup></td></tr><tr><td>ζ<sup>el</sup></td><td>4.964 × 10<sup>-9</sup></td></tr><tr><td>δ<sup>el</sup></td><td>-4.869 × 10<sup>-7</sup></td></tr><tr><td>ν<sup>el</sup></td><td>7.984 × 10<sup>-5</sup></td></tr><tr><td>χ</td><td>1.073 × 10<sup>-3</sup></td></tr><tr><td>η</td><td>1.940 × 10<sup>-3</sup></td></tr><tr><td>λ</td><td>0</td></tr><tr><td><strong>Anelastic contributions to the compliances</strong></td><td></td></tr><tr><td>ΔK<sup>AE</sup></td><td>(13 + i 5) × 10<sup>-6</sup></td></tr><tr><td>Δγ<sup>AE</sup></td><td>(22 + i 9) × 10<sup>-6</sup></td></tr><tr><td>Δθ<sup>AE</sup></td><td>(3.7 + i 1.5) × 10<sup>-8</sup></td></tr><tr><td>Δξ<sup>AE</sup></td><td>(2.5 + i 1.0) × 10<sup>-6</sup></td></tr><tr><td>Δβ<sup>AE</sup></td><td>(7.4 + i 3.0) × 10<sup>-6</sup></td></tr><tr><td>Δα<sup>AE</sup></td><td>(1 + i 0.4) × 10<sup>-8</sup></td></tr><tr><td>Δζ<sup>AE</sup></td><td>(25.1 + i 1.2) × 10<sup>-11</sup></td></tr><tr><td>Δδ<sup>AE</sup></td><td>0</td></tr><tr><td>Δν<sup>AE</sup></td><td>0</td></tr><tr><td><strong>Parameter τ</strong></td><td></td></tr><tr><td>τ</td><td>3.481 × 10<sup>-3</sup></td></tr><tr><td><strong>Parameters α<sub>i</sub></strong></td><td></td></tr><tr><td>α<sub>1</sub></td><td>0.9463</td></tr><tr><td>α<sub>2</sub></td><td>0.8294</td></tr><tr><td>α<sub>3</sub></td><td>0.0537</td></tr><tr><td><strong>Precession</strong></td><td></td></tr><tr><td>H<sub>R</sub></td><td>0.0032737668</td></tr><tr><td>P<sub>R</sub></td><td>50 384.565 mas a<sup>-1</sup></td></tr><tr><td>dP<sub>R</sub>/dH<sub>R</sub></td><td>15 397 060 mas a<sup>-1</sup></td></tr><tr><td>P<sub>NR</sub></td><td>-0.2015 mas a<sup>-1</sup></td></tr></tbody></table>

et al. 1996). The parameters that are to be estimated are then $e_f$ + Re(K<sub>CMB</sub>), Im(K<sub>CMB</sub>), Re(K<sub>ICB</sub>) and Im(K<sub>ICB</sub>). Their values are directly dependent on the values adopted for the fixed parameters. For example, if $e_s$ departs slightly from its hydrostatic equilibrium value, this departure is absorbed in the parameter Re(K<sub>ICB</sub>). Similarly, errors in the computation of the compliances $\beta$ and $\nu$ directly affect the estimate of the coupling constants.

The numerical values of A, A<sub>f</sub>, A<sub>s</sub> and e<sub>s</sub>, as well as the parameters $\alpha_i$ are taken as those computed by Mathews et al. (1991b) from the PREM interior model (Dziewonski & Anderson 1981). They are given in Table 2.

### 2.4.2 Secular terms parameters

The secular component of the motion, given by eq. (9), depends on the parameters P, Δe, c<sub>ψ</sub> and c<sub>ε</sub>. The latter three are new parameters to estimate and are listed in Table 1. The precession rate P can be expressed in terms of the precession constant H in the form

$$P(H) = P_{R}(H) + P_{NR}, \qquad (22)$$

where P<sub>R</sub>(H) is the precession rate of the rigid Earth, which depends directly on H, and P<sub>NR</sub> is the contribution to the precession rate due to non-rigid effects and is independent of H. The value of H is close to that estimated from rigid Earth models (H<sub>R</sub>) so that P<sub>R</sub>(H) can be expanded as

$$P_{R}(H) = P_{R}(H_{R}) + \frac{dP_{R}}{dH_{R}}(H - H_{R}). \qquad (23)$$

Numerical values for H<sub>R</sub>, P<sub>R</sub>(H<sub>R</sub>) and dP<sub>R</sub>/dH<sub>R</sub> from Bretagnon et al. (1998) are listed in Table 2 along with the value of P<sub>NR</sub> computed by Mathews et al. (2002).

### 2.4.3 Free mutation parameters

The amplitudes a<sub>FCN</sub> = {a<sub>FCN</sub>(t<sub>i</sub>)<sub>i=1</sub><sup>6</sup> of the FCN mode at the time of the selected time nodes {t<sub>i</sub>}<sub>i=1</sub><sup>6</sup> are parameters of the model to be estimated from the observations.

### 2.4.4 Atmospheric contributions parameters

Variations in the AAM, mainly due to the diurnal cycle in the heating of the atmosphere by the Sun, also induce forced nutations whose amplitudes are typically five orders of magnitude smaller than those driven by the tidal potential. In principle, these contributions could be computed from AAM models (Bizouard et al. 1998) but, as shown by Yseboodt et al. (2002), the results strongly depend on the choice of the atmospheric model. The estimation from AAM models is thus not precise enough to be added to our model. However, atmospheric effects influence mainly the annual prograde nutation (Bizouard et al. 1998); we thus follow Mathews et al. (2002) and incorporate them as a correction to the annual prograde nutation. This correction a<sub>ap</sub> is an additional parameter of the model to be estimated.

## 3 NUTATION DATA

Nutations are computed from VLBI observations. This technique consists in observing very distant galactic objects (quasars) from an array of Earth based large radio-telescopes. After each observing session, the signals recorded at each telescope are compared and differences between arrival times (delays) are inferred by correlations. These observed delays can be modelled in terms of several parameters such as the geographic location of the telescopes, the position of the quasars, the propagation delays due to the atmosphere, and the orientation of the Earth in space (and thus, in particular, the precession/nutation motion). The delays are then inverted to obtain an estimation of these parameters.

VLBI-delay data are processed by several analysis centres in order to get estimations of these parameters, including the precession/nutation angles of interest in our study. Each centre uses different analysis strategies and softwares to compute the nutation offsets. In this study, we use nutation angle time-series produced by three such centres: the NASA Goddard Space Flight Center (GSFC), the Paris Observatory (OPA) and the Institute of Applied Astronomy (IAA). We choose these centres because they present important differences in their analysis strategies. While both the GSFC and OPA centres use the Calc/Solv software, the IAA centre uses the OCCAM software. One of the characteristics of the OPA solution is that they use a different celestial reference frame that is defined by 247 stable sources selected by Feissel-Vernier et al. (2006), while the other centres use the conventional 212 sources that defined the
International Celestial Reference Frame (ICRF). All three nutation time-series cover the time interval from 1979 up to the present. The specific data sets used here are available on the International VLBI Service (IVS) web site: http://ivscc.gsfc.nasa.gov and labelled ‘gsf2008a.eops’, ‘opa2009b.eops’, and ‘iaa2007a.eops’.

# 4 INVERSION STRATEGY

The strategy that we use for the inversion of the nutation observations has been described in details in Koot et al. (2008). We summarize here its main features.

## 4.1 Inversion in the time-domain

As the gravitational forcing is mainly periodic, nutation models are often developed in the frequency domain. However, nutation observations are available as time-series. In order to compare the model to the data, there are two possibilities: either (1) to estimate the amplitudes of the most important periodic terms in the time-domain data and compare these amplitudes to the frequency-domain nutation model or (2) to expand the nutation model in the time domain and compare it directly to the time-domain observations. The first method was used by Herring et al. (2002) while in Koot et al. (2008) we developed the second one.

The main advantage of using the time-domain method is that it allows us to use all the available data, which is not the case in the frequency domain where only the amplitudes of some 21 periodic terms can be extracted from the data (Herring et al. 2002).

The time-domain nutation model that we use is composed of the secular (eq. 9), periodic (eqs 10, 13), Poisson (eq. 18) and free FCN (eq. 20) motions,

$$ \tilde{\eta}(t, \theta) = \tilde{\eta}^{\text{sec}}(t, \theta) + \tilde{\eta}^{\text{per}}(t, \theta) + \tilde{\eta}^{P}(t, \theta) + \tilde{\eta}^{\text{fcn}}(t, \theta), \quad (24) $$

where the symbol $\theta$ represents the parameters of the model that are to be estimated. The time-domain model given by (24) is more complete than the one used in Koot et al. (2008) where the Poisson terms were not taken into account. Also, in Koot et al. (2008), the time-variable amplitude of the FCN mode was estimated separately from a preliminary analysis and was removed from the data before the final inversion. Here, we estimate it jointly with the other geophysical parameters, a more self-consistent approach.

## 4.2 Bayesian inversion method

### 4.2.1 The Bayesian framework

In a Bayesian inversion framework (e.g. Gregory 2005; Tarantola 2005), knowledge of the model parameters ($\theta$) is inferred from independent probability density functions (pdf) that describe the information provided by the observed data ($\mathbf{d}$), the data-independent prior information on the parameters ($\pi$) and the model predicted data ($\mathbf{d}^*$).

The probability for the parameters given the observed data and the prior information is called the ‘posterior’ pdf and can be written as

$$ p(\theta|\mathbf{d}) \propto L(\mathbf{d}|\theta) \pi(\theta), \quad (25) $$

where $L(\mathbf{d}|\theta)$ represents the likelihood of observing the data $\mathbf{d}$ given the parameter values $\theta$. We construct $L$ by assuming that the observed data are related to the model prediction by

$$ \mathbf{d} = \mathbf{d}^* + \mathbf{e}_d, \quad (26) $$

where $\mathbf{e}_d$ represents the uncertainty associated with the observed data and is chosen to be a normal distribution of zero mean and variance $\sigma^2$. The model predictions $\mathbf{d}^*$ are related to the parameters by

$$ \mathbf{d}^* = \mathbf{M}(\theta) + \mathbf{e}_M, \quad (27) $$

where $\mathbf{M}(\theta)$ is the model, also called the ‘forward problem’, and $\mathbf{e}_M$ characterizes model uncertainty, independent of data. We choose the pdf representing $\mathbf{e}_M$ to be a normal distribution of zero mean and variance $\sigma_M^2$. In this case, it can be shown (e.g. Gregory 2005) that the likelihood for observing data $\mathbf{d}$ given parameter values $\theta$ is

$$ L(\mathbf{d}|\theta) \propto \left[ \prod_{i=1}^{N} (\sigma_i^2 + \sigma_M^2)^{-1/2} \right] \exp \left\{ -\frac{1}{2} \sum_{i=1}^{N} \frac{[d_i - M_i(\theta)]^2}{(\sigma_i^2 + \sigma_M^2)} \right\}. \quad (28) $$

From the posterior pdf, we can get the pdf for one specific parameter $\theta_i$ by computing the marginal of the posterior defined by

$$ p(\theta_i|\mathbf{d}) = \int p(\theta|\mathbf{d}) d\theta_1 \dots d\theta_{i-1} d\theta_{i+1} \dots d\theta_n. \quad (29) $$

From this marginal, an estimated value for the parameter $\theta_i$ can be obtained by computing its mean or the value for which the probability is the highest (the maximum a posteriori, MAP, solution). Confidence intervals (CI) corresponding to a given probability can also be estimated from the marginal in order to get an estimation of the error on the parameter.

### 4.2.2 Posterior used in our study

In our case, the data are two time series: $\mathbf{d}^R = \{\Delta\psi_i \sin\epsilon_0\}_{i=1}^N$ for the nutation in longitude and $\mathbf{d}^I = \{\Delta\epsilon_i\}_{i=1}^N$ for the nutation in obliquity, corresponding to the time indices $\{t_i\}_{i=1}^N$. We assume that the uncertainties in the measurement of each are independent and that they are given by normal distributions of zero mean and standard deviations $\sigma^R = \{\sigma_i^\psi \sin\epsilon_0\}_{i=1}^N$ and $\sigma^I = \{\sigma_i^\epsilon\}_{i=1}^N$, respectively. In this case we have

$$ L(\mathbf{d}|\theta) = L(\mathbf{d}^R|\theta) L(\mathbf{d}^I|\theta) \quad (30) $$

and the likelihood functions $L(\mathbf{d}^R|\theta)$ and $L(\mathbf{d}^I|\theta)$ can both be computed from eq. (28), which gives

$$ L(\mathbf{d}|\theta) \propto \left( \prod_{i=1}^{N} \left[ (\sigma_i^R)^2 + \sigma_M^2 \right]^{-1/2} \left[ (\sigma_i^I)^2 + \sigma_M^2 \right]^{-1/2} \right) \times \exp \left\{ -\frac{1}{2} \sum_{i=1}^{N} \left[ \frac{(d_i^R - M_i^R(\theta))^2}{(\sigma_i^R)^2 + \sigma_M^2} + \frac{(d_i^I - M_i^I(\theta))^2}{(\sigma_i^I)^2 + \sigma_M^2} \right] \right\}, \quad (31) $$

where $M_i^R(\theta) = \text{Re}[\tilde{\eta}(t_i, \theta)]$ and $M_i^I(\theta) = \text{Im}[\tilde{\eta}(t_i, \theta)]$, with $\tilde{\eta}(t_i, \theta)$ given by eq. (24).

In writing eq. (31), we make the hypothesis that the standard deviation characterizing the model uncertainty $\sigma_M$ is the same for nutation in longitude and obliquity. As this parameter is unknown, it is considered as an additional parameter to be estimated. From eq. (25), the posterior pdf is then given by

$$ p(\theta, \sigma_M|\mathbf{d}) \propto L(\mathbf{d}|\theta) \pi(\theta, \sigma_M), \quad (32) $$

where $L(\mathbf{d}|\theta)$ is given by eq. (31) and $\pi(\theta, \sigma_M)$ is the prior pdf on the parameters $\theta$ and $\sigma_M$. Detail on the prior that we use are given in Koot et al. (2008). Note that, by our choice of prior, we impose the following constraints on the coupling constants:
<table><caption>**Table 3.** Estimated values and errors for the nutation model parameters computed from the GSFC, OPA and IAA nutation time-series, as well as for the joint inversion.</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Units</th><th colspan="6">Data sets</th></tr><tr><th>Joint inversion<br>GSFC, OPA, IAA</th><th>GSFC</th><th>OPA</th><th>IAA</th><th>GSFC1<br>(1990-2009.3)</th><th>MHB</th></tr></thead><tbody><tr><td>e</td><td>10<sup>-3</sup></td><td>3.2845482 ± 1</td><td>3.2845481 ± 2</td><td>3.2845477 ± 2</td><td>3.2845488 ± 2</td><td>3.2845484 ± 2</td><td>3.2845479 ± 13</td></tr><tr><td>e<sub>f</sub> + Re(K<sub>CMB</sub>)</td><td>10<sup>-3</sup></td><td>2.6753 ± 2</td><td>2.6751 ± 4</td><td>2.6755 ± 4</td><td>2.6753 ± 4</td><td>2.6751 ± 3</td><td>2.6680 ± 21</td></tr><tr><td>κ<sup>el</sup></td><td>10<sup>-3</sup></td><td>1.038 ± 2</td><td>1.038 ± 3</td><td>1.035 ± 3</td><td>1.042 ± 3</td><td>1.038 ± 4</td><td>1.034 ± 10</td></tr><tr><td>γ<sup>el</sup></td><td>10<sup>-3</sup></td><td>1.9649 ± 3</td><td>1.9647 ± 5</td><td>1.9652 ± 5</td><td>1.9646 ± 5</td><td>1.9647 ± 5</td><td>1.966 ± 2</td></tr><tr><td>Im(K<sub>CMB</sub>)</td><td>10<sup>-5</sup></td><td>-1.78 ± 2</td><td>-1.80 ± 3</td><td>-1.78 ± 3</td><td>-1.76 ± 4</td><td>-1.79 ± 3</td><td>-1.85 ± 15</td></tr><tr><td>Re(K<sub>ICB</sub>)</td><td>10<sup>-3</sup></td><td>1.01 ± 2</td><td>0.99 ± 4</td><td>1.03 ± 3</td><td>1.00 ± 4</td><td>1.02 ± 4</td><td>1.11 ± 11</td></tr><tr><td>Im(K<sub>ICB</sub>)</td><td>10<sup>-3</sup></td><td>-1.09 ± 3</td><td>-1.06 ± 5</td><td>-1.09 ± 6</td><td>-1.10 ± 6</td><td>-1.04 ± 5</td><td>-0.78 ± 14</td></tr><tr><td>Re(a<sub>ap</sub>)</td><td>μ<sub>as</sub></td><td>0 ± 4</td><td>-10 ± 6</td><td>12 ± 7</td><td>-1 ± 8</td><td>-9 ± 7</td><td>-10</td></tr><tr><td>Im(a<sub>ap</sub>)</td><td>μ<sub>as</sub></td><td>107 ± 4</td><td>104 ± 6</td><td>108 ± 7</td><td>109 ± 9</td><td>106 ± 7</td><td>108</td></tr><tr><td>Δe</td><td>mas a<sup>-1</sup></td><td>-0.2572 ± 7</td><td>-0.2584 ± 12</td><td>-0.2530 ± 13</td><td>-0.2602 ± 14</td><td>-0.2615 ± 15</td><td>—</td></tr><tr><td>c<sub>ε</sub></td><td>mas</td><td>-6.880 ± 5</td><td>-6.914 ± 7</td><td>-6.845 ± 7</td><td>-6.879 ± 9</td><td>-6.910 ± 7</td><td>—</td></tr><tr><td>c<sub>ψ</sub></td><td>mas</td><td>-41.68 ± 1</td><td>-41.70 ± 2</td><td>-41.64 ± 2</td><td>-41.71 ± 2</td><td>-41.71 ± 2</td><td>—</td></tr><tr><td>σ<sub>M</sub></td><td>mas</td><td>0.121 ± 3</td><td>0.106 ± 5</td><td>0.107 ± 5</td><td>0.132 ± 5</td><td>0.097 ± 5</td><td>—</td></tr></tbody></table>

Notes: Results for the GSFC time-series are shown for the whole period (1979–2009.3) and also for the period 1990–2009.3 (labeled GSFC1). The errors, corresponding to the 99.7 per cent CI, refer to the last written digits. The last column shows the results obtained by Mathews *et al.* (2002) (MHB).

<table><caption>**Table 4.** Estimated values and 99.7 per cent CI for the time-variable amplitude of the FCN mode computed from the inversion of the GSFC data set (units: mas).</caption><thead><tr><th>Time</th><th>Re(a<sub>FCN</sub>)</th><th>Im(a<sub>FCN</sub>)</th></tr></thead><tbody><tr><td>1979</td><td>-0.14 ± 0.33</td><td>-0.33 ± 0.32</td></tr><tr><td>1984</td><td>0.20 ± 0.06</td><td>-0.32 ± 0.05</td></tr><tr><td>1990</td><td>0.03 ± 0.02</td><td>-0.15 ± 0.02</td></tr><tr><td>1997</td><td>-0.07 ± 0.02</td><td>-0.10 ± 0.02</td></tr><tr><td>2001</td><td>0.16 ± 0.02</td><td>0.11 ± 0.02</td></tr><tr><td>2010</td><td>-0.21 ± 0.02</td><td>0.20 ± 0.02</td></tr></tbody></table>

$Im(K_{CMB}) < 0$, $Re(K_{ICB}) > 0$ and $Im(K_{ICB}) < 0$, in accordance with theoretical models of the couplings (Buffett *et al.* 2002; Mathews & Guo 2005). Samples of the posterior pdf are obtained from the Metropolis-Hastings algorithm (Metropolis *et al.* 1953; Hastings 1970).

We note finally that the parameter $\sigma_M$ appearing in eq. (31) cannot be distinguished from a constant additive correction to the measurement errors $\sigma_i^R$ and $\sigma_i^I$. Therefore, our inversion takes into account a possible underestimation of VLBI nutation angles errors, which have been shown by Herring *et al.* (2002) to be too optimistic.

# 5 INVERSION RESULTS

We have carried out individual inversions for each of the three nutation time-series independently, as well as an inversion with all the times series together. As the quality of the data largely improves with time, we have also performed an inversion of the GSFC time-series restricted to the last 20 yr, that is, from 1990 to 2009.3.

## 5.1 Estimated parameters

The estimated values and the 99.7 per cent CI for all parameters, for each inversion, are presented in Table 3. The time-dependent amplitude of the FCN is given in Table 4 for the GSFC data set and shown in Fig. 1 for the three data sets. As expected, for all the parameters, the estimated values obtained from the joint inversion are practically equal to the mean of the individual inversions. The associated errors of the joint inversion are smaller than for

the individual inversions, consistent with the fact that more data were used.

Figs 2 and 3 show the marginal posteriors for $Im(K_{CMB})$, $Re(K_{ICB})$ and $Im(K_{ICB})$. The marginal posteriors for other parameters are similar. Figs 2 and 3 also illustrate that while the errors for the inversions are small, the fact that the pdfs of individual inversions do not completely overlap indicates that these errors are likely underestimated.

In general, our retrieved parameters are in good agreement with the values estimated in Mathews *et al.* (2002). The only discrepancies are for the estimates of $Re(K_{ICB})$, $Im(K_{ICB})$ and $e_f + Re(K_{CMB})$. The differences for the former two and their implication are the focus of Section 6. The discrepancy for $e_f + Re(K_{CMB})$ does not come from our different inversion strategy as we get the same estimation when inverting the frequency data used by Mathews *et al.* (2002). It comes probably from a different value of $Re(\Delta\beta^{AE})$.

Differences between the results presented here and those in Koot *et al.* (2008) are from the inclusion of the Poisson terms which alter the amplitude of the long periods nutations. Other differences are from the contribution of mantle anelasticity on all of the compliances (rather than only on $\kappa$ and $\gamma$). The inclusion of the anelastic contribution to the compliance $\beta$ affects directly our estimates of $e_f + Re(K_{CMB})$ and $Im(K_{CMB})$ because, as already noted in Section 2.4.1, it is the sum ($e_f + K_{CMB} - \beta$) that enters the expression of the FCN frequency (16b).

## 5.2 Residuals

The residuals between the data and the model can be computed by the following estimator:

$$ \bar{r}_i = \int (d_i - \tilde{\eta}(t_i, \boldsymbol{\theta})) p(\boldsymbol{\theta} | \mathbf{d}) \, d\boldsymbol{\theta}. \quad (33) $$

The estimated residuals are shown on Fig. 4 for the GSFC time-series. The residuals for the other data sets are very similar. The weighted root-mean-squares (WRMS) of the residuals, defined by

$$ \text{WRMS}(\mathbf{r}, \boldsymbol{\sigma}) = \sqrt{\sum_i w_i r_i^2}, \quad \text{where } w_i = \frac{1/\sigma_i^2}{\sum_j 1/\sigma_j^2} \quad (34) $$
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Estimated values and 99.7 per cent CI for the time-variable amplitude of the free FCN mode computed from the GSFC, OPA and IAA data sets.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Marginal posterior pdf of Im($K_{CMB}$) for the individual inversion of the three time-series (GSFC, OPA and IAA), the joint inversion (Joint) and for the GSFC data truncated to the period 1990–2009.3 (GSFC1). The blue box is the 99.7 per cent CI obtained by Mathews *et al.* (2002) (MHB).</figcaption></figure>

are 0.12 mas for both the $\Delta\psi \sin\epsilon_0$ and the $\Delta\epsilon$ components, for the GSFC data set, and 0.13 mas for the OPA and IAA data sets. These are smaller than the WRMS obtained by Herring *et al.* (2002). They are also smaller than in Koot *et al.* (2008) because the Poisson terms are no longer part of the residuals. This is also why the parameter $\sigma_M$, characterizing model uncertainty, is smaller here than in Koot

*et al.* (2008). The direct consequence of a smaller $\sigma_M$ is smaller errors on the parameters.

## 5.3 Periods and Q of the normal modes

From the samples of the posterior pdf given by eq. (32), we can infer the pdf of the complex frequencies ($\sigma$) of the rotational normal modes that are in the frequency band of nutation, namely the FCN and FICN. These frequencies correspond to values of $\sigma$ for which Det[$M$]($\sigma$] = 0. The estimations for the corresponding periods T and quality factors Q, defined by: $\sigma + 1 = T^{-1}[1 - (i/2Q)]$, are given in Table 5 for the individual inversions.

# 6 COUPLING CONSTANTS AND PHYSICAL PROPERTIES OF THE CMB AND ICB

We now focus our attention on the coupling parameters. At the CMB, as Re($K_{CMB}$) cannot be separated from $\epsilon_f$, only Im($K_{CMB}$) can be used to infer information on the coupling. We are thus interested in the parameters Im($K_{CMB}$), Re($K_{ICB}$) and Im($K_{ICB}$). The numerical values of the latter two parameters are significantly different from previous nutation inversions. We explore below the reasons for these differences and investigate the consequence of these new values on the coupling at both the CMB and ICB.

## 6.1 Coupling constants

The estimates of Im($K_{CMB}$), Re($K_{ICB}$) and Im($K_{ICB}$) obtained from the OPA and IAA time-series are in close agreement (see Table 3).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Marginal posterior pdf of Re($K_{ICB}$) and Im($K_{ICB}$) for the individual inversion of the three time-series (GSFC, OPA and IAA), the joint inversion (Joint) and for the GSFC data truncated to the period 1990–2009.3 (GSFC1). The blue boxes are the 99.7 per cent CI obtained by Mathews *et al.* (2002) (MHB).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Residuals for the GSFC nutation data set in longitude (upper) and obliquity (lower).</figcaption></figure>

The estimates from the GSFC time-series are somewhat different but are still consistent with the two other time-series as the 99.7 per cent CI always have an intersection. The inversion of the GSFC1 time-series gives mean values slightly different from that of the complete GSFC time-series but again the 99.7 per cent CI widely overlap one another.

Our estimate of Im($K_{CMB}$) from the joint inversion is ($-1.78 \pm 0.02 \times 10^{-5}$). As mentioned above, and shown on Fig. 2, the pdf's of Im($K_{CMB}$) from the different inversions do not overlap one another exactly. The smallest (resp. largest) absolute value of Im($K_{CMB}$), within the 99.7 per cent CI, from all our inversions is Im($K_{CMB}$) = $-1.72 \times 10^{-5}$ from the IAA time-series [resp. Im($K_{CMB}$) = $-1.83 \times 10^{-5}$, from GSFC]. These upper and lower bounds on Im($K_{CMB}$) are perhaps more representative of the uncertainties than the very small errors obtained in each separate inversions. We use these values in our analysis below.

Similarly, our estimates of Re($K_{ICB}$) and Im($K_{ICB}$) from the joint inversion are (see Tab. 3) ($1.01 \pm 0.02 \times 10^{-3}$ and ($-1.09 \pm 0.03 \times 10^{-3}$, respectively. Upper and lower bounds within the

99.7 per cent CI from individual time-series are $1.06 \times 10^{-3}$ and $0.95 \times 10^{-3}$, and $-0.99 \times 10^{-3}$ and $-1.16 \times 10^{-3}$, respectively.

Although Re($K_{ICB}$) and Im($K_{ICB}$) are significantly different, our estimate of Im($K_{CMB}$) is in agreement with that of Mathews *et al.* (2002). This can be explained as follows. The parameter $K_{ICB}$ affects directly the complex frequency of the FICN (see eq. 16c). The period of this mode is approximately 1000 d in the celestial frame (Mathews *et al.* 2002, and Table 5) so that the forced nutations that

<table><caption>**Table 5.** Periods and quality factors of the FCN and FICN. The errors correspond to the 99.7 per cent CI.</caption><thead><tr><th rowspan="2">Data set</th><th colspan="2">FCN</th><th colspan="2">FICN</th></tr><tr><th>T (days)</th><th>Q</th><th>T (days)</th><th>Q</th></tr></thead><tbody><tr><td>GSFC</td><td>-429.09 ± 0.07</td><td>19641 ± 270</td><td>904 ± 29</td><td>467 ± 24</td></tr><tr><td>OPA</td><td>-429.00 ± 0.07</td><td>19716 ± 288</td><td>945 ± 30</td><td>455 ± 23</td></tr><tr><td>IAA</td><td>-429.05 ± 0.08</td><td>19886 ± 328</td><td>919 ± 33</td><td>453 ± 24</td></tr><tr><td>GSFC1</td><td>-429.07 ± 0.07</td><td>19721 ± 277</td><td>929 ± 31</td><td>480 ± 24</td></tr></tbody></table>
are the most affected by resonance are the ones with long periods and in particular the nutation of 18.6 yr. The amplitude of this nutation is known with less precision than the other dominant forced nutations because only 30 yr of data are available. The difference in $K_{ICB}$ comes then mainly from a better account of the long-period terms: all the long-period periodic and Poisson terms are part of our model (not only the 18.6 yr periodic term). A part of the difference in $K_{ICB}$ also originates simply from the 10 yr of additional nutation observations that are part of our inversion. In contrast, $\text{Im}(K_{CMB})$ is mainly determined from the amplitude of the retrograde annual nutation, the forced nutation with the closest frequency to that of the FCN. Since the amplitude of the retrograde annual nutation could already be well determined with data prior to 2000, our estimate of $\text{Im}(K_{CMB})$ is not different from that of Mathews et al. (2002).

## 6.2 Coupling model

Our new estimations of the coupling constants have implications for the coupling at the CMB and ICB and in turn, on the physical parameters on which these depend. Our goal here is not to develop new models of coupling. We are merely interested in determining the numerical values of the physical parameters that are consistent with our coupling constants. We restrict our attention to EM and viscous coupling, which have been the focus of the most recent studies. We adopt the viscomagnetic coupling model of Mathews & Guo (2005), in which the coupling constant at one boundary $b$ is given by

$$K_b = \frac{\pi a_b^4}{2i \Omega_0 A_s} \int_0^\pi (\mu_0^{-1} I_b^{\text{em}} + \rho \nu I_b^{\text{vis}}) \sin \theta \, d\theta, \quad (35)$$

where $a_b$ is the boundary radius, $\mu_0$ is the magnetic permeability of free space, and $\theta$ is colatitude. $I_b^{\text{em}}$ and $I_b^{\text{vis}}$ include the contributions from EM and viscous forces, respectively, to the boundary torque. They depend on the radial magnetic field at the boundary ($B_r$), on the electrical conductivities of the fluid core ($\sigma_f$) and the solid side of the boundary ($\sigma_m$ for the mantle and $\sigma_s$ for the inner core), on the kinematic viscosity of the fluid core close to the boundary ($\nu$), and on the density of the fluid core. The complete expressions for $I_b^{\text{em}}$ and $I_b^{\text{vis}}$, at the CMB and ICB, can be found in Mathews & Guo (2005). In the limit of zero viscosity, the second term of the integrand in (35) vanishes and the purely EM model of Buffett et al. (2002) is retrieved.

We fix the electrical conductivities of the outer and inner cores to that of iron under core conditions, namely $\sigma_f = \sigma_s = 5 \times 10^5 \text{ S m}^{-1}$ (Stacey & Anderson 2001). We take the densities of the fluid core at the CMB and ICB as given by PREM (Dziewonski & Anderson 1981), 9903.4 and 12166.3 kg m$^{-3}$, respectively.

The radial magnetic field at the CMB and ICB that enters $I_b^{\text{em}}$ is decomposed into a dipolar ($B_r^D$) and a uniform ($B_r^{\text{ND}}$) component, the latter representing an effective contribution to the EM torque from all field components other than the axial dipole. The viscous coupling part depends on the kinematic viscosity of the fluid core. Since it is poorly known, it is considered as an unknown parameter.

## 6.3 Coupling at the CMB

Our estimate of $\text{Im}(K_{CMB})$ is in agreement with that of Mathews et al. (2002) used in the study of Mathews & Guo (2005). Thus, we agree with the interpretation of $\text{Im}(K_{CMB})$ in terms of physical parameters presented in this latter study. Nevertheless, for completeness, we do include here an analysis of the coupling at the

CMB on the basis of updated values of $B_r^D$ and $B_r^{\text{ND}}$ based on the recent model CHAOS-2s (Olsen et al. 2009) constructed from satellite observations of the magnetic field. Also, in our study, we give the errors on the physical parameters that arise from the errors on $\text{Im}(K_{CMB})$, which was not done by Mathews & Guo (2005).

Assuming a perfectly insulating mantle, the long-wavelength components (spherical harmonic degrees smaller or equal to 13) of the magnetic field at the CMB can be determined from a downward continuation of the magnetic field observed at the Earth's surface. The amplitude of the axial dipole at the CMB is thus well determined though it is currently decreasing at a rate of approximately 150 nT per year (e.g Gubbins et al. 2006). Since this corresponds to a very small change over the time-span of the nutation observations, we simply assume a fixed rms value of $\tilde{B}_r^D = 0.209 \text{ mT}$ based on CHAOS-2s (Olsen et al. 2009) evaluated at J2000. As for rms strength of $B_r^{\text{ND}}$ at the CMB, though its long wavelength contribution is well determined, the short-wavelength contribution cannot be obtained from surface observations because it is masked by the crustal field. The true rms strength of $B_r^{\text{ND}}$ is thus unknown. However, the power spectrum of the field between degrees 2 and 13 is found to follow a log-linear trend. Assuming this trend continues to higher degrees, an estimate of $\tilde{B}_r^{\text{ND}}$ can be constructed by extrapolation. Based on the mean power spectrum of CHAOS-2s between 2000 and 2008, we get $\tilde{B}_r^{\text{ND}} = 0.281 \text{ mT}$. We refer to this value of $\tilde{B}_r^{\text{ND}}$ below as that inferred from magnetic field observations, though this remains an unknown quantity.

We consider first the case of a purely EM coupling. We use $\tilde{B}_r^D = 0.209 \text{ mT}$ and treat $\tilde{B}_r^{\text{ND}}$ as an unknown parameter. The other unknown parameter is $\sigma_m$. These two parameters cannot be both estimated from the single parameter $\text{Im}(K_{CMB})$. We choose to estimate $\tilde{B}_r^{\text{ND}}$ for different choices of $\sigma_m$. Following Buffett et al. (2002), we first consider the case of a very high lowermost mantle conductivity taken to be equal to that of the core, $\sigma_m = 5 \times 10^5 \text{ S m}^{-1}$.

Fig. 5 shows the real and imaginary parts of $K_{CMB}$ as a function of $\tilde{B}_r^{\text{ND}}$. The black solid horizontal line corresponds to $\text{Im}(K_{CMB}) = -1.78 \times 10^{-5}$, our estimate from the joint inversion (see Table 3). The dashed (resp. dotted) horizontal line corresponds to the upper (resp. lower) bound reported above: $\text{Im}(K_{CMB}) = -1.72 \times 10^{-5}$ [resp. $\text{Im}(K_{CMB}) = -1.83 \times 10^{-5}$]. The value $\text{Im}(K_{CMB}) = -1.78 \times 10^{-5}$ corresponds to $\tilde{B}_r^{\text{ND}} = 0.631 \text{ mT}$, while the two extreme values $\text{Im}(K_{CMB}) = -1.72 \times 10^{-5}$ and $-1.83 \times 10^{-5}$ give $\tilde{B}_r^{\text{ND}} = 0.617$ and 0.642 mT, respectively. These values, as well as the corresponding values of the rms of the total radial field $\tilde{B}_r^{\text{Tot}} = \sqrt{(\tilde{B}_r^D)^2 + (\tilde{B}_r^{\text{ND}})^2}$, are reported in Table 6.

Since $\sigma_m$ cannot be larger than the electrical conductivity of the fluid core, the retrieved value of $\tilde{B}_r^{\text{ND}} = 0.617 \text{ mT}$ under the above scenario corresponds to a lower bound on $\tilde{B}_r^{\text{ND}}$. Smaller values of $\sigma_m$ require larger values of $\tilde{B}_r^{\text{ND}}$: as reported in Table 6, with $\sigma_m = 5 \times 10^4 \text{ Sm}^{-1}$, $\tilde{B}_r^{\text{ND}}$ must be approximately 0.91 mT. These estimates of $\tilde{B}_r^{\text{ND}}$ are much larger than our above quoted value of 0.281 mT based on magnetic field observations. This has been interpreted to indicate that, if no other coupling than EM participates, the small wavelength component of the radial field at the CMB has much larger amplitude than that expected from an extrapolation of its long wavelength component (Buffett et al. 2002; Mathews et al. 2002; Buffett & Christensen 2007).

As an alternative to the large $\tilde{B}_r^{\text{ND}}$ required in the purely EM coupling scenario, a part of $\text{Im}(K_{CMB})$ may be explained by viscous coupling. In Fig. 5, we show how adding viscous coupling affects the real and imaginary parts of $K_{CMB}$ and how a particular value of $\text{Im}(K_{CMB})$ can be explained by different combinations of $\nu$ and $\tilde{B}_r^{\text{ND}}$.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Re($K_{CMB}$) and Im($K_{CMB}$) as a function of $\tilde{B}_r^{ND}$ at the CMB, for $\sigma_m = 5 \times 10^5$ S m$^{-1}$ and $\tilde{B}_r^D = 0.209$ mT. The results for three different choices of the fluid core viscosity are shown: $\nu = 0$ (orange), $\nu = 10^{-2}$ m$^2$ s$^{-1}$ (purple), and $\nu = 2.6 \times 10^{-2}$ m$^2$ s$^{-1}$ (blue). The horizontal lines are the values of Im($K_{CMB}$) estimated from nutation observations and the three black vertical lines are the corresponding values of $\tilde{B}_r^{ND}$ in the case of the purely EM model. The vertical red line correspond to $\tilde{B}_r^{ND} = 0.281$ mT inferred from magnetic field observation.</figcaption></figure>

<table><caption>**Table 6.** The rms of the radial magnetic field at the CMB required to match the estimated coupling constant for purely EM coupling. $\tilde{B}_r^D$ is fixed to 0.209 mT.</caption><thead><tr><th>σ<sub>m</sub><br>(S m<sup>-1</sup>)</th><th>Im(K<sub>CMB</sub>)<br>(10<sup>-5</sup>)</th><th>B&#773;<sub>r</sub><sup>ND</sup><br>(mT)</th><th>B&#773;<sub>r</sub><sup>Tot</sup><br>(mT)</th></tr></thead><tbody><tr><td rowspan="3">5 × 10<sup>5</sup></td><td>-1.78 ± 0.02</td><td>0.631 ± 0.004</td><td>0.664 ± 0.004</td></tr><tr><td>-1.72</td><td>0.617</td><td>0.651</td></tr><tr><td>-1.83</td><td>0.642</td><td>0.675</td></tr><tr><td rowspan="3">5 × 10<sup>4</sup></td><td>-1.78 ± 0.02</td><td>0.911 ± 0.006</td><td>0.935 ± 0.006</td></tr><tr><td>-1.72</td><td>0.893</td><td>0.917</td></tr><tr><td>-1.83</td><td>0.926</td><td>0.949</td></tr></tbody></table>

Fig. 6 shows the relation between $\nu$ and $\tilde{B}_r^{ND}$ in agreement with our Im($K_{CMB}$) estimates. Our curve for $\sigma_m = 5 \times 10^5$ S m$^{-1}$ is similar to that obtained by Mathews & Guo (2005); we also show the curves of $\nu$ for other choices of $\sigma_m$. As expected, the viscosity required to explain the estimated coupling constant increases as $\tilde{B}_r^{ND}$ decreases. In particular, assuming $\sigma_m = 5 \times 10^5$ S m$^{-1}$, in order to reconcile $\tilde{B}_r^{ND}$ with 0.281 mT (the lowest value shown on Fig. 6), the viscosity of the fluid core has to be between 0.023 and 0.028 m$^2$ s$^{-1}$. Smaller $\sigma_m$ require larger values of $\nu$, though the increase is modest: even for $\sigma_m$ reduced by four orders of magnitude, the required increase in $\nu$ is less than a factor 2. These values are upper bounds on the viscosity of the fluid core at the CMB and are reported in Table 7.

## 6.4 Coupling at the ICB

Our estimate of $K_{ICB}$ is not in agreement with that of Mathews et al. (2002). Both the real and imaginary parts are outside of the 99.7 per cent CI estimated by these authors. Moreover, for each

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Fluid core viscosity and radial uniform magnetic field at the CMB, for different values of the conductivity of the lowermost mantle $\sigma_m$. The solid lines correspond to Im($K_{CMB}$) obtained from the joint inversion, the dashed (resp. dotted) lines correspond to the lowest (resp. highest) absolute value of Im($K_{CMB}$) allowed by our various inversions.</figcaption></figure>

<table><caption>**Table 7.** The fluid core viscosity $\nu$ at the CMB required to match the estimated coupling constant when the magnetic field at the CMB is fixed to $\tilde{B}_r^D = 0.209$ mT and $\tilde{B}_r^{ND} = 0.281$ mT.</caption><thead><tr><th>σ<sub>m</sub><br>(S m<sup>-1</sup>)</th><th>Im(K<sub>CMB</sub>)<br>(10<sup>-5</sup>)</th><th>ν<br>(10<sup>-2</sup> m<sup>2</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="3">5 × 10<sup>5</sup></td><td>-1.78 ± 0.02</td><td>2.6 ± 0.1</td></tr><tr><td>-1.72</td><td>2.3</td></tr><tr><td>-1.83</td><td>2.8</td></tr><tr><td rowspan="3">5 × 10<sup>4</sup></td><td>-1.78 ± 0.02</td><td>3.3 ± 0.1</td></tr><tr><td>-1.72</td><td>3.0</td></tr><tr><td>-1.83</td><td>3.5</td></tr><tr><td rowspan="3">10</td><td>-1.78 ± 0.04</td><td>3.9 ± 0.1</td></tr><tr><td>-1.72</td><td>3.6</td></tr><tr><td>-1.83</td><td>4.1</td></tr></tbody></table>

of our inversions we find that Re($K_{ICB}$)/|Im($K_{ICB}$)| < 1, while in Mathews et al. (2002) found that Re($K_{ICB}$)/|Im($K_{ICB}$)| > 1. This ratio is important because it provides information on the nature of the coupling at the ICB. Fig. 7 shows, for a purely EM coupling model, how the ratio Re($K_{ICB}$)/|Im($K_{ICB}$)| varies as a function of $\tilde{B}_r^{Tot}$ at the ICB and for different partitions between $\tilde{B}_r^D$ and $\tilde{B}_r^{ND}$. Regardless of the value of $\tilde{B}_r^{Tot}$ or of the relative proportion between the dipolar and uniform component, the ratio Re($K_{ICB}$)/|Im($K_{ICB}$)| is always larger than 1. Consequently, our new value of $K_{ICB}$ appears incompatible with a purely EM coupling, at least on the basis of the model developed by Buffett et al. (2002).

It is important to note that there remains a degree of uncertainty in the ratio Re($K_{ICB}$)/|Im($K_{ICB}$)| that we infer. Since Re($K_{ICB}$) is determined by the resonance of the nutations with the FICN (see eq. 16c), larger values of Re($K_{ICB}$) would be obtained for larger value of the combination $\alpha_2 e_s + \nu^{el}$. The numerical values of these three parameters is given in Table 2. An increase of 5 per cent in $e_s$,
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The ratio Re($K_{ICB}$)/|Im($K_{ICB}$)| as a function of $\tilde{B}_r^{\text{Tot}}$ at the ICB. The different curves show different partitions of the field between $\tilde{B}_r^D$ and $\tilde{B}_r^{\text{ND}}$, from a purely uniform field (lowest curve) to a purely dipolar field (upper curve).</figcaption></figure>

corresponding for instance to a non-hydrostatic contribution to the flattening, would lead to an increase in Re($K_{ICB}$) of approximately $10^{-4}$ and to Re($K_{ICB}$)/|Im($K_{ICB}$)| ≈ 1. A larger increase in $e_s$, however, is incompatible with the upper bound of the non-hydrostatic contribution on $e_f$ (Mathews et al. 2002). Thus, though the exact ratio Re($K_{ICB}$)/|Im($K_{ICB}$)| is difficult to pin down, it seems unlikely that Re($K_{ICB}$)/|Im($K_{ICB}$)| can be much higher than about 1, and thus remains incompatible with the values suggested in Fig. 7.

For a purely viscous coupling, we expect Re($K_{ICB}$)/|Im($K_{ICB}$)| ≈ 0.1 (Busse 1968; Rochester 1976; Mathews & Guo 2005). It is thus possible to match our new estimates of $K_{ICB}$ if both viscous and EM forces contribute to the coupling at the ICB.

We use again the model of Mathews & Guo (2005) given by eq. (35), in which we prescribe $\sigma_f = \sigma_s = 5 \times 10^5 \text{ S m}^{-1}$. $\tilde{B}_r^D$ and $\tilde{B}_r^{\text{ND}}$ at the ICB cannot be determined from magnetic field observations. These two parameters, together with $\nu$ close to the ICB, are the three unknowns of our model. Only two parameters can be derived from nutation observations for constraining the coupling: Re($K_{ICB}$) and Im($K_{ICB}$). Because the coupling model has one more free parameter, we choose to fix the value of $\tilde{B}_r^D$ and we estimate the values of $\tilde{B}_r^{\text{ND}}$ and $\nu$ required to match Re($K_{ICB}$) and Im($K_{ICB}$).

We take $\tilde{B}_r^D$ to vary between 0 and the maximum value for which we can find a solution. For each chosen value of $\tilde{B}_r^D$, we find $\tilde{B}_r^{\text{ND}}$ and $\nu$ by inverting Re($K_{ICB}$) and Im($K_{ICB}$), using once more the Bayesian inversion method with MCMC sampling described in Section 4.2. The parameters that we want to estimate are $\mathbf{p} = \{\tilde{B}_r^{\text{ND}}, \nu\}$. The forward problem is the viscomagnetic coupling model $K_{\text{ICB}}^{\text{model}}(\mathbf{p})$ (eq. 35). The ‘data’ are our estimated values of Re($K_{ICB}$) and Im($K_{ICB}$). Because the marginal pdf obtained for these parameters are close to Gaussian distributions, we consider our estimate of those parameters as Gaussian ‘measurements’. The estimates of Re($K_{ICB}$) and Im($K_{ICB}$) are not correlated, so they are considered as independent measurements. The likelihood function

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Estimates of $\tilde{B}_r^{\text{ND}}$ and $\nu$ (solid lines) obtained from the inversion of the coupling constant at the ICB for different choices of $\tilde{B}_r^D$. The dashed lines correspond to the 99.7 per cent CI.</figcaption></figure>

(eq. 28), without model uncertainty, can be written as

$$L(\mathbf{d}|\mathbf{p}) \propto \frac{1}{\sigma_1 \sigma_2} \exp \left\{ -\frac{1}{2} \left( \frac{d_1 - \text{Re}(K_{\text{ICB}}^{\text{model}}(\mathbf{p}))}{\sigma_1} \right)^2 - \frac{1}{2} \left( \frac{d_2 - \text{Im}(K_{\text{ICB}}^{\text{model}}(\mathbf{p}))}{\sigma_2} \right)^2 \right\}, \quad (36)$$

where $\mathbf{d}$ is the ‘data’, that is, $d_1$ is Re($K_{ICB}$) and $d_2$ is Im($K_{ICB}$) taken from Table 3, and $\sigma_1$, $\sigma_2$ are the corresponding standard deviations.

The prior pdf for $\nu$ is chosen to be uniform on the interval $0–100 \text{ m}^2 \text{ s}^{-1}$ and zero elsewhere. We assume $\tilde{B}_r^{\text{ND}}$ has a Jeffreys prior [i.e. a prior pdf in $1/x$, see e.g. Gregory (2005)] on the interval $0–20 \text{ mT}$ and zero elsewhere.

From the posterior pdf, we compute the MAP solutions and the 99.7 per cent CI. Fig. 8 and Table 8 show the results obtained when $K_{ICB}$ is taken from the joint inversion of the three time-series, namely for Re($K_{ICB}$) = (1.01 ± 0.02) 10⁻³ and Im($K_{ICB}$) = (-1.09 ± 0.03) 10⁻³. Results for $K_{ICB}$ obtained from the individual inversions are very similar and are not shown here.

In Fig. 8 and Table 8, we also show the corresponding values and CI of $\tilde{B}_r^{\text{Tot}}$. Samples of the pdf for $\tilde{B}_r^{\text{Tot}}$ are obtained from the samples of the marginal of $\tilde{B}_r^{\text{ND}}$ by computing, for each sample, $\tilde{B}_r^{\text{Tot}} = \sqrt{(\tilde{B}_r^{\text{ND}})^2 + (\tilde{B}_r^D)^2}$. From these samples, we compute the MAP values and 99.7 per cent CI.

As it can be seen from Table 8 and Fig. 8, $\tilde{B}_r^D$ varies from 0 to 6.1 mT and the corresponding value of $\tilde{B}_r^{\text{ND}}$ varies from 6.5 mT to 0. The value of $\tilde{B}_r^{\text{ND}}$ decreases as $\tilde{B}_r^D$ increases. The rms of the total radial magnetic field remains almost constant and approximately equal to 6.5 mT, regardless of the individual values of $\tilde{B}_r^D$ and $\tilde{B}_r^{\text{ND}}$. The viscosity of the fluid core at the ICB is almost constant for values of $\tilde{B}_r^D$ between 0 and 3 mT. This value, around 10 m² s⁻¹, is a lower bound on the value of this physical parameter. For $\tilde{B}_r^D$ larger than 3 mT, the viscosity begins to increase and reaches a maximum around 30 m² s⁻¹, corresponding to the maximum value of $\tilde{B}_r^D$, namely 6.1 mT, and the minimum value of $\tilde{B}_r^{\text{ND}}$, namely 0 mT.
<table><caption>**Table 8.** Estimates of $v$, $\tilde{B}_r^{ND}$ and $\tilde{B}_r^{Tot}$ for specific values of $\tilde{B}_r^D$.</caption><thead><tr><th>$\tilde{B}_r^D$</th><th>$v$</th><th>$\tilde{B}_r^{ND}$</th><th>$\tilde{B}_r^{Tot}$</th></tr></thead><tbody><tr><td>0</td><td>9.97<br>(8.39, 11.54)</td><td>6.51<br>(6.4, 6.63)</td><td>6.51<br>(6.4, 6.63)</td></tr><tr><td>1</td><td>9.88<br>(8.41, 11.36)</td><td>6.49<br>(6.38, 6.6)</td><td>6.57<br>(6.46, 6.67)</td></tr><tr><td>2</td><td>9.72<br>(8.23, 11.22)</td><td>6.39<br>(6.29, 6.5)</td><td>6.7<br>(6.6, 6.8)</td></tr><tr><td>3</td><td>9.86<br>(8.37, 11.36)</td><td>6.16<br>(6.04, 6.28)</td><td>6.85<br>(6.74, 6.96)</td></tr><tr><td>4</td><td>12.09<br>(10.55, 13.63)</td><td>5.58<br>(5.43, 5.73)</td><td>6.87<br>(6.74, 6.99)</td></tr><tr><td>5</td><td>18.46<br>(16.85, 20.06)</td><td>4.35<br>(4.13, 4.58)</td><td>6.63<br>(6.48, 6.78)</td></tr><tr><td>5.4</td><td>22.04<br>(20.48, 23.61)</td><td>3.55<br>(3.27, 3.83)</td><td>6.47<br>(6.31, 6.62)</td></tr><tr><td>6</td><td>28.54<br>(26.96, 30.12)</td><td>1.23<br>(0, 1.9)</td><td>6.11<br>(6., 6.3)</td></tr><tr><td>6.1</td><td>29.52<br>(28.00, 31.04)</td><td>0<br>(0, 1.1)</td><td>6.1<br>(6.1, 6.2)</td></tr></tbody></table>

*Notes:* The first line gives the estimated value whereas the values in parentheses on the second line correspond to the 99.7 per cent CI. The results are given for $K_{ICB}$ based on the joint inversion of GSFC, OPA and IAA time-series. Magnetic field values are in mT and kinematic viscosities are in $m^2 s^{-1}$.

We note finally that it is not possible to explain our estimated value of $K_{ICB}$ by viscous coupling alone. This is because the ratio Re($K_{ICB}$)/|Im($K_{ICB}$)| does not permit it, regardless of how large the viscosity near the ICB can be. Based on our results, $30 m^2 s^{-1}$ corresponds to an approximate upper bound for the kinematic viscosity near the ICB. Values as high as $10^7 m^2 s^{-1}$ that have been proposed (e.g. Smylie 1999) are incompatible with nutation observations.

## 7 DISCUSSION AND CONCLUSION

The coupling constant at the CMB that we have inferred from nutation observations, is in agreement with the results previously obtained by Mathews et al. (2002). We find $Im(K_{CMB}) = (-1.78 \pm 0.02) 10^{-5}$ for an inversion that includes three different nutation time-series. The errors correspond to 99.7 per cent confidence interval, though these are likely underestimates of the true uncertainties in $Im(K_{CMB})$. The highest and lowest value of $Im(K_{CMB})$ within the 99.7 per cent CI that we find from inversions of individual time-series are $-1.72 \times 10^{-5}$ and $-1.83 \times 10^{-5}$, respectively.

The physical properties of the CMB that we derive on the basis of $Im(K_{CMB})$ are thus in agreement with those of Buffett et al. (2002) and Mathews & Guo (2005). If the coupling is purely of EM nature, the non-dipolar part of the rms radial field at the CMB, $\tilde{B}_r^{ND}$, must be at least $0.631 \pm 0.004$ mT or larger in order to explain the strength of the coupling. This is much larger than $\tilde{B}_r^{ND} = 0.281$ mT inferred on the basis of an extrapolation of the long wavelength power spectrum of the field model CHAOS-2s. However, it is not inconceivable that the small wavelength part of the field at the CMB contains more energy than based on such a simple extrapolation (Buffett & Christensen 2007). If this is the case, then $\tilde{B}_r^{ND} = 0.631$ mT may be taken as a constraint from nutation observations.

If instead one takes $\tilde{B}_r^{ND} = 0.281$ mT to represent an upper bound, then an additional coupling mechanism must be present. If this coupling is from viscous forces, the required kinematic viscosity of the fluid near the CMB is between $0.023$ and $0.041 m^2 s^{-1}$, depending on the electrical conductivity of the lower mantle and

the precise value of $Im(K_{CMB})$. These values are upper bounds on the viscosity of the fluid core at the CMB.

For the real part of the coupling constant at the ICB, we find $Re(K_{ICB}) = (1.01 \pm 0.02) 10^{-3}$ for the inversion that includes all three nutation time-series. Upper and lower bounds within the 99.7 per cent CI from individual time-series are, $1.06 \times 10^{-3}$ and $0.95 \times 10^{-3}$, respectively. These values are lower, but remain compatible within error bars with the value $Re(K_{ICB}) = (1.11 \pm 0.11) 10^{-3}$ obtained by Mathews et al. (2002). However, the imaginary part of $K_{ICB}$ that we find is significantly different: $Im(K_{ICB}) = (-1.09 \pm 0.03) 10^{-3}$, for the joint inversion, with upper and lower bounds of $-0.99 \times 10^{-3}$ and $-1.16 \times 10^{-3}$, compared to $Im(K_{ICB}) = (-0.78 \pm 0.14) 10^{-3}$ obtained by Mathews et al. (2002).

Our new values of $Re(K_{ICB})$ and $Im(K_{ICB})$ have profound implications for the nature of the coupling at the ICB. The previous estimate of these parameters by Mathews et al. (2002) could be explained by EM coupling alone: this required a total rms radial field of 7.17 mT at the ICB. This is a much larger value than the typical magnetic field amplitudes in the core of 2–3 mT expected on the basis of numerical simulations of the geodynamo (Christensen & Aubert 2006), but could be taken as an observational constraint from nutations. However, with our new values of $Re(K_{ICB})$ and $Im(K_{ICB})$, the coupling at the ICB cannot be explained in terms of EM coupling alone, at least on the basis of existing models, regardless of the amplitude of the field. This is because the ratio $Re(K_{ICB})/|Im(K_{ICB})|$ that we find is smaller than 1, whereas the EM coupling model can only accommodate a ratio larger than 1 (see Fig. 7).

An additional coupling mechanism must then be present. If this coupling is of viscous nature, kinematic viscosities between 10 and $30 m^2 s^{-1}$ are required, depending on the partition between the dipolar and non-dipolar part of the radial magnetic field at the ICB. These correspond to approximate upper bounds on the kinematic viscosity near the ICB. Higher values are incompatible with nutation observations. Importantly, we note that viscous coupling alone cannot explain our values of $Re(K_{ICB})$ and $Im(K_{ICB})$: a part of the coupling must still be accomplished by EM forces. Our estimate of the rms of the total radial field is between 6 and 7 mT, a value that remains much larger than that inferred from geodynamo simulations.

The kinematic viscosities quoted above, at both the CMB and ICB, are very large. Typical values inferred from laboratory measurements (Rutter et al. 2002b; Rutter et al. 2002a) and ‘ab initio’ computations (Alfè et al. 2000) suggest that the molecular viscosity of iron at temperatures and pressures corresponding to that of the core should be of the order of $10^{-6} m^2 s^{-1}$. Our much larger viscosity estimates must then be taken as representing an effective viscosity $v_e$ caused by turbulent motion transport in the fluid (e.g. Deleplace & Cardin 2006), similar to the eddy viscosity concept adopted when modelling boundary layers in the ocean and atmosphere (Pedlosky 1987).

Such a concept may seem appropriate at first glance. Adopting a simple isotropic turbulence framework, the effective kinematic viscosity should be approximately equal to other diffusivities in the system. The largest molecular diffusivity in the Earth’s core is the magnetic diffusivity $\eta = 1/\mu_0 \sigma_f \approx 2 m^2 s^{-1}$, a value approximately compatible with the effective viscosities that we infer. However, turbulence in the Earth’s core is expected to be highly anisotropic on account of Coriolis and Lorentz forces (Braginsky & Meytlis 1990). It is then not clear whether taking $v_e \approx \eta$ is justified.

The eddy viscosity interpretation also has additional shortcomings. An effective viscosity is a property of the flow, not a property
of the fluid. The size and turn-over timescale of the eddies that act as an apparent viscosity must be compatible with those expected to be present in the flow. A discussion on these issues is presented in Buffett & Christensen (2007). We shall not repeat this discussion here, but it is worth emphasizing a few important points. The eddy viscosity should scale as $v_e \approx vl$, where $v$ and $l$ are typical velocity and eddy size ($l$ is sometimes referred to as a mixing length). An upper bound limit on $l$ is the thickness $\delta$ of the effective viscous boundary layer, otherwise the whole concept of eddy transport within the boundary layer falls appart. $\delta$ should scale as $(v_e/\Omega_o)^{1/2}$ based on an Ekman boundary layer thickness (e.g. Greenspan 1968). This leads to $l \sim \delta \sim 300$ m when $v_e = 10 \text{ m}^2 \text{ s}^{-1}$, and requiring $v$ to be of the order of 0.03 m s$^{-1}$. This is a velocity two orders of magnitude larger than the typical velocities of the large scale flows in the core. Realistically, $l$ should be much smaller than $\delta \sim 300$ m for eddies to act like an apparent viscosity, in which case even higher eddy velocities are required. The same exercise with a viscosity typical of those that we retrieve near the CMB, $v_e = 0.03 \text{ m}^2 \text{ s}^{-1}$, leads to $l \sim \delta \sim 20$ m and eddy velocities of at least 0.0015 m s$^{-1}$, still an order of magnitude larger than typical core flows.

Whether typical eddy velocities higher than 0.03 m s⁻¹ and mixing lengths smaller than ~100 m are compatible with the energetics and the force balance in the fluid core is clearly too vast a question to be tackled here. Though these would be the typical values that are required to explain the coupling at boundaries inferred from nutation observations.

We thus reach the conclusion that, at the ICB, not only a large radial magnetic field is required to explain the coupling constant, a large fluid viscosity difficult to reconcile with molecular or eddy viscosity values is also required. As these large values for both the magnetic field and fluid viscosity appear to be incompatible with other geophysical constraints, this may indicate that the present-day models of EM and viscous coupling have serious inadequacies. For instance, the assumption that the differential rotation of the fluid core in the characterization of the free modes is adequately represented by a rigid motion may be incorrect. Using a Largrangian approach, Rogister & Valette (2009) has shown that indeed, for both the FCN and FICN, departures from rigid motion may be important. These non-rigid motions need to be taken into account in the coupling between the flow and the boundaries. Alternatively, other unmod- elled physical processes may also be at play. For instance, viscous dissipation in the inner core may be important on diurnal timescales and could explain a significant part of $K_{\text{ICB}}$ (e.g. Greff-Lefftz et al. 2000). Some of these ideas are currently under investigation.

## ACKNOWLEDGMENTS

Part of the work of L. Koot was performed under a research fel- lowship of the Belgian National Fund for Scientific Research. L. Koot and A. Rivoldini acknowledge the support of the Belgian Sci- ence Policy Office (Action 1 Projects). M. Dumberry is currently supported by a NSERC/CRSNG Discovery grant. The contribution of O. de Viron to this paper is IPGP contribution 3008. We thank T. A. Herring, the editor B. A. Buffett, and an anonymous referee for their useful comments that helped us to improve the paper.

## REFERENCES

Alfè, D., Kresse, G. & Gillan, M.J., 2000. Structure and dynamics of liquid iron under Earth's core conditions, *Phys. Rev. B*, **61**, 132–142.

Bizouard, C., Brzeziński, A. & Petrov, S., 1998. Diurnal atmospheric forcing and temporal variations of the nutation amplitudes, *J. Geod.*, **72**, 561–577.

Braginsky, S.I. & Meytlis, V.P., 1990. Local turbulence in the Earth's core, *Geophys. Astrophys. Fluid Dyn.*, **55**, 71–87.

Bretagnon, P., Francou, G., Rocher, P. & Simon, J.L., 1998. SMART97: a new solution for the rotation of the rigid Earth, *Astron. Astrophys.*, **329**, 329–338.

Buffett, B.A. & Christensen, U.R., 2007. Magnetic and viscous coupling at the core-mantle boundary: inferences from observations of the Earth's nutations, *Geophys. J. Int.*, **171**, 145–152.

Buffett, B.A., Mathews, P.M. & Herring, T.A., 2002. Modeling of nutation and precession: effects of electromagnetic coupling, *J. geophys. Res.*, **107**, 2070, doi:10.1029/2000JB000056.

Busse, F.H., 1968. Steady fluid flow in a precessing spheroidal shell, *J. Fluid Mech.*, **33**, 739–751.

Chao, B.F., Ray, R.D., Gipson, J.M., Egbert, G.D. & Ma, C., 1996. Diurnal/semidiurnal polar motion excited by oceanic tidal angular momentum, *J. geophys. Res.*, **101**, 20151–20164.

Christensen, U.R. & Aubert, J., 2006. Scaling properties of convection- driven dynamos in rotating spherical shells and application to planetary magnetic fields, *Geophys. J. Int.*, **166**, 97–114.

Defraigne, P., Dehant, V. & Wahr, J.M., 1996. Internal loading of an inhomogeneous compressible earth with phase boundaries, *Geophys. J. Int.*, **125**, 173–192.

Deleplace, B. & Cardin, P., 2006. Viscomagnetic torque at the core mantle boundary, *Geophys. J. Int.*, **167**, 557–566.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth model, *Phys. Earth planet. Int.*, **25**, 297–356.

Feissel-Vernier, M., Ma, C., Gontier, A. & Barache, C., 2006. Analysis strategy issues for the maintenance of the ICRF axes, *Astron. Astrophys.*, **452**, 1107–1112.

Folgueira, M., Dehant, V., Lambert, S.B. & Rambaux, N., 2007. Impact of tidal Poisson terms on nonrigid Earth rotation, *Astron. Astrophys.*, **469**, 1197–1202.

Fukushima, T., 1991. Geodesic nutation, *Astron. Astrophys.*, **244**, L11+.

Greenspan, H.P., 1968. *The Theory of Rotating Fluids*, Cambridge University Press, New York.

Greff-Lefftz, M., Legros, H. & Dehant, V., 2000. Influence of the inner core viscosity on the rotational eigenmodes of the Earth, *Phys. Earth planet. Int.*, **122**, 187–204.

Gregory, P., 2005. *Bayesian Logical Data Analysis for the Physical Sciences*, Cambridge University Press, Cambridge.

Gubbins, D., Jones, A.L. & Finlay, C.C., 2006. Fall in Earth's magnetic field is erratic, *Science*, **312**, 900–902.

Hastings, W.K., 1970. Monte carlo sampling methods using markov chains and their applications., *Biometrika*, **57**, 97–109.

Herring, T.A., Mathews, P.M. & Buffett, B.A., 2002. Modeling of nutation- precession: very long baseline interferometry results, *J. geophys. Res.*, **107**, 2069, doi:10.1029/2001JB000165.

Koot, L., Rivoldini, A., de Viron, O. & Dehant, V., 2008. Estimation of Earth interior parameters from a Bayesian inversion of very long baseline interferometry nutation time series, *J. geophys. Res.*, **113**, 8414, doi:10.1029/2007JB005409.

Lambert, S.B., 2006. Atmospheric excitation of the Earth's free core nutation, *Astron. Astrophys.*, **457**, 717–720.

Lambert, S.B. & Mathews, P.M., 2006. Second-order torque on the tidal redistribution and the Earth's rotation, *Astron. Astrophys.*, **453**, 363–369.

Mathews, P.M. & Bretagnon, P., 2003. Polar motions equivalent to high frequency nutations for a nonrigid Earth with anelastic mantle, *Astron. Astrophys.*, **400**, 1113–1128.

Mathews, P.M. & Guo, J.Y., 2005. Viscoelectromagnetic coupling in precession-nutation theory, *J. geophys. Res.*, **110**, 2402, doi:10.1029/2003JB002915.

Mathews, P.M., Buffett, B.A., Herring, T.A. & Shapiro, I.I., 1991a. Forced nutations of the earth: influence of inner core dynamics. 1. Theory, *J. geophys. Res.*, **96**, 8219–8242.
Mathews, P.M., Buffett, B.A., Herring, T.A. & Shapiro, I.I., 1991b. Forced nutations of the Earth: influence of inner core dynamics. 2. Numerical results and comparisons, *J. geophys. Res.*, **96**, 8243–8258.

Mathews, P.M., Herring, T.A. & Buffett, B.A., 2002. Modeling of nutation and precession: new nutation series for nonrigid Earth and insights into the Earth's interior, *J. geophys. Res.*, **107**, 2068, doi:10.1029/2001JB000390.

Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H. & Teller, E., 1953. Equations of state calculations by fast computing machines, *J. Chem. Phys.*, **21**, 1087–1091.

Moritz, H. & Mueller, I.I., 1987. *Earth Rotation: Theory and Observation*, Ungar Publishing Company, New York.

Olsen, N., Mandea, M., Sabaka, T.J. & Tøffner-Clausen, L., 2009. CHAOS-2-a geomagnetic field model derived from one decade of continuous satellite data, *Geophys. J. Int.*, **179**, 1477–1487.

Pedlosky, J., 1987. *Geophysical Fluid Dynamics*, 2nd edn, Springer-Verlag, New York.

Rochester, M.G., 1976. The secular decrease of obliquity due to dissipative core-mantle coupling, *Geophys. J. R. astr. Soc.*, **46**, 109–126.

Rogister, Y. & Valette, B., 2009. Influence of liquid core dynamics on rotational modes, *Geophys. J. Int.*, **176**, 368–388.

Roosbeek, F. & Dehant, V., 1998. RDAN97: an analytical development of rigid earth nutation series using the torque approach, *Celest. Mech. Dyn. Astron.*, **70**, 215–253.

Rutter, M.D., Secco, R.A., Liu, H., Uchida, T., Rivers, M.L., Sutton, S.R. & Wang, Y., 2002a. Viscosity of liquid Fe at high pressure, *Phys. Rev. B*, **66**(6), 060102, doi:10.1029/2001GL014392.

Rutter, M.D., Secco, R.A., Uchida, T., Liu, H., Wang, Y., Rivers, M.L. &

Sutton, S.R., 2002b. Towards evaluating the viscosity of the Earth's outer core: an experimental high pressure study of liquid Fe-S (8.5 wt. per cent S), *Geophys. Res. Lett.*, **29**(8), 080000–1.

Sasao, T. & Wahr, J.M., 1981. An excitation mechanism for the free 'core nutation', *Geophys. J. R. astr. Soc.*, **64**, 729–746.

Sasao, T., Okubo, S. & Saito, M., 1980. A simple theory on dynamical effects of stratified fluid core upon nutational motion of the earth, in *Nutation and the Earth's Rotation*, IAU Symp. Vol. 78, pp. 165, eds Federov, E.P., Smith, M.L. & Bender, P.L., D. Reidel Publishing Co., Norwell, MA.

Smylie, D.E., 1999. Viscosity near Earth's solid inner core, *Science*, **284**, 461–463.

Souchay, J., Loysel, B., Kinoshita, H. & Folgueira, M., 1999. Corrections and new developments in rigid earth nutation theory. III. Final tables "REN-2000" including crossed-nutation and spin-orbit coupling effects, *Astron. Astrophys. Suppl.*, **135**, 111–131.

Stacey, F.D. & Anderson, O.L., 2001. Electrical and thermal conductivities of Fe-Ni-Si alloy under core conditions, *Phys. Earth planet. Int.*, **124**, 153–162.

Stacey, F.D. & Davis, P.M., 2008. *Physics of the Earth*, Cambridge University Press, Cambridge, ISBN: 9780521873628.

Tarantola, A., 2005. *Inverse Problem Theory and Methods for Model Parameter Estimation*, Society for Industrial and Applied Mathematics.

Wahr, J.M. & Sasao, T., 1981. A diurnal resonance in the ocean tide and in the earth's load response due to the resonant free 'core nutation', *Geophys. J. R. astr. Soc.*, **64**, 747–765.

Yseboodt, M., de Viron, O., Chin, T.M. & Dehant, V., 2002. Atmospheric excitation of the Earth's nutation: comparison of different atmospheric models, *J. geophys. Res.*, **107**, 2036, doi:10.1029/2000JB000042.