# Length of day variations due to mantle dynamics at geological timescale

Marianne Greff-Lefftz

► To cite this version:

Marianne Greff-Lefftz. Length of day variations due to mantle dynamics at geological timescale. Geophysical Journal International, 2011, 187, pp.595-612. 10.1111/j.1365-246X.2011.05169.x. insu-03606431

HAL Id: insu-03606431

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

Submitted on 12 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Length of day variations due to mantle dynamics at geological timescale

Marianne Greff-Lefftz

Institut de Physique du Globe de Paris 1 rue Jussieu, 75238 Paris Cedex 05, France. E-mail: greff@ipgp.fr

Accepted 2011 July 25. Received 2011 July 21; in original form 2011 March 31

## SUMMARY

The geological evolution of length of day (LOD) variations is mainly controlled by the frictional tidal torque responsible for the secular slowdown of the Earth's rotation and for the receding of the Moon. Superimposed on this variation, which has existed since the early history of the planet, there are, at shorter timescales (less than 1 Myr), LOD perturbations induced by the glaciation–deglaciation cycles. In this paper, we investigate the influence of mantle dynamics on LOD at the geological timescale. We use the complete non-linear equations to compute the influence of mantle density heterogeneities on the angular velocity of the rotation, that is to say on the LOD. We discuss the degree zero coefficient of the spherical harmonic expansion of the mantle mass anomaly, which is strongly dependent on the conservation of mass of the Earth.

We first compute the effects induced by upwelling domes and subducted plates sinking into the mantle, which are known to induce geological variations in the orientation of the rotation axis with respect to a fixed terrestrial frame (the so-called True Polar Wander). We find that the time-variable mantle density heterogeneities associated with the large-scale pattern of plate tectonic motions since 120 Ma and with the upwelling domes can perturb LOD by about 0.4 µs per year, that is, with an order of magnitude smaller than the effects induced by the last deglaciation. Superimposed on this linear trend, we show that there are fluctuations around 0.1–0.2 s Myr⁻¹ on a timescale of a few tens of millions of years.

We then combine a spherical model of mantle circulation with solutions for the equations governing the rotation of a viscous planet, to improve the constraint of mantle mass conservation.

Finally, we compare our results with other effects.

**Key words:** Earth rotation variations; Dynamics: convection currents, and mantle plumes; Planetary interiors.

## 1 INTRODUCTION

The length of day (LOD) varies under the action of external torques and also when mass moves on the surface or inside the Earth, affecting its angular momentum on different timescales.

(i) On a seasonal scale, both changes in the atmosphere and in ocean currents induce fluctuations of the LOD (e.g. Chao 1994; Chen 2005). On this short timescale, other factors affect the LOD: large-scale groundwater flux and the shape of the ocean surface in response to the movement of air masses above it.

(ii) On the decade scale, a slow drift seems to be related to motions within the Earth's fluid core (Jault *et al.* 1988; Jault & Le Mouël 1991), and more recently Dumberry & Bloxham (2006), and on a secular scale, the glaciation–deglaciation cycles involve secular acceleration or deceleration of the LOD (Nakiboglu & Lambeck 1980; Sabadini & Peltier 1981; Wahr *et al.* 1993).

(iii) On the geological timescale, the braking torque due to lunar tidal forces acting on the Earth has a long-term effect of slowing the Earth down (by about 2 ms per century) and speeding up the orbital rate of the Moon (which therefore recedes from Earth by a few centimetres per year, Lambeck 1980). Superimposed on this variation, which has existed since the early history of the planet, there are possible variations induced by dynamic redistribution of mantle mass.

The main features of mantle dynamics which can perturb the rotation (i.e. at very large length-scales) are downwelling plates and upwelling large plumes or superswells which are observed in tomography studies (e.g. Forte & Mitrovica 2001; Romanowicz & Gung 2002). In this paper, we compute the influence of upwelling domes and sinking plates on LOD, using very simple models of mantle dynamics.
The paper is organized as follows. In the first part, we present the non-linear equations governing the rotational dynamics of a viscous planet, taking into account simultaneously the LOD variation and rotational polar wander. The mantle density heterogeneities due to upwelling domes and to sinking plates, as well as the associated viscoelasto-gravitational deformations are described in Section 3. In Section 4, we compute the LOD changes induced by two stationary equatorial domes superimposed on realistic mantle density anomalies associated with the plates which have been subducted in the mantle since 120 Ma as derived by Ricard et al. (1993b); the conservation of mass and the resulting heating of the mantle are discussed. A mantle circulation model derived by Monnereau & Quere (2001) is also used to improve our results. Finally, we compare our results with the observable variations of LOD generated by other effects and we conclude.

## 2 THEORY OF THE EARTH'S ROTATION ON GEOLOGICAL TIMESCALES

To investigate if there is perturbation of LOD induced by degree 2 mantle mass anomaly, it is necessary to consider the non-linear Liouville equations. As a matter of fact, at geological timescale, the axis of rotation moves with respect to the reference frame of the mantle and this motion has a large amplitude (more than a few tens of degrees) and cannot be linearized. Taking into account simultaneously the degree 0 and 2 of the mantle density heterogeneities allows us to point out if there could be fluctuations of LOD on a timescale ranging from few tens of millions of years up to few hundred of millions of years.

We first recall the equations governing the rotational dynamics of a viscous planet. The temporal evolution of the rotation axis of the Earth is governed by the theorem of angular momentum $\vec{H}$. In a reference frame corotating with the mantle,

$$ \frac{d\vec{H}}{dt} + \vec{\omega} \wedge \vec{H} = \vec{L}, \qquad (1) $$

where $\vec{L}$ is an external torque and $\vec{H} = J \cdot \vec{\omega}$. J is the inertia tensor of the Earth and $\vec{\omega}$ is the angular velocity of the reference frame. The time-dependent inertia tensor $J(t)$ may be written as the sum of three parts;

$$ J_{ij} = I_o \delta_{ij} + C_{ij} + I_{ij}, \qquad (2) $$

where $\delta_{ij}$ the Kronecker symbol.

The first term $I_o$ is the inertia tensor of a spherical non-rotating Earth ($I_o = 0.33Ma^2$ for the actual Earth with a mass M and a radius a). The second term $C_{ij}$ is the inertia perturbation associated with the response of the Earth to the changing centrifugal potential: the difference between the associated principal moments of inertia represents the rotational ‘bulge’. The last term $I_{ij}$ describes changes in inertia due to internal mass redistribution.

$C_{ij}$ is equal to the time convolution of the degree 2 tidal Love number $k^T(t)$ with the time history of the changes in the centrifugal potential (e.g. Lefftz et al. 1991; Ricard et al. 1993a)

$$ C_{ij}(t) = \frac{a^5}{3G} k^T(t) * \left[ \omega_i(t)\omega_j(t) - \frac{1}{3}\omega^2(t)\delta_{ij} \right], \qquad (3) $$

where G is the gravitational constant and * denotes time convolution.

For a viscoelastic body, the tidal Love number $k^T$ may be written in the Laplace domain s and in the time domain as (e.g. Peltier 1974; Spada et al. 1990)

$$ k^T(s) = k^e + \sum_{j=1}^{NM} \frac{r_j}{s + 1/\tau_j} \qquad k^T(t) = k^e \delta(t) + \sum_{j=1}^{NM} r_j e^{-t/\tau_j} H(t), \qquad (4) $$

where $k^e$ is the tidal elastic Love number, $\tau_j$ are the NM viscoelastic relaxation times associated with the earth model, $r_j$ the viscous tidal coefficients, and $H(t)$ is the Heaviside step function.

For the long timescales associated with mantle convection, $t/\tau_j \gg 1$. Consequently, in the Laplace domain, we have $s\tau_j \ll 1$ and we neglect, in the Love numbers, the terms in $(s\tau_j)^2$. This approximation is equivalent to neglecting second time-derivative terms in the time domain. The Love number $k^T$ is then approximately equal to the quasi-fluid Love number (Munk & MacDonald 1960; Ricard et al. 1993a)

$$ k^T(s) \approx k^e + \sum_{j=1}^{NM} r_j \tau_j (1 - s\tau_j) \qquad k^T(s) \approx k^f [1 - T_1 s] \quad \text{with} \quad T_1 = \frac{1}{k^f} \sum_{j=1}^{NM} r_j \tau_j^2, \qquad (5) $$

where $k^f = k^e + \sum_{j=1}^{NM} r_j \tau_j$ is the tidal fluid Love number. In the temporal domain

$$ k^T(t) \approx k^f [\delta(t) - T_1 \dot{\delta}(t)], \qquad (6) $$

where the dot denotes the time-derivative. $T_1$ is the characteristic time of the delayed readjustment of the equatorial bulge of a rotating planet toward the fluid position (i.e. perpendicular to the rotation axis). With this quasi-fluid approximation, the inertia perturbation due to the
centrifugal potential becomes (Ricard et al. 1993a)

$$C_{ij}(t) = \frac{C_o}{\Omega^2} \left[ \omega_i \omega_j - \frac{\omega^2}{3} \delta_{ij} - T_1 \left( \dot{\omega}_i \omega_j + \omega_i \dot{\omega}_j - \frac{2}{3} \dot{\omega}_k \omega_k \delta_{ij} \right) \right] \quad (7)$$

with $C_o = \frac{\omega_i(t)}{\Omega^5} \Omega^2 k^f$, $\Omega$ being the sidereal rotation rate.

We introduce non-dimensional variables $m_i(t) = \frac{\omega_i(t)}{\Omega}$ and we define $m = \sqrt{m_1^2 + m_2^2 + m_3^2}$.

On geological timescales, $\frac{d\tilde{H}}{dt} \ll \tilde{\omega} \wedge \tilde{H}$: the ‘geological’ time is large in comparison with the ‘diurnal’ time related to the sidereal rotation. From the conservation of the kinetic energy of the Earth, we can show that neglecting $\frac{d\tilde{H}}{dt}$ forbids any temporal variation of the amplitude of the rotation vector, implying $\omega^2 = \Omega^2$: the LOD cannot vary in this approach. This hypothesis is usually made in studies relative to the wander of the rotational axis. Here, we are interested in the LOD variation, and consequently we have to keep this term in the conservation of angular momentum.

To be consistent with a quasi-fluid approximation of the tidal Love number, we have to neglect the second time derivative [i.e. terms in $\frac{d^2m_1}{dt^2}$, $\frac{dm_1}{dt}$, $\frac{d^2m_2}{dt^2}$, ... in eq. (1)].

The equations governing the rotational dynamics of a viscous planet may be written, from the conservation of angular momentum

$$A \begin{pmatrix} \frac{d m_1}{\Omega dt} \\ \frac{d m_2}{\Omega dt} \\ \frac{d m_3}{\Omega dt} \end{pmatrix} = \begin{pmatrix} F_1 \\ F_2 \\ F_3 \end{pmatrix} + \begin{pmatrix} \frac{L_1}{\Omega^2} \\ \frac{L_2}{\Omega^2} \\ \frac{L_3}{\Omega^2} \end{pmatrix}, \quad (8)$$

where

$$A = \begin{pmatrix} I_o + I_{11} + \frac{2}{3} C_o (m^2 + 2m_1^2) & I_{12} + \frac{4}{3} C_o m_1 m_2 + \Omega T_1 C_o m^2 m_3 & I_{13} + \frac{4}{3} C_o m_1 m_3 - \Omega T_1 C_o m^2 m_2 \\ I_{12} + \frac{4}{3} C_o m_1 m_2 - \Omega T_1 C_o m^2 m_3 & I_o + I_{22} + \frac{2}{3} C_o (m^2 + 2m_2^2) & I_{23} + \frac{4}{3} C_o m_2 m_3 + \Omega T_1 C_o m^2 m_1 \\ I_{13} + \frac{4}{3} C_o m_1 m_3 + \Omega T_1 C_o m^2 m_2 & I_{23} + \frac{4}{3} C_o m_2 m_3 - \Omega T_1 C_o m^2 m_1 & I_o + I_{33} + \frac{2}{3} C_o (m^2 + 2m_3^2) \end{pmatrix} \quad (9)$$

and

$$\left\{ \begin{aligned} F_1 &= (m_3^2 - m_2^2)I_{23} + m_1 m_3 I_{12} - m_1 m_2 I_{13} + m_3 m_2 (I_{22} - I_{33}) \\ &\quad - m_1 \frac{dI_{11}}{\Omega dt} - m_2 \frac{dI_{12}}{\Omega dt} - m_3 \frac{dI_{13}}{\Omega dt} \\ F_2 &= (m_1^2 - m_3^2)I_{13} + m_1 m_2 I_{23} - m_3 m_2 I_{12} + m_1 m_3 (I_{33} - I_{11}) \\ &\quad - m_1 \frac{dI_{12}}{\Omega dt} - m_2 \frac{dI_{22}}{\Omega dt} - m_3 \frac{dI_{23}}{\Omega dt} \\ F_3 &= (m_2^2 - m_1^2)I_{12} + m_2 m_3 I_{13} - m_1 m_3 I_{23} + m_1 m_2 (I_{11} - I_{22}) \\ &\quad - m_1 \frac{dI_{13}}{\Omega dt} - m_2 \frac{dI_{23}}{\Omega dt} - m_3 \frac{dI_{33}}{\Omega dt} \end{aligned} \right. \quad (10)$$

Note that the matrix A may be separated into symmetric and antisymmetric parts: the latter is proportional to $\Omega T_1$ and is responsible for polar wander.

We first neglect the external torque, which could be responsible for the tidal secular deceleration of the axial rotation of the Earth, as we will see in Section 4.6.

Hereafter, the exciting sources are the perturbations of the inertia tensor induced by the internal load $I_{ij}$; the unknowns are the coordinates of the rotation vector $m_1(t)$, $m_2(t)$ and $m_3(t)$.

# 3 VISCOELASTO-GRAVITATIONAL DEFORMATIONS AND MANTLE DENSITY HETEROGENEITIES

This section is focused on the sources of excitation $I_{ij}$ which are the components of the inertia tensor perturbations induced by internal loads within the mantle.

## 3.1 Mantle density heterogeneities and inertia tensor perturbations: theory

The anomaly in mantle mass may be expanded in spherical harmonics, at each radius $r$

$$\Delta\rho(r, \theta, \varphi) = \sum_{n=0}^{\infty} \sum_{m=0}^{n} \left[ \delta\rho_{n}^{cm}(r) \cos m\varphi + \delta\rho_{n}^{sm}(r) \sin m\varphi \right] P_{n}^{m}[\cos(\theta)],$$
where $P_n^m(x)$ are the associated Legendre functions, $\theta$ is colatitude and $\phi$ is longitude. We define $Y_n^{mc} = \cos m\phi P_n^m(\cos\theta)$ and $Y_n^{ms} = \sin m\phi P_n^m(\cos\theta)$ the non-normalized spherical harmonics.

The inertia tensor perturbations $I_{ij}$ depend only on the mantle mass anomalies of degree 0 and 2. However, we must also take into account a delayed effect due to the viscoelastic deformations in our earth model. Based on the MacCullagh theorem (e.g. Munk & MacDonald 1960), we may relate the inertia tensor components to the degree 2 spherical harmonics coefficients of the surface gravitational potential (direct gravitational effect of the load + viscoelastic mass redistribution potential). The inertia tensor perturbations $I_{ij}$ may thus be written, taking into account viscoelastic deformations

$$ 
\left\{
\begin{aligned}
I_{11}(t) &= \frac{4\pi}{5} \int_b^a r^4 \left[ \delta(t) + \left(\frac{a}{r}\right)^2 k_2^2(r, t) \right] * \left[ \frac{1}{3} \delta\rho_2^{c0}(r, t) - \delta\rho_2^{c2}(r, t) \right] dr + \delta I_o(t) \\
I_{22}(t) &= \frac{4\pi}{5} \int_b^a r^4 \left[ \delta(t) + \left(\frac{a}{r}\right)^2 k_2^2(r, t) \right] * \left[ \frac{1}{3} \delta\rho_2^{c0}(r, t) + \delta\rho_2^{c2}(r, t) \right] dr + \delta I_o(t) \\
I_{33}(t) &= -\frac{8\pi}{15} \int_b^a r^4 \left[ \delta(t) + \left(\frac{a}{r}\right)^2 k_2^2(r, t) \right] * \delta\rho_2^{c0}(r, t) dr + \delta I_o(t) \\
I_{12}(t) &= -\frac{8\pi}{5} \int_b^a r^4 \left[ \delta(t) + \left(\frac{a}{r}\right)^2 k_2^2(r, t) \right] * \delta\rho_2^{c2}(r, t) dr \\
I_{13}(t) &= -\frac{4\pi}{5} \int_b^a r^4 \left[ \delta(t) + \left(\frac{a}{r}\right)^2 k_2^2(r, t) \right] * \delta\rho_2^{c1}(r, t) dr \\
I_{23}(t) &= -\frac{4\pi}{5} \int_b^a r^4 \left[ \delta(t) + \left(\frac{a}{r}\right)^2 k_2^2(r, t) \right] * \delta\rho_2^{c1}(r, t) dr
\end{aligned}
\right.
\quad (11) $$

with

$$ \delta I_o(t) = \frac{8\pi}{3} \int_b^a r^4 \left[ \delta(t) + \left(\frac{a}{r}\right)^2 k_0^2(r, t) \right] * \delta\rho_0^{c0}(r, t) dr, \quad (12) $$

where a and b are, respectively, the surface and core radii. Inertia changes due to mantle density heterogeneities act directly on the planet through the $\delta(t)$ Dirac function. However, they also have a delayed effect due to viscoelastic deformations: these effects are taken into account through the time convolution with degree 0 and degree 2 internal loading Love number $k_0^2(r, t)$ and $k_2^2(r, t)$. These Love numbers express the contribution of viscoelastic deformations induced by internal loads within the Earth to the inertia tensor perturbations. We define the viscous Love numbers $k_0^f(r)$ and $k_2^f(r)$ as the fluid limit of the viscoelastic Love numbers: this is a good approximation when the timescale of internal loads greatly exceeds the viscoelastic relaxation times. In this case, the temporal convolution in (11) disappears and the inertia tensor perturbations become, for example for $I_{13}$ and $\delta I_o$

$$ 
\left\{
\begin{aligned}
I_{13}(t) &= -\frac{4\pi}{5} \int_b^a r^4 \left[ 1 + \left(\frac{a}{r}\right)^2 k_2^f(r) \right] \delta\rho_2^{c1}(r, t) dr \\
\delta I_o(t) &= \frac{8\pi}{3} \int_b^a r^4 \left[ 1 + \left(\frac{a}{r}\right)^2 k_0^f(r) \right] \delta\rho_0^{c0}(r, t) dr
\end{aligned}
\right.
\quad (13) $$

The internal loading degree 2 viscous Love number has been widely investigated in the literature on the incompressible Earth. Degree 0 is related to mass conservation and consequently, we cannot assume anymore an incompressible earth model: as a matter of fact, this assumption would forbid any volume variation. However such variations are necessary to ensure mass conservation for the Earth. The computation of these internal loading Love numbers for the degrees 0 and 2 is detailed in Appendix A.

For degree $n > 1$ the viscous geoid kernels are not significantly influenced by the compressibility of the mantle (e.g. Forte & Peltier 1991; Corrieu *et al.* 1995, and more recently Cambiotti *et al.* 2009). So we will assume, for simplicity an incompressible planet when investigating degree 2 deformations. The stratification of the mantle plays an important role in the computation of the degree 2 geoid kernel. Geodynamic and convection studies establish that an isoviscous homogeneous earth model is far from reality (Ricard *et al.* 1993a; Richards *et al.* 1999). So we will use hereafter a simple four-layered model consisting of a viscoelastic lithosphere, a viscoelastic upper and lower mantle (we assume that the lithosphere and the whole mantle have the same density $\rho^m$), and an inviscid fluid core (with a density $\rho^c$). The geometrical and physical parameters of our four-layered earth model are given in Table 1. VR is the ratio of lower to upper mantle viscosities. The associated degree zero and two kernels are plotted in Fig. 1, for different VR ranging from 5 to 100. The thickest solid line is for VR = 40.

Our degree 2 kernel for VR = 40 (solid line) is close to the degree 2 geoid kernel computed in studies based on the best fit between the observed geoid and the one induced by mantle density heterogeneities derived from tomography and/or geodynamic models (e.g. Ricard *et al.* 1993a; Steinberger 2000; Marquart *et al.* 2005). So, we will use this viscosity ratio for our reference model. For the degree $n = 2$, and

<table><caption>**Table 1.** Geometrical and physical parameters of the four-layer reference Earth model (PREM-averaged values).</caption><thead><tr><th>Radius (km)</th><th>Density (kg m<sup>-3</sup>)</th><th>Rigidity (Pa)</th><th>Viscosity (Pa s)</th></tr></thead><tbody><tr><td>6271 &lt; r &lt; a</td><td>4414</td><td>1.66 &times; 10<sup>11</sup></td><td>1.1 &times; 10<sup>22</sup></td></tr><tr><td>5701 &lt; r &lt; 6271</td><td>4414</td><td>1.66 &times; 10<sup>11</sup></td><td>10<sup>21</sup></td></tr><tr><td>3480 &lt; r &lt; 5701</td><td>4414</td><td>1.66 &times; 10<sup>11</sup></td><td>VR &times; 10<sup>21</sup></td></tr><tr><td>0 &lt; r &lt; 3480</td><td>12 420</td><td>0</td><td>0</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Degree zero (dashed line) and two (solid line) kernels [$\frac{r_o^2}{a^2} + k'_n(r_o)$] as a function of $r_o$, for different viscosity ratio VR ranging from 5 to 100. The thickest solid line is for VR = 40.</figcaption></figure>

<table><caption>**Table 2.** Rotational time $T_1$ in kyr for different viscosity ratio between upper and lower mantle.</caption><thead><tr><th>Viscosity ratio VR</th><th>5</th><th>10</th><th>20</th><th>30</th><th>40</th><th>50</th><th>60</th><th>70</th><th>80</th><th>90</th><th>100</th></tr></thead><tbody><tr><th>T<sub>1</sub> (kyr)</th><td>3.29</td><td>6.10</td><td>10.98</td><td>15.08</td><td>18.59</td><td>21.63</td><td>24.29</td><td>26.66</td><td>28.78</td><td>30.69</td><td>32.42</td></tr></tbody></table>

$VR = 40$, a positive mass anomaly located in the upper part of the mantle will involve a positive geoid whereas a positive mass anomaly in the lower mantle will involve a negative geoid; the degree 2 geoid kernel changes sign at a depth 1200 km.

The degree 0 kernel does not depend on the viscosity ratio (see Appendix A). For the degree $n = 0$, a positive mass anomaly within the mantle will involve a positive perturbation of the moment of inertia, whatever the depth of the load.

We have also computed for these earth models the degree 2 tidal Love number, and especially the rotational time $T_1$ which increases with the viscosity ratio (see Ricard et al. 1993a) for more details): its value, for incompressible earth model, ranges from 3 to 32 kyr (Table 2) so that the quantity $\Omega T_1$ is around $10^7$ that is to say much larger than 1. In matrix A, the terms in $\Omega T_1$ corresponding to the non-symmetric part and directly related to polar wander are dominant.

In the next section, we test a simple model for mantle mass anomalies: a cold blob sinking vertically within the mantle in order to compare the degree 0 and 2 contributions of the inertia tensor perturbations.

## 3.2 Plates sinking into the mantle

In this part, we want to investigate a mantle mass anomaly (with a mass $\delta M$) which remains constant $\frac{\delta M}{\delta t} = 0$ but which moves within the mantle: for example, a cold blob sinking vertically into the mantle. Fig. 2 is a simple cartoon where the blob is modelled as a spherical cap with radius $L$ and thickness $H$, sinking vertically within the mantle along an axis with colatitude $\theta_o$ and longitude $\phi_o$. This shape is highly unrealistic but because it involves simple solutions for the mantle mass anomaly and the associated inertia perturbations, we use it as a pedagogical example. Assuming that $L \ll a$, the mass of this blob is $\delta M = \Delta\rho\pi L^2 H$, where the density contrast $\Delta\rho$ between the cold blob and the surrounding mantle is assumed to be constant with time.

The inertia tensor perturbation induced by the degree 2 mass anomaly associated with the blob at radius $r_o$ may be written

$$ 
\left\{
\begin{aligned}
I_{11} &= \frac{1}{3}\delta M a^2 \left(Y_2^{0c}(\theta_o, \varphi_o) - \frac{1}{4}Y_2^{2c}(\theta_o, \varphi_o)\right) \left[\frac{r_o^2}{a^2} + k'_2(r_o)\right] + \delta I_o \\
I_{22} &= \frac{1}{3}\delta M a^2 \left(Y_2^{0c}(\theta_o, \varphi_o)r + \frac{1}{4}Y_2^{2c}(\theta_o, \varphi_o)\right) \left[\frac{r_o^2}{a^2} + k'_2(r_o)\right] + \delta I_o \\
I_{33} &= -\frac{2}{3}\delta M a^2 Y_2^{0c}(\theta_o, \varphi_o) \left[\frac{r_o^2}{a^2} + k'_2(r_o)\right] + \delta I_o \\
I_{12} &= -\frac{1}{6}\delta M a^2 Y_2^{2s}(\theta_o, \varphi_o) \left[\frac{r_o^2}{a^2} + k'_2(r_o)\right] \\
I_{13} &= -\frac{1}{3}\delta M a^2 Y_2^{1c}(\theta_o, \varphi_o) \left[\frac{r_o^2}{a^2} + k'_2(r_o)\right] \\
I_{23} &= -\frac{1}{3}\delta M a^2 Y_2^{1s}(\theta_o, \varphi_o) \left[\frac{r_o^2}{a^2} + k'_2(r_o)\right]
\end{aligned}
\right.
\quad (14) $$

with $\delta I_o = \frac{2}{3}\delta M a^2 \left[\frac{r_o^2}{a^2} + k'_o(r_o)\right]$.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Cartoon where the blob is modelled as a spherical cap with a radius *L* and a thickness *H* sinking vertically within the mantle along the axis with a colatitude $\theta_o$ and a longitude $\phi_o$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Moment of inertia $I_{33}/I_o$ (left-hand panel) and $\delta I_o/I_o$ (right-hand panel). The solid thick line is for VR = 40.</figcaption></figure>

At $t=0$ the blob is located at the top of the lithosphere and then it sinks within the upper mantle at a constant velocity about 5 cm yr$^{-1}$. When the blob crosses the 670 km discontinuity, its downward velocity is assumed to decrease by a factor that may be written as the logarithm of the viscosity ratio between lower and upper mantle (Ricard *et al.* 1993b).

We use $L = 1000$ km, $H = 100$ km, and $\Delta\rho = 80$ kg m$^{-3}$. We assume that the blob sinks along the axis ($\theta_o = 105^\circ, \phi_o = 208^\circ$), and we compute $I_{33}$ and $\delta I_o = \frac{2}{3}\delta Ma^2[\frac{L_o^2-a^2}{a^2} + k'_0(r_o) - k'_0(a)]$, for viscosity ratios VR ranging from 5 to 100 (Fig. 3).

Note that $\delta I_o$ is larger than $I_{33}$. The degree 0 and 2 mantle density anomalies have the same order of magnitude, but because of the deformations at the CMB and at the surface, the degree 2 geoid kernel is smaller than the degree 0 kernel. This point implies that LOD will essentially depend on degree 0 mantle mass anomalies, whereas TPW will depend on degree 2.

We now compute how the rotational axis is perturbed by this sinking blob, solving the set of eq. (8). We first plot in Fig. 4 the temporal evolution of the latitude and the longitude of the pole and of course obtain the well-known result that the rotational axis tends to be in the direction of the maximal principal inertia axis (PIA): for $t < 52.5$ Myr, the blob is in the upper part of the mantle at a depth $< 1200$ km where the degree 2 geoid kernel is positive, and the maximal PIA is in a plane, with a colatitude about $15^\circ$, perpendicular to the axis of the blob. For $t > 52.5$ Myr, the blob is at a depth where the geoid kernel is negative and consequently the maximal PIA is the direction of the axis of the blob. At $t = 184$ Myr, the blob arrives at the CMB where it is isostatically compensated and there is no more drift of the pole of rotation. The cuspidal point in the latitude is associated to a change by $180^\circ$ in longitude: this expresses the fact that the pole crosses the geographical North axis.

The LOD variation with respect to the initial value (at $t=0$) is plotted in Fig. 5 (left-hand column) for various viscosity ratios VR.

Note that for $VR \ge 50$ the LOD variation is close to $\frac{2\pi}{\Omega}\frac{\delta I_o(t)}{I_o}$ (right-hand column) that is to say is mainly due to the degree 0 mantle mass anomaly.

For $VR < 50$, there are fluctuations induced by the degree 2 mantle mass anomaly, although the dominant term results from the degree zero mantle density heterogeneities.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Latitude (left-hand panel) and longitude (right-hand panel) of the pole of rotation (i.e. the intersection between the axis of rotation and the surface of the planet) for various viscosity ratio. In thick solid line, VR = 40.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** LOD variation (left-hand panel) and $\frac{2\pi \delta I_o(t)}{\Omega I_o}$ (right-hand panel) with respect to the initial value (at $t = 0$) induced by the sinking blob, for various viscosity ratio VR ranging from 10 to 100: in thick solid line, VR = 40, in dashed line VR = 50 and in long-dashed line VR = 30.</figcaption></figure>

There is a geological acceleration of the angular velocity of the Earth, and thus a geological decrease in the LOD with an order of magnitude about 0.01 – 0.02 µs yr⁻¹.

Before investigating the perturbation of the rotation vector induced by inertia tensor perturbations for a realistic mantle model with mass anomalies related to subduction and to large-scale upwellings, let us recall the simple problem of the secular cooling of the Earth.

### 3.3 Secular cooling

Mantle convection cools the Earth and the amount of heat that escapes can be measured by studying the heat flux through the ocean floor and continents. The total heat loss of the Earth is about $40 \times 10^{12}$ W (TW). Subtracting radiogenic heat production (about 20 TW) from this total heat loss yields 20 TW due to secular cooling (that lowers mantle temperature and extracts heat from the core through the CMB).

We define T the temperature.

Secular cooling amounts to about $\frac{dT}{dt} = -100$ K Gyr⁻¹. As a consequence, the average density of the Earth increases by about $\frac{d\rho}{dt} = -\alpha\rho\frac{dT}{dt} = 16.5 \text{ kg m}^{-3} \text{ Gyr}^{-1}$ with $\alpha = 3 \times 10^{-5} \text{ K}^{-1}$ the average thermal expansion coefficient of the mantle (for a review, see Jaupart et al. 2009).

To ensure the conservation of the mass of the Earth, a contraction of the planet results, that is, a decrease of its radius $\frac{da}{dt}$ such that

$$ \frac{da}{dt} = \frac{\alpha}{3} a \frac{dT}{dt} = -6371 \text{ m Gyr}^{-1} < 0. $$

We assume the planet to be homogeneous, in this simple model. The variation of its moment of inertia associated with such cooling will be

$$ \begin{aligned} \frac{d\delta I_o}{dt} &= \frac{8\pi}{15} \left[ \frac{d\rho}{dt} a^5 + 5a^4 \rho \frac{da}{dt} \right] \\ &= I_o \left[ -\alpha \frac{dT}{dt} + \frac{5}{3}\alpha \frac{dT}{dt} \right] \\ &= \frac{2}{3}\alpha I_o \frac{dT}{dt} < 0. \end{aligned} \qquad (15) $$
The effect of the decrease in radius is larger than the effect of the increase in density, so that the moment of inertia decreases with an order of magnitude about: $\frac{\delta I_o}{dr} = -2 \times 10^{-3} I_o$ per Gyr.

Because of the conservation of the axial angular momentum of the Earth there is an acceleration of the angular velocity associated with this decrease of moment of inertia and the associated variation of the LOD is

$$ \frac{2\pi}{\Omega} \frac{\delta I_o(t)}{I_o} = -172 \, \text{s Gyr}^{-1} $$

it is to say about $-0.2 \, \text{s Myr}^{-1}$.

In the next sections, we will compute the perturbation of the rotation vector induced by inertia tensor perturbations corresponding to mantle mass anomalies related to subduction and to large-scale upwellings, on a timescale of a hundred million years. Secular cooling will be included in this computation as well as variations at shorter timescale (a few tens of millions years).

# 4 GEOPHYSICAL APPLICATION

## 4.1 Some remarks about the conservation of the mass of the mantle

The ocean floor is formed at ridges where the upwelling mantle is partially molten. As it moves away from the ridge, the oceanic plate cools and becomes denser. At a certain distance which depends on the velocity and age of the plate, density becomes too large for the plate to be stable at the surface, and the plate then sinks into the mantle at a subduction zone. Extraction of hot mantle at ridges, plate cooling, injection of the cold plate at subduction zones forms a cycle in which heat is removed from the Earth (see previous section).

We define $\delta m^S(t)$ the mass injected at subduction zones.

The above cycle conserves the total mass of the Earth, so that the mass extracted at ridges is equal to $-\delta m^S(t)$.

Sinking plates are colder than the surrounding mantle and consequently induce a positive excess of mass: the density contrast of the slab with the surrounding mantle is assumed to remain constant with depth and equal to $80 \, \text{kg m}^{-3}$ (Ricard et al. 1993b), until the depth where the plates stop sinking. We define $\delta m^m(t)$ the total mass in the mantle

$$ \delta m^m(t) = 4\pi \int_b^a \delta\rho_0^{c0}(r, t)r^2 dr. \qquad (16) $$

We define $\delta m^C(t)$ the total mass which has been introduced into the Earth by subduction and has reached the CMB, at a given time $t$. When the slab arrives at the bottom of the mantle, it cannot stay since additional matter is arriving from above. Mass conservation requires that for each element of mass arriving at a given depth, an equal amount must move up to compensate for it. Consequently, we do not know how the cold mass is distributed at the CMB, but it stays there for an undetermined amount of time before to disappearing due to heating (because of radioactive elements). We assume that this cold mass cools the whole Earth inducing an increase in the mean density of about $\frac{3}{4\pi} \frac{\delta m^C(t)}{a^3}$.

Let us now compute the perturbations of the moment of inertia associated with these three sources of degree 0 mass anomalies

(i) The perturbation of the moment of inertia induced by the extraction of hot mantle at the ridges is

$$ \delta I_o^S(t) = -\frac{2}{3} \delta m^S(t) a^2. \qquad (17) $$

It is exactly compensated by the inertia perturbation induced by the injection of cold plates at subduction zones.

(ii) The perturbation of inertia due to mantle mass anomalies resulting from the sinking plate $\delta I_o^m$ may be written

$$ \delta I_o^m(t) = \frac{8\pi}{3} \int_b^a r^4 \left[ \delta(t) + \left(\frac{a}{r}\right)^2 k_0'(r, t) \right] * \delta\rho_0^{c0}(r, t) dr. \qquad (18) $$

(iii) The perturbation of inertia due to mass anomaly resulting from the cooling of the mantle $\delta I_o^c$

$$ \delta I_o^c(t) = \frac{2}{5} \delta m^C(t) a^2. \qquad (19) $$

We next investigate the influence of upwelling domes and subducted plates sinking into the mantle on LOD, using very simple models of mantle dynamics.

## 4.2 A simple model for the large-scale pattern of mantle dynamics

We first present our simple model for the large-scale pattern of mantle dynamics. This model is described in detail in (Rouby et al. 2010). The model is able to explain the observed present-day geoid and the wander of the rotation axis since 120 Ma, based on the distribution of mantle density heterogeneities. In this section, we just point out its main features.

### 4.2.1 Domes

Laboratory study (Davaille 1999) of thermochemical convection in a fluid with stratified density and viscosity has shown that a 'doming' regime may be observed for a range of parameters plausible for the Earth. In this regime, the domes oscillate vertically within the mantle
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Cartoon representing the construction of a density heterogeneities model for the Earth based on subduction history from Lithgow-Bertelloni & Richards (1998). The reference frame used in this cartoon is related to the surface of the oceanic plate and consequently the sinking of the plate does not look vertical. In a reference frame related to the mantle the plate sinks vertically with the velocity $V_c$.</figcaption></figure>

with very long pseudo-periodicities ranging from some hundred million years to 1 Gyr. Such domes could be responsible for the superswells observed at the Earth's surface (e.g. Lithgow-Bertelloni & Silver 1998; Gurnis et al. 2000; Tan et al. 2002). From experimental study as well as palaeomagnetic observations (Davaille et al. 2005), the characteristic timescale of the doming regime seems larger than the timescale of subduction. Consequently, we have assumed that the spatial distribution of mass anomalies in the lower mantle associated with these large-scale upwellings have been constant in time for the past 120 Myr, and we have modelled them using the present-day tomography of the lower mantle, assuming a constant density contrast $\Delta\rho = -50 \text{ kg m}^{-3}$ between the 'hot' domes and the surrounding mantle.

We introduce inertia tensor perturbations, $I^D$, in order to take into account this non-hydrostatic internal load which varies at a timescale longer than that of subduction.

### 4.2.2 Subductions

Subducted lithosphere is a major component of mantle density heterogeneities. We start with the model of mantle density heterogeneities derived by Ricard et al. (1993b), which uses plate-motion reconstructions under the assumption that subducted slabs sink vertically into the mantle (Fig. 6), first with the observed surface velocity. When they cross the upper/lower mantle interface, the slabs are assumed to be folded or thickened and the velocity is decreased by a factor $\ln VR$ of about 4.

Lithgow-Bertelloni & Richards (1998) modelled the temporal evolution of mantle density heterogeneities associated with the large-scale pattern of plate tectonic motions since 120 Ma, assuming that the plates sink down to the CMB. However, the tomographic images show that subducted plates do not reach systematically the CMB (Fukao et al. 2001) and we have investigated the influence to the geoid of the depth where the plates stop to sink (Rouby et al. 2010).

### 4.2.3 Our simple model

Finally, our 'preferred' model of mantle density heterogeneities, i.e. the model for which we obtain the best variance between the present-day computed and observed geoid (for degrees 2 to 12), assumes time-independent large-scale upwellings within the lower mantle and that all plates sink down to the CMB except for North America that stops at 2250 km and South America at about 800 km (Rouby et al. 2010).

We use this model to compute the time evolution of these internal loads; taking into account viscoelastic deformations, we compute the inertia tensor perturbations and then the LOD variations.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Rate of mass injected in the subduction zones, $\frac{d\delta m^S(t)}{dt}$ since 200 Ma, in $10^{18}$ kg Myr$^{-1}$.</figcaption></figure>

## 4.3 Inertia tensor perturbations

### 4.3.1 Degree 0

We first compute the rate of mass injected into subduction zones, $\frac{d\delta m^S(t)}{dt}$, for the 12 different stages of the plate-motion reconstruction model of Lithgow-Bertelloni & Richards (1998) since 200 Ma (Fig. 7). The average rate is about $27 \times 10^{18}$ kg Myr$^{-1}$ in the past 200 Myr.

In order to be self-consistent over the period under study, the model must consider the material subducted during a previous period that amounts to the characteristic time needed for a plate to sink through the whole mantle. The time evolution of mantle density heterogeneities associated with the large-scale pattern of plate tectonic motions since 120 Myr is then modelled assuming that plates have been sinking down to the CMB for 200 Myr. Initial time is then taken to be 120 Ma. We plot, in Fig. 8(a), the total mass subducted within the mantle $\delta m^S(t)$, the mass of the slabs sinking within the mantle $\delta m^m(t)$ and the mass $\delta m^C(t)$. This cold mass anomaly $\delta m^C(t)$ is equivalent to an increase of the mean mantle density plotted in Fig. 8(b); it is about $0.020 \text{ kg m}^{-3} \text{ Myr}^{-1}$, that is to say a cooling of the mantle by about

$$ \frac{dT}{dt} = -\frac{1}{\alpha \rho_o} \frac{d\rho}{dt} \approx -0.12 \text{ K Myr}^{-1} $$

a value which is consistent with the one given in Section 3.3 from Jaupart et al. (2009).

Using the eqs (17)–(19), we now compute and plot, in Fig. 9, the moments of inertia $\delta I_o^S$, $\delta I_o^m$, $\delta I_o^C$ and the sum of the three contributions (solid lines) induced by the degree 0 mass anomalies.

The order of magnitude of the total moment of inertia perturbation is about $10^{-4}$, with a rate of change of about $-4 \times 10^{-6}$ per Myr.

We must add to this time-variable moment a constant term induced by the domes: for our simple model, its value is about $-2.6 \times 10^{-5}$.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Left-hand panel: total mass subducted within the mantle $\delta m^S(t)$ (solid line), mass of the slabs sinking within the mantle $\delta m^m(t)$ (dashed line) and mass $\delta m^C(t)$ (long dashed line). Right-hand panel: increase of the mean mantle density in kg m$^{-3}$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Perturbation of the moment of inertia $\frac{\delta I_o^S}{I_o}$ (dotted line), $\frac{\delta I_o^M}{I_o}$ (dashed line), $\frac{\delta I_o^C}{I_o}$ (dot-dashed line) and the sum of the three contributions (solid line).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Six components of the inertia tensor perturbations induced by the degree 2 mantle density heterogeneities.</figcaption></figure>

### 4.3.2 Degree 2

We now compute, and plot in Fig. 10, the components of the inertia tensor perturbations induced by the degree 2 mantle density heterogeneities, taking into account both domes and subductions. Their orders of magnitude are about $10^{-5}$.

## 4.4 True polar wander and length of day variations

To compute how the rotational axis is perturbed by mantle dynamics, we solve the set of eq. (8).

As already found in Rouby *et al.* (2010), we show that there is a large motion of the pole of rotation at the Earth's surface during the last 120 Ma towards the present North pole (there is a large filled circle every 10 Myr), that is plotted in Fig. 11(a). Because of the errors in modelling the mantle convection inertia tensor, the present-day rotation axis fails to coincide with the North pole by about 3°.

The LOD variation with respect to the initial value (at $t = -120$ Ma) is plotted in Fig. 11(b), since 120 Ma: it is about a few tens of seconds. Because the moment of inertia has been decreasing for the past 120 Ma (Fig. 9), there is an acceleration of the angular velocity of the Earth and consequently a decrease in the LOD.

We plot in Fig. 12 the time-derivative of LOD $-\frac{2\pi}{\Omega}\frac{\dot{m}}{m^2}$ since 120 Ma. The rate of LOD variation induced by our simple mantle density heterogeneities is about $-0.4 \text{ s Myr}^{-1}$. Superimposed on this linear trend, there are fluctuations around $0.1-0.2 \text{ s Myr}^{-1}$ on a timescale of a few tens of millions of years. In conclusion, there are fluctuations induced by the degree 2 mantle mass anomaly, although the dominant term results from the degree zero mantle density heterogeneities.

## 4.5 Mantle circulation models and LOD variations

Before comparing these results with the others sources of perturbations of LOD, we would like to improve the problem of degree 0 mantle mass anomaly (i.e. the problem of mass conservation), with the use of a 3-D spherical Mantle Circulation Model (MCM) with consistency between downwelling and up-welling flows.

We use a numerical simulation from the 'Tectosphère' code developed by Monnereau & Quere (2001). For this test, the model has no plate and boundary conditions are free-slip at the top and at the bottom of the mantle. There is 50 per cent internal heating. The viscosity ratio between lower and upper mantle is about 30.

This model assumes a thermal equilibrium and there is no long-term evolution of the temperature. As shown in Section 3.3, the fluid heats up or cools down out of thermal equilibrium on the long-term scale, and there would be a change in volume (the Earth would contract) in order to keep the total mass of the Earth constant and an associated LOD variation of about $-0.2 \text{ s Myr}^{-1}$.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Left-hand panel: polar wander using an orthographic projection, with the longitude varying from 0° to 360° and the latitude from 0° to 90°. The centre of the figure is the North pole and the parallels of latitude are plotted at 15° intervals; Right-hand panel: length of day variation since 120 Ma.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Time-derivative of the Length of Day variation in s Myr⁻¹ since 120 Ma.</figcaption></figure>

Here we investigate the mantle dynamics on a timescale of a few hundred millions years. In MCMs, all the extensive variables (mass and energy) are conserved. But, in the principle of the equations, it is still assumed that, within the uncertainties of the numerical treatment, all the variations in density are neglected except for the term [noted $\delta\rho_0^c(r, t)$] in the Navier-Stokes equations (the so-called Boussinesq approximation). This means that the volume change over time is neglected, because mass does not change. Consequently, for the calculation of moment of inertia induced by $\delta\rho_0^c(r, t)$, we must add a term related to the associated variation of the radius of the planet that ensures mass conservation.

If $\delta m(t) = 4\pi \int_0^a \delta\rho_0^c(r, t) r^2 dr$ is the long-term excess or default of mass, we compute the variation of the Earth's radius $\delta a = -\frac{\delta m(t)}{4\pi a^2 \rho_0(a)}$ which ensures the conservation of total mass and we add to the computation of the moment of inertia a corrected term $\frac{8\pi}{3} \rho_0(a) \delta a a^4 = -\frac{2}{3} \delta m(t) a^2$.

We find a decrease of the mass $\frac{\delta m(t)}{M}$, about $-1.32 \times 10^{-7} \text{ Myr}^{-1}$, that is to say an increase of the average temperature of the mantle by about $4.4 \times 10^{-3} \text{ K Myr}^{-1}$. The planet therefore undergoes some dilatation and its radius increases by about 0.3 m Myr⁻¹. Note that this effect is several orders of magnitude smaller than the one induced by secular cooling.

The perturbations of the moment of inertia and of the inertia tensor components are plotted respectively in Figs 13(a) and (b). The moment of inertia induced by degree 0 mantle mass anomalies is larger than the inertia tensor components induced by degree 2 mass anomalies. It grows with a rate of about $8.9 \times 10^{-8} \text{ Myr}^{-1}$.

Solving the equations for the conservation of angular momentum, we compute the polar wander and its velocity (Figs 13c and d) and the LOD variation and its time derivative (Figs 13e and f). The amplitude of polar wander is large and may reach 2 degrees per Myr, an order of magnitude larger than the observed one.

The LOD variation is mainly created by the moment of inertia and is about ten ms Myr⁻¹. Superimposed on the long-term LOD variations, there are fluctuations at the smaller timescale (less than 100 Myr) mainly induced by degree 2 mantle mass anomalies, with a similar order of magnitude.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Mantle circulation models and LOD variations. (a) Temporal evolution of the moment of inertia $\delta I_o$ normalized to $I_o$. (b) Temporal evolution of the 6 components of the inertia tensor $I_{ij}$ normalized to $I_o$. (c) Temporal evolution of the latitude (solid line) and longitude (dotted-dashed line) of the rotational axis. (d) LOD variations as a function of time. (e) Velocity of the polar wander. (f) Time derivative of LOD in seconde per Myr.</figcaption></figure>

Our approach is different from the study of Machetel & Thomassot (2003) who investigated the influence of the catastrophic mantle mass transfers due to avalanches at the 670 km depth discontinuity on LOD, taking into account only the degree 0 coefficient in the spherical harmonic expansion of mantle density heterogeneities, that is to say only the radial perturbations. In this case, since there is no possible polar wander associated with degree 0 mantle density heterogeneities, the authors used the linearized form of the angular momentum theorem (Liouville equations) to compute the LOD variations. In their paper they did not discuss the conservation of mass and they did not take into account the degree 2 coefficients of the spherical harmonic expansion of the mass anomaly. These induce smaller LOD variations than degree 0 but nevertheless involve non-negligible fluctuations around a linear trend.

## 4.6 Frictional tidal torque between the Moon and the Earth

So far, we have neglected external torque that might act on the Earth, taking only into account changes in the inertia tensor of the Earth. In this section, we would like to take these two effects simultaneously into account, using our non-linear approach, in order to calculate the associated variation in LOD.

The secular deceleration of the Earth's axial rotation and the evolution of the distance between the Moon and the Earth result from the braking torque due to lunar tidal forces acting on the Earth. The dissipated energy is due to the viscosity of the mantle (Zschau 1978) and to the friction at the bottom of the oceans (Lambeck 1980). These two causes of dissipation have existed since the beginning of the evolution of the Earth-Moon system. Nevertheless, it is difficult to quantify their relative magnitudes, since they depend not only on the time evolution of the viscosity of the Earth but also on the friction at the bottom of the oceans related to palaeogeographical modifications.

On a secular timescale, the effect of such a frictional torque ($L_3^f$) on LOD variations is usually studied with the axial linearized Liouville equation

$$C\dot{\omega}_3 = L_3^f, \qquad (20)$$

where C is the axial moment of inertia, leading to a secular deceleration of the Earth's axial rotation by about $\frac{L_3^f}{C} = -5.98 \times 10^{-22} \text{ rad s}^{-2}$, that is to say a LOD variation about $2 \times 10^{-3} \text{ s per century}$.

The angular deceleration rate deduced from studies of bivalves and stromatolites during the last 500 Myr is more than $-5.2 \times 10^{-22}$ and less than $-6.6 \times 10^{-22} \text{ rad s}^{-2}$ (Lambeck 1978; Bursa 1982, 1990), that-is-to say seems similar to the present rate. The amplitude of the tidal torque is inversely proportional to the cube of the Earth-Moon distance, and thus was only -4 per cent larger 120 Ma. Consequently, at the timescale of 100 Myr, we assume that the amplitude of this frictional torque remains constant.
Its direction is always aligned with the Earth's rotation axis and consequently, in our present-day reference frame, it may be written as

$$ \left\{ \begin{aligned} L_1 &= \frac{\omega_1}{\sqrt{\omega_1^2 + \omega_2^2 + \omega_3^2}} L_3^f \\ L_2 &= \frac{\omega_2}{\sqrt{\omega_1^2 + \omega_2^2 + \omega_3^2}} L_3^f \\ L_3 &= \frac{\omega_3}{\sqrt{\omega_1^2 + \omega_2^2 + \omega_3^2}} L_3^f \end{aligned} \right. \qquad (21) $$

Equatorial components of the torque appear because the rotation axis varies at the surface of the planet in our reference frame (TPW).

We solve the non-linear set of eqs (8), taking into account the mantle density heterogeneities described in the previous section, together with the frictional tidal torque, in order to investigate if the latter may induce polar wander due to possible coupling between equatorial and axial components.

We find that:

(i) The polar wander differs only by about 2 per cent depending on whether one takes into account or not the frictional torque.

(ii) The dominant term in the LOD variation is a secular term with an amplitude about 20 s Myr⁻¹.

## 5 COMPARISON WITH OTHER EFFECTS

The LOD varies under the action of external torques as well as of time evolution of mass redistribution on or in the Earth, on different timescales. We compare our results with other effects varying at decadal or secular timescale. Note that the present level of accuracy of LOD observations, using Very Long Baseline Interferometry, is about 10 µs.

Which sources are able to generate observable variations of LOD on the secular timescale (Chao 1994; De Viron et al. 2002)?

(i) Secular deceleration of the Earth's axial rotation and the evolution of the distance between the Moon and the Earth resulting from the braking torque due to lunar tidal forces acting on the Earth: 20 µs yr⁻¹ that-is-to-say 20 s Myr⁻¹ (Lambeck 1980; Christodoulidis et al. 1988).

(ii) Secular acceleration related to the postglacial rebound: a few µs yr⁻¹. This value strongly depends on models for the last deglaciation and on the viscosity profile within the mantle (e.g. Mitrovica & Peltier 1993).

(iii) Effect of melting glaciers on the Earth's rotation: using yearly mass balance data for 85 glaciers in 13 of the 31 mountain glacier systems designated in (Meier 1984)'s paper, the contribution of glaciers to the total excess LOD signal is about 0.4 µs yr⁻¹ (Trupin et al. 1992).

(iv) Worldwide sea-level rise: 0.4 µs yr⁻¹

(v) Decadal and secular LOD fluctuations attributed to exchanges of angular momentum between the core and the mantle: from recent models of secular variation of the magnetic field, there are variations about a few ms in a few years (Greff-Lefftz et al. 2004).

(vi) Atmospheric and oceanic excitation of decadal-scale Earth orientation variations during 1949–2002 studied using atmospheric and oceanic angular momentum series computed from the NCEP/NCAR reanalysis system is found to be about 14 per cent of that observed, that is to say with peak-to-peak variations of about 0.3 ms (Gross et al. 2005).

## 6 CONCLUSION

We have quantified the length of day variations induced by mantle density heterogeneities, and especially discussed the influence of the degree 0 versus degree 2 coefficients of the spherical harmonic expansion of mantle mass anomalies. The problem of LOD is strongly related to the conservation of mass within the mantle.

We have found that the time-variable mantle density heterogeneities associated with large-scale patterns of plate tectonic motions since 120 Ma (Lithgow-Bertelloni & Richards 1998) can perturb the LOD by about 0.4 µs yr⁻¹, that is to say with an order of magnitude comparable with the effect of the secular cooling and the associated contraction of the Earth, but one order of magnitude smaller than the effect of the last deglaciation. Superimposed on this linear trend, we have shown that there are fluctuations around 0.1–0.2 s Myr⁻¹ on a timescale of a few tens of millions of years.

On the geological scale, this perturbation is too small to be constrained by studies of bivalves and stromatolites as done during the last 500 Myr to compute the tidal secular deceleration (Lambeck 1978; Bursa 1990; Varga 2006).

## ACKNOWLEDGMENTS

This study is an IPGP contribution. I thank Marc Monnereau for providing numerical simulation results and for discussion on the mantle circulation model. I gratefully acknowledge Anne Davaille for helpful discussions and Vincent Courtillot for remarks on the original manuscript.
REFERENCES

Alterman, Z., Jarosch, H. & Pekeris, C.H., 1959. Oscillation of the Earth, *Proc. R. Soc. Lond.*, **A252**, 80-95.

Bursa, M., 1982. On some topical problems of the dynamics of the Earth-Moon system, in *Tidal Friction and the Earth's Rotation*, Vol. II, pp. 19-28, eds Brosche, P. & Sundermann, J., Springer-Verlag, Berlin.

Bursa, M., 1990. The variation in *J*₂ and in the moments of inertia: satellite results and consequences for the angular momentum budget of the Earth-Moon-Sun system, in *Earth's Rotation from Eons to Days*, pp. 52-57, eds Brosche, P. & Sundermann, J., Springer-Verlag, Berlin.

Cambiotti, G., Barletta, V.R., Bordoni, A. & Sabadini, R., 2009. A comparative analysis of the solutions for a Maxwell Earth: the role of the advection and buoyancy force, *Geophys. J. Int.*, **176**(3), 995-1006.

Chao, B.F., 1994. The Geoid and Earth's rotation, in *Geoid and its Geophysical Interpretations*, pp. 285-298, eds Vanicek, P. & Christou, N.T., CRC Press, Boca Raton, FL.

Chen, J., 2005. Global mass balance and the length-of-day variation, *J. geophys. Res.*, **110**, B08404, doi:10.1029/2004JB003474.

Christodoulidis, D.C., Smith, D.E., Williamson, R.G. & Klosko, S.M., 1988. Observed tidal braking in the Earth-Moon-Sun System, *Geophys. J. R.*, **93**, 6216-6236.

Corrieu, V., Thoraval, C. & Ricard, Y., 1995. Mantle dynamics and geoid Green functions, *Geophys. J. Int.*, **120**(2), 516-523.

Dahlen, F.A. & Tromp, J., 1998. *Theoretical Global Seismology*, Princeton University Press.

Davaille, A., 1999. Simultaneous generation of hotspots and superswells by convection in a heterogeneous planetary mantle, *Nature*, **402**, 756-760.

Davaille, A., Stutzmann, E., Silveira, G., Besse, J. & Courtillot, V., 2005. Convective patterns under the Indo Atlantic ‘box’, *Earth planet. Sci. Lett.*, **239**(3-4r), 233-252.

De Viron, O., Dehant, V., Goosse, H., Crucifix, M. & Participing CMIP Model, 2002. Effect of global warming on the length-of-day, *Geophys. Res. Lett.*, **29**, 1146, doi:10.1029/2001GL013672.

Dumberry, M. & Bloxham, J., 2006. Azimuthal flows in the Earth's core and changes in the length of day at millennial timescales, *Geophys. J. Int.*, **165**, 32-46.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary Reference Earth Model PREM, *Phys. Earth Planet. Int.*, **25**, 297-356.

Forte, A.M. & Peltier, W.R., 1991. Viscous-flow models of global geophysical observables. 1. Forward problems, *J. G. R., Solid Earth*, **96**(B12), 20131-20159.

Forte, A.M. & Mitrovica, J.X., 2001. Deep-mantle high-viscosity flow and thermochemical structure inferred from seismic and geodynamic data, *Nature*, **410**(6832), 1049-1056.

Fukao, Y., Widiyantoro, S. & Obayashi, M., 2001. Stagnant slabs in the upper and lower mantle transition region. *Rev. Geophys.*, **39**(3), 291-323.

Greff-Lefftz, M., Pais, A. & Le Mouel, J.L., 2004. Surface geopotential and topography induced by fluid core motions, *J. Geod.*, **78**(6), 386-392. doi:10.1007/s00190-004-0418-x.

Greff-Lefftz M., Métivier, L. & Besse, J., 2010. Dynamic heterogeneities and global geodetic observables, *Geophys. J. Int.*, **180**, 1080-1094.

Gross, R.S., Fukumori, I. & Menemenlis, D., 2005. Atmospheric and oceanic excitation of decadal-scale Earth orientation variations, *J.G.R., Solid Earth*, **110**, B9, doi:10.1029/2004JB003565.

Gurnis, M., Mitrovica, J.X, Ritsema, J. & van Heijst, H.-J., 2000. Constraining mantle density structure using geological evidence of surface uplift rates: The case of the African Superplume, *Geochem. Geophys. Geosyst.*, **1**, 1020, doi:10.1029/1999GC000035.

Jault, D., Gire, C. & Le Mouël, J.L., 1988. Westward drift, core motions and exchanges of angular momentum between core and mantle, *Nature*, **333**, 353-356.

Jault, D. & Le Mouël, J.L., 1991. Exchange of angular momentum between the core and the mantle, *J. Geomag. Geoelectr.*, **43**, 111-129.

Jaupart, C., Labrosse, S. & Mareschal, J.C., 2009, Temperatures, heat and energy in the mantle of the Earth, in *Treatise on Geophysics*, Vol. 7: Mantle dynamics, pp. 253-303, eds Bercovici, D. & Schubert, G., Elsevier, Amsterdam.

Lambeck, K., 1978. The Earth's paleorotation, in *Tidal Friction and the Earth's Rotation*, pp. 145-152, eds Brosche, P. & Sundermann, J. Springer-Verlag, Berlin.

Lambeck, K., 1980. *The Earth's Variable Rotation*. Cambridge University Press, Cambridge, p. 449.

Lefftz, M., Legros, H. & Hinderer, J., 1991. Non-linear equations for the rotation of a visco-elastic planet taking into account the influence of a liquid core, *Celest. Mech. Dyn. Astron.*, **52**, 14-43.

Lithgow-Bertelloni, C. & Richards, M., 1998. The dynamics of cenozoic and mesozoic plate motions, *Rev. Geophys.*, **36**, 27-78.

Lithgow-Bertelloni, C. & Silver, P.G., 1998. Dynamic topography, plate driving forces and the African superswell, *Nature*, **395**(6699), 269-272.

Longman, I.M., 1962. Greens function for determining deformation of Earth under surface mass loads. 1. Theory, *J. geophys. Res.*, **67**(2), doi:10.1029/JZ067i00845.

Machetel, P. & Thomassot, E., 2003. Cretaceous length of day perturbation by mantle avalanche, *EPSL*, **202**(2), 379-386.

Marquart, G., Steinberger, B. & Niehuus, K., 2005. On the effect of a low viscosity asthenosphere on temporal change of the geoid—a challenge for future gravity missions, *J. Geodyn.*, **39**, 493-511.

Meier, M.F., 1984. Contribution of small glaciers to global sea-level, *Science*, **226**, 1418-1421.

Mitzovica, J.X. & Peltier, W.R., 1993. Present-day secular variations in the zonal harmonics of Earths geopotential, *J. G. R.-Solid Earth*, **98**(B3), 4509-4526.

Monnereau, M. & Quere, S., 2001. Spherical shell models of mantle convection with tectonic plates, *E.P.S.L.*, **184**, 575-587.

Munk, W.H. & MacDonald, G.J.F., 1960. *The Rotation of the Earth*, Cambridge University Press, New-York.

Nakiboglu, S.M. & Lambeck, K., 1980. Deglaciation effects upon the rotation of the Earth, *Geophys. J. R. astr. Soc.*, **62**, 49-58.

Peltier, W.R., 1974. Impulse response of a Maxwell Earth, *Rev. Geophys. Space Phys.*, **12**, 649-669.

Ricard, Y., Spada, G., & Sabadini, R., 1993a. Polar Wandering of a dynamic Earth, *Geophys. J. Int.*, **113**, 284-298.

Ricard, Y., Richards, M., Lithgow-Bertelloni, C. & Le Stunff, Y., 1993b. A geodynamic model of mantle density heterogeneity, *J. G. R.*, **B98**(12), 21895-218909.

Richards, M.A., Bunge, H.-P., Ricard, Y. & Baumgardner, J.R., 1999. Polar wandering in mantle convection models, *Geophys. Res. Lett.*, **26**, 1777-1780.

Romanowicz, B. & Gung, Y., 2002. Superplumes from the core-mantle boundary to the base of the lithosphere, *Science*, **296**, 513-516.

Rouby, H., Greff-Lefftz, M. & Besse, J., 2010. Mantle dynamics, geoid, inertia and TPW since 120 Myr, *E.P.S.L.*, **292**(3-4), 301-311, doi:10.1016/j.epsl.2010.01.033

Sabadini, R. & Peltier, W.R., 1981. Pleistocene deglaciation and the Earth's rotation implications for mantle viscosity, *Geophys. J. R. astr. Soc.*, **66**, 553-578.

Spada, G., Yuen, D., Sabadini, R., Morin, P. & Gasperini, P., 1990. A computed-aided, algebraic approach to the post-glacial rebound problem, *Math. J.*, **1**(2), 65-69.

Steinberger, B., 2000. Slabs in the lower mantle - results of dynamic modelling compared with tomographic images and the geoid, *PEPI*, **118**, 241-257.

Tan, E., Gurnis, M. & Han, L.J., 2002. Slabs in the lower mantle and their modulation of plume formation, *Geochem. Geophys. Geosyst.*, **3**, No: 1067.

Trupin, A., Meier, M.F. & Wahr, J.M., 1992. Effect of melting glaciers in the Earth's rotation and gravitational field: 1965-1984, *Geophys. J. Int.*, **108**(1), 1-15.

Varga, P., 2006. Temporal variation of geodynamical properties due to tidal friction, *J. Geodyn.*, **41**(1-3), 140-146.

Wahr, J., Han, D.Z., Trupin, A. & Lindqvist, V., 1993. Secular changes in rotation and gravity—evidence of postglacial rebound or of changes in polar ice, *Adv. Space Res.*, **13**(11), 257-269.
# APPENDIX A: DEGREE 0 VISCOELASTO-GRAVITATIONAL DEFORMATION

In this part, we review the viscoelasto-gravitational theory we use to compute the deformations induced by internal loads varying at geological timescale.

## A1 Equations

This theory is a linear first order theory

(i) At zero order, the planet is assumed to be radially stratified [with a density $ρ_o(r)$, $r$ being the radius] and in hydrostatic equilibrium.

(ii) At first order (perturbed order), the equations governing the elastic deformations within a hydrostatically pre-stressed planet are the momentum equation, the conservation of mass and the Poisson equation (for details see Dahlen & Tromp 1998). A rheological law is necessary to relate the stress to the strain: the Hookean law.

The unknowns are the displacement vector $\vec{u}$ and the mass redistribution potential $\Phi_1$.

In the classical elastogravitational theory, in a frame related to the centre of mass of the Earth, the displacement vector field $\vec{u}$ and the traction $\vec{T}$ are expanded in spherical spheroidal vectors of degree $n$ and order $m$, with the use of six radial function $y_i(r)$ (Alterman et al. 1959).

$$ \vec{u} = \sum_{n=1}^{\infty} \sum_{m=0}^{n} y_{1n}(r) Y_{n}^{m}(\theta, \varphi) \frac{\vec{r}}{r} + r y_{3n}(r) \vec{\nabla} Y_{n}^{m}(\theta, \varphi) \quad (A1) $$

$$ \vec{T} = \sum_{n=1}^{\infty} \sum_{m=0}^{n} y_{2n}(r) Y_{n}^{m}(\theta, \varphi) \frac{\vec{r}}{r} + r y_{4n}(r) \vec{\nabla} Y_{n}^{m}(\theta, \varphi). \quad (A2) $$

We define for the potential (first order mass redistribution potential + exciting internal loading potential $S^{\text{int}}$): $\Phi_1 + S^{\text{int}} = \sum_{n=1}^{\infty} \sum_{m=0}^{n} y_{5n}(r) Y_{n}^{m}(\theta, \varphi)$

and introduce for the radial derivative of the potential (Longman 1962) a function defined by

$$ y_{6n}(r) = \frac{dy_{5n}(r)}{dr} - 4\pi G \rho y_{1n}(r). \quad (A3) $$

For each degree $n$, the elastogravitational system can be then written, as a first order differential equations in a form given by Alterman et al. (1959)

$$ \frac{dy_{1n}(r)}{dr} = A_{ij}^{n} y_{jn}(r), \quad (A4) $$

where $A_{ij}^{n}$ is a 6 × 6 matrix whose elements are function of the compressibility $K(r)$, the rigidity $\mu(r)$, the density $\rho_o(r)$ and the gravity $g_o(r)$. We assume that the deformations are static, that is to say we neglect the inertial acceleration in the momentum equation.

For the degree 0, the deformations do not dependent on the latitude and longitude and consequently the displacement is radial ($y_3(r) = 0$) and there is no tangential stress [$y_4(r) = 0$].

## A2 Boundary conditions

To solve this differential system, we have to add boundary conditions which depend on the excitation sources (Greff-Lefftz et al. 2010). Here, the source is a mass anomaly $\delta M$, within the mantle, located at the radius $r_o$. This mass anomaly induces an internal loading potential $S^{\text{int}}$ such as, for the degree 0

$$ S^{\text{int}} = S_o \begin{cases} \frac{r_o}{r} & \text{if } r_o \le r \\ 1 & \text{if } r_o \ge r \end{cases} \quad \text{with} \quad S_o = \frac{G \delta M}{r_o}. \quad (A5) $$

The unknowns radial functions $y_1(r)$, $y_2(r)$, $y_5(r)$ and $y_6(r)$, are continuous within the Earth excepted at the radius $r_o$ where is located the load: at this interface, the continuity of the radial stress and that of the gravity will depend on this excitation source

$$ y_2(r_o^+) - y_2(r_o^-) = \frac{1}{3} \frac{g_o(r_o)}{g_o(a)} \frac{a}{r_o} \tilde{\rho} S_o; \quad y_6(r_o^+) - y_6(r_o^-) = -\frac{a}{r_o} \frac{S_o}{a}, \quad (A6) $$

$\tilde{\rho}$ is the mean density of the Earth.

At the the Earth surface, we have

$$ y_2(a) = -\frac{1}{3} \tilde{\rho} S_o(a); \quad \text{and} \quad y_6(a) + \frac{y_5(a)}{a} = \frac{S_o(a)}{a}. \quad (A7) $$

We solve the system (A4) with these boundary conditions: finally, we know, in the whole Earth, the radial displacement $y_1(r)$ and mass redistribution potential $\Phi_1(r)$ induced by a mass anomaly within the mantle.
A3 Moment of inertia induced by the degree 0 mantle mass anomaly

Let us first assume that the planet is rigid and then compute the direct effect of such internal mass, located at the radius $r_o$, on the perturbation of the moment of inertia $\delta I_o^R$.

$$ \delta I_o^R = \frac{2}{3} \delta M r_o^2. \qquad (A8) $$

Because of the internal loading, there are viscoelasto-gravitational deformations within the planet. The degree 0 internal load will induce a radial displacement $\vec{u} = y_1(r)\vec{e}_r$ and a redistribution of mass leading to a perturbation of the density $\rho_1 = -\text{div}(\rho_o\vec{u})$, that ensures the conservation of the mass.

The moment of inertia of a radially stratified planet may be computed from

$$ I_o + \delta I_o = \frac{8\pi}{3} \int_0^{a+y_1(a)} (\rho_o + \rho_1)r^4 dr + \frac{2}{3}\delta M r_o^2 \qquad (A9) $$

which may be written in a first order theory

$$ I_o + \delta I_o = \frac{8\pi}{3} \int_0^a \rho_o r^4 dr + \frac{8\pi}{3} \int_a^{a+y_1(a)} \rho_o r^4 dr + \frac{8\pi}{3} \int_0^a \rho_1 r^4 dr + \frac{2}{3} \delta M r_o^2. \qquad (A10) $$

The first term of the right member is $I_o$ the moment of inertia of our reference sphere. The three last terms are the perturbation of the moment of inertia $\delta I_o$ which is composed with $\delta I_o^R = \frac{2}{3}\delta M r_o^2$, the direct (rigid) effect induced by the internal load, and with the second and third terms related to the surface deformation and the mass redistribution. In spherical coordinates, we have $\text{div}(\rho_o\vec{u}) = \frac{1}{r^2}\frac{\partial}{\partial r}[\rho_o r^2 y_1(r)]$ and consequently, $\delta I_o$ may be simply written

$$ \delta I_o = \frac{2}{3} \delta M r_o^2 + \frac{16\pi}{3} \int_0^a \rho_o(r) r^3 y_1(r) dr. \qquad (A11) $$

Using the solution $y_1(r)$ of the elasto-gravitational differential system, $\delta I_o$ may be computed. We introduce an 'inertia' degree 0 loading Love number $k'_0(r_o)$ such as

$$ \delta I_o = \frac{2}{3} \delta M a^2 \left[ \frac{r_o^2}{a^2} + k'_0(r_o) \right]. $$

The kernel $[r_o^2 + k'_0(r_o)]$ is plotted in Fig. A1, for the PREM model (Dziewonski & Anderson 1981), as a function of $r_o$ in solid line. We plot also this kernel for an incompressible planet ($K > \infty$ in PREM) in dotted line, and for a fluid planet ($\mu > 0$ in PREM) in dashed line.

Note that if the planet is incompressible or rigid, this loading Love number $k'_o$ will be exactly zero. The fluid and elastic compressible cases are very similar. It means that for the degree 0, the deformation is essentially dependent on the compressibility. The rigidity is related to the shear which does not appear in radial deformations.

Now, we want to extend these results to a viscoelastic mantle. In the Laplace domain, this stress-strain relation is the Hookean law but the rigidity $\mu^v$ and the compressibility $K^v$ are now function of the frequency $s$

$$ \mu^v(s) = \mu \frac{s}{s + \frac{\mu}{v}} \quad \text{and} \quad K^v(s) = K. \qquad (A12) $$

In the fluid limit ($s > 0$), the viscoelastic rigidity tends to zero, that is to say the planet has the behaviour of an inviscid fluid, whereas the compressibility remains constant and equal to the elastic one.

<figure><img src="image_15.png" /><figcaption>Figure A1. Degree zero kernel $[r_o^2 + k'_0(r_o)]$ as a function of $r_o$ for an incompressible PREM-Earth (dotted line), an elastic compressible Earth (solid line) and a fluid compressible Earth (dashed line), for the PREM model.</figcaption></figure>
In Laplace domain, the viscoelastic equations and the boundary conditions are the same than those for an elastic body with the same geometry. Consequently, we may use the correspondence principle: we solve the elastic problem for different frequencies in order to build the viscoelastic solutions (Peltier 1974). We find that the Love numbers $k_0'(r_o)$ is not significantly dependent on the frequency. This is due to the fact that, for the elastic PREM model, the compressibility is always stronger than the rigidity and that for radial deformation, the rigidity (the shear modulus) is negligible. This result was already appeared in the previous Fig. A1 where we saw that the elastic and fluid kernel are quasi-identical. Consequently, to compute the perturbation of the moment of inertia induced by degree 0 internal load, we do not need to investigate the viscoelastic relaxation modes but we can simply put $\mu = 0$ in our elasto-gravitational system.