# Dynamic mantle density heterogeneities and global geodetic observables

Marianne Greff-Lefftz, Laurent Métivier, Jean Besse

► To cite this version:

Marianne Greff-Lefftz, Laurent Métivier, Jean Besse. Dynamic mantle density heterogeneities and global geodetic observables. Geophysical Journal International, 2010, 180, pp.1080-1094. 10.1111/j.1365-246X.2009.04490.x. insu-03605340

HAL Id: insu-03605340

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamic mantle density heterogeneities and global geodetic observables

Marianne Greff-Lefftz,<fnref n="1" /> Laurent Métivier<fnref n="1" />,<fnref n="2" /> and Jean Besse<fnref n="1" />

Accepted 2009 December 16. Received 2009 October 30; in original form 2009 July 9

## SUMMARY

We investigate the influence of mantle dynamics on low degree deformations of the Earth at geological timescale. We first compute surface deformations, and discuss the analytical form of the tangential surface displacement induced by internal loads, in a reference frame related to the centre of mass of the planet. We use the theoretical Love numbers formalism since the Earth has a viscoelastic behaviour at geological timescale.

Then we quantify degree-one and degree-two deformations induced by upwelling domes and subducted plates sinking into the mantle. We use a simple model in which the slabs are modelled as blobs diving vertically through the mantle, and in which the domes are assumed to be stable over the last 120 Ma. Their location is modelled from seismic tomography within the lower mantle. The temporal evolutions of the $J_2$ gravitational potential coefficient and of the geocentre motion are plotted since 120 Ma.

We find that:

(1) The mantle density heterogeneities within the mantle can explain the present-day non-hydrostatic flattening of the Earth. However they vary at a too slow timescale to significantly perturb the $J_2$ coefficient.

(2) Although there is a significant discrepancy of about a few hundred metres between the centre of figure and the centre of mass of the Earth, the secular variation of the geocentre motion is one order of magnitude smaller than the one induced by surface loads.

**Key words:** Time variable gravity; Global change from geodesy; Dynamics of lithosphere and mantle.

## 1 INTRODUCTION

The recent space geodesy measurements [Doppler Orbitography of Radiopositioning Integrated by Satellite (DORIS), Global Positioning System (GPS) and Satellite Laser Ranging (SLR)] give the variations of centre of figure of the outer surface of the Earth with respect to a fixed reference frame, the so-called geocentre (Blewith 2003). These observations (since 1993) have annual and semi-annual components with amplitudes ~10 mm. Since very recent time, one observed values of the time-derivative of the geocentre coordinates. The up-to-date realization of the International Terrestrial Reference System (ITRF2005) presents a particularly large translation rate of 1.8 mm yr⁻<fnref n="1" /> on the Z-component with respect to the ITRF2000. It is not sure that such rate can be interpreted as a secular motion of the geocentre: as a matter of fact, it may probably be linked to the inhomogeneous shape of the station measurement network and other technique systematic errors (Collilieux *et al.* 2009). It nevertheless allows to quantify the order of magnitude of the time-derivative of the geocentre coordinates: less than a few millimetre per year.

The geocentre motion, defined as the figure centre of the deformed outer surface, is related to spherical harmonics degree 1 deformations of the Earth. Some geophysical studies have computed the displacement vector between the centre of figure of the outer deformed surface of the Earth and the centre of mass induced by the atmospheric loading and the ground water storage (Dong *et al.* 1997; Greff-Lefftz & Legros 1997; Chen *et al.* 1999) in order to explain its annual and semi-annual variations. At decadal timescale, the magnetic pressure within the fluid core acting at the CMB deforms the mantle and involves weak perturbation of the geocentre motion, at a level of 0.1 mm yr⁻<fnref n="1" /> (Greff-Lefftz & Legros 2007). At secular timescale, the geocentre motion induced by postglacial rebound has been shown to be at the level of -0.4 to 0.2 mm yr⁻<fnref n="1" /> (Greff-Lefftz 2000). A recent study (Métivier *et al.* 2009) found that the recent climate changes and the associated ice melting and sea level rise can contribute, essentially in the Z coordinates, to a rate of secular geocentre motion of about a few tenths of a millimetres per year. Geocentre motions are due to the combined effect of mass re-distributions within the planet and changes in the Earth's shape. For

<footnote n="1">¹Institut de Physique du Globe de Paris 4 place Jussieu, 75252 Paris 05, France. E-mail: greff@ipgp.jussieu.fr</footnote>

<footnote n="2">²Institut Géographique National, LAREG/GRGS, Champs-sur-Marne, France</footnote>
these reasons, at secular and geological timescales, the geocentre position should be affected by the mantle dynamics, which generate lateral variations of density within the Earth, and motions of the surface plates and continents, which change the Earth's shape. These contributions on the geocentre motion are today poorly known.

According to space geodetic observations over the past 25 yr, the zonal degree 2 coefficient of the Earth's gravitational potential $J_2$-related to the Earth's dynamic oblateness, has undergone a secular decrease with superimposed decadal variations (Cox & Chao 2002) and, with the well-known superimposed seasonal variation induced by the mass-redistribution within the oceans and the atmosphere. Several mechanisms which could explain the decadal observations have been discussed: the melting of glaciers or polar ice caps (Cox & Chao 2002; Dickey et al. 2002), mass moving in the fluid outer core (Greff-Lefftz et al. 2004) or in oceans (Cazenave & Nerem 2002).

The secular decrease is usually associated with the post glacial rebound and is used in a non-linear inverse problem approach to give information about the internal viscoelastic structure of the planet (Peltier 1985; Mitrovica & Peltier 1991, 1993; Bills & James 1997; Vermeersen et al. 1998) and also on the history of the deglaciation (Tushingham & Peltier 1991). These approaches assumed that the secular variation of the gravitational field is only induced by the last deglaciation. Forte & Mitrovica (1996) performed joint inversions for mantle viscosity of geophysical observables associated with both mantle convection and glacial isostatic adjustment (GIA) and found that both data sets may be reconciled using a single profile of mantle viscosity.

The geological evolution of the Earth's rotational axis is most likely controlled by internal mass redistribution within the mantle (e.g. Richards et al. 1997; Steinberger & O'Connell 1997, 2002). Palaeomagnetic evidence from recent Earth history (Besse & Courtillot 2002) suggests that this motion has been small with a rate about 0.1°–0.4° per Myrs, that is to say about a few cm yr$^{-1}$. Consequently, changes in mantle density heterogeneities induce a polar wander at geological timescale with the same order of magnitude than the one induced by the last deglaciation at secular timescale.

The time-variable mantle density heterogeneities also affect global geodetic observables such as the time derivative of the degree 2 coefficient of the gravitational field or the geocentre motion. How mantle density heterogeneities do affect the present-day secular variation of the Earth's geocentre and $J_2$ coefficient is the aim of this paper.

In Section 2, we review the viscoelasto-gravitational theory and we present the analytical form of the Love numbers related to surface displacement induced by internal loads. In Section 3, we compute the present-day geoid, the $J_2$ coefficient and its time-derivative, the present-day surface topography and the geocentre variation for a simple model of time-dependent mass anomalies within the mantle. The values of these parameters and their time-derivatives produce independent pieces of information which in some cases allow to separate the contributions induced by surface loads and the ones induced by internal loads. The effects of plate tectonics are also discussed. The conclusions are presented in Section 4.

## 2 VISCOELASTIC LOVE NUMBERS AND INTERNAL LOADS

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

### 2.1 Equations

The equations governing the elastic deformations within a hydrostatically pre-stressed planet are the momentum equation, the conservation of mass and the Poisson equation. A rheological law is necessary to relate the stress to the strain. 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 vector 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 (1) $$

$$ \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 (2) $$

where r is the radius. We note for the potential $U = \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 (3) $$

where $Y_n^m$ are the spherical harmonics, $\theta$ the colatitude and $\varphi$ the longitude. The spherical coordinates are defined in a geographical reference frame centred on the centre of mass of the planet and in which the x (respectively the y) axis is the intersection between the equatorial plane and the Greenwich meridian (respectively the 90°E meridian) and the z-axis is the axis perpendicular to the equatorial plane.

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

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

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

Our earth model consists of N homogeneous incompressible layers: a lithosphere (layer 1), a stratified mantle (layer 2 to $N-2$), an inviscid fluid core (layer $N-1$) and an inner core (layer $N$).

The analytical solutions of the elastogravitational system are described in the Appendix for such an earth model, when the excitation source is a surface load or an internal load.

### 2.2 Viscoelastic mantle and lithosphere

At secular and geological timescales, we have to take into account the viscous response of the Earth. To do that we have chosen as rheological law, a linear viscoelastic Maxwell model of rheology. In Laplace domain, the stress–strain relation for an incompressible Maxwell body is the Hookean law, but the rigidity is function of the frequency s:

$$ \mu^i(s) = \mu_{el}^i \frac{s}{s + \frac{\mu_{el}^i}{v^i}}, $$

where $v^i$ is the viscosity of the layer i and $\mu_{el}^i$ its elastic shear modulus.

In Laplace domain, the viscoelastic equations and the boundary conditions are the same than those for an elastic body with
<table><caption>**Table 1.** Geometrical and physical parameters of our five-layered earth model (PREM-averaged values).</caption><thead><tr><th>Layer (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; 6371</td><td>4414</td><td>1.66 &times; 10<sup>11</sup></td><td>11 &times; 10<sup>21</sup></td></tr><tr><td>5701 &lt; r &lt; 6271</td><td>4414</td><td>1.66 &times; 10<sup>11</sup></td><td>1 &times; 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>1225.5 &lt; r &lt; 3480</td><td>12420</td><td>0</td><td>0</td></tr><tr><td>0 &lt; r &lt; 1225.5</td><td>12420</td><td>1.64 &times; 10<sup>11</sup></td><td>1 &times; 10<sup>13</sup></td></tr></tbody></table>

Note: VR is the viscosity ratio between lower and upper mantle.

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 compute the determinant of the viscoelasto-gravitational system of boundary conditions, and we find a frequency dependent polynomial whose zeros are called relaxation modes. For the deformation of degree $n > 1$, it is well known that these relaxation modes are generated at each interface by the discontinuity of physical parameters. For example, for a five 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$), an in- viscid fluid core (with a density $\rho^c$) and a viscoelastic inner core, we have seven relaxation modes (Peltier 1974): $M_o$ due to the surface discontinuity at $r = a$, $C$ and $G$ because of the density discontinuity between the fluid core and the lower mantle and between the fluid core and the inner core, $L_1$ and $L_2$ due to the discontinuity in the Maxwell time between the lithosphere and the upper mantle, $T_1$ and $T_2$ due to a jump in the Maxwell time between the upper and the lower mantle and called transition modes because they relax rapidly and are weakly excited.

For the degree one, the initial system of boundary conditions degenerates, and consequently the number of modes is different. We find that the $M_o$ mode disappears and that there is only one transition mode instead of two at each interface with a viscosity or rigidity jump (Greff-Lefftz & Legros 1997).

The geometrical and physical parameters of our five-layered earth model are given in the Table 1 where VR is the viscosity ratio between lower and upper mantle. As done by Lithgow-Bertelloni & Richards (1998), we introduce an effective viscosity for the lithosphere of about $11 \times 10^{21}$ Pa s in order to take into account plate tectonics.

We plot, on Fig. 1, the viscoelastic relaxation times for the degree $n = 1$ (left-hand side) and the degree $n = 2$ (right-hand side) in Kyr as a function of the viscosity ratio VR ranging from 5 to 100.

Note that the transition times (related to the viscosity contrast between upper mantle and lithosphere and to viscosity contrast between lower and upper mantle) do not depend on VR, for $VR > 10$, whereas the relaxation time associated with the mode C, for $n = 1$, or with the modes $M_o$ and $C$, for $n = 2$, strongly increases with VR. The order of magnitude of these relaxation times are smaller than 100 Kyr, whatever VR, that is to say are small in comparison with the timescale of the mantle convection (about the million years).

## 2.3 Viscoelastic Love numbers

The solutions of the degree $n$ viscoelasto-gravitational deformations are usually written using viscoelastic degree $n$ Love numbers, noted $k'_n(s, r_o)$ for the mass redistribution potential, $h'_n(s, r_o)$ for the radial displacement or $l'_n(s, r_o)$ for the tangential displacement. For surface ($S^e$) or internal ($S^{\text{int}}$) loading at radius $r_o$, the gravitational potential and the radial and tangential displacements may be written, at the Earth's surface ($r = a$).

(i) In the Laplace transform domain

$$y_5(a, s) = [1 + k'_n(s)]S^e(s) + \left[ \left(\frac{r_o}{a}\right)^{n+1} + k'_n(s, r_o) \right] S^{\text{int}}(s, r_o) \quad (5)$$

$$y_1(a, s) = h'_n(s) \frac{S^e(s)}{g_o} + h'_n(s, r_o) \frac{S^{\text{int}}(s, r_o)}{g_o} \quad (6)$$

$$y_3(a, s) = l'_n(s) \frac{S^e(s)}{g_o} + l'_n(s, r_o) \frac{S^{\text{int}}(s, r_o)}{g_o}. \quad (7)$$

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Viscoelastic relaxation times for the degree $n = 1$ (left-hand panel) and the degree $n = 2$ (right-hand panel) in Kyr as a function of a viscosity ratio VR ranging from 5 up to 100. Circle symbol: transition times associated with the viscosity jump between lower and upper mantle. Square symbol: transition times associated with the viscosity jump between upper mantle and lithosphere. Triangle up symbol: time associated with C mode induced by density jump at the CMB. Star: time associated with $M_o$ mode induced by density jump at the surface.</figcaption></figure>
(ii) In the temporal domain

$$y_5(a, t) = [\delta(t) + k'_n(t)] * S^e(t) + \left[ \left(\frac{r_o}{a}\right)^{n+1} \delta(t) + k'_n(t, r_o) \right] * S^{\text{int}}(t, r_o) \quad (8)$$

$$y_1(a, t) = h'_n(t) * \frac{S^e(t)}{g_o} + h'_n(t, r_o) * \frac{S^{\text{int}}(t, r_o)}{g_o} \quad (9)$$

$$y_3(a, s) = l'_n(t) * \frac{S^e(t)}{g_o} + l'_n(t, r_o) * \frac{S^{\text{int}}(t, r_o)}{g_o}, \quad (10)$$

where * denotes the temporal convolution, that-is-to-say the viscous memory of the planet, and $\delta(t)$ is the Dirac function, which expresses the instantaneous deformation induced by the source. $g_o$ is the surface gravity.

### 2.3.1 Surface load viscoelastic Love numbers

These Love numbers have already been investigated for surface loading. For a Maxwell model of rheology, Peltier (1974) (see also Spada et al. 1990) has shown that the viscoelastic Love numbers have the following form, whatever the degree $n$, where the subscript $n$ has been omitted for simplicity

$$h'(s) = h^e + \sum_{j=1}^{M} \frac{h_j \tau_j}{1 + s \tau_j} \quad \text{or} \quad h'(t) = h^e \delta(t) + \sum_{j=1}^{M} h_j e^{-t/\tau_j} H(t)$$

$$l'(s) = l^e + \sum_{j=1}^{M} \frac{l_j \tau_j}{1 + s \tau_j} \quad \text{or} \quad l'(t) = l^e \delta(t) + \sum_{j=1}^{M} l_j e^{-t/\tau_j} H(t)$$

$$k'(s) = k^e + \sum_{j=1}^{M} \frac{k_j \tau_j}{1 + s \tau_j} \quad \text{or} \quad k'(t) = k^e \delta(t) + \sum_{j=1}^{M} k_j e^{-t/\tau_j} H(t),$$

where $H(t)$ is the Heaviside distribution.

The first term in the right member $h^e$ is the instantaneous elastic Love number, $\tau_j$ and $h_j$ are the relaxation times and residues of the $M$ modes of the earth model. $h_j$ depend on the excitation sources.

For the degree one, the form is conserved and only the number of relaxation modes changes.

### 2.3.2 Internal load viscoelastic Love numbers

#### Analytical form

For internal load, we have analytically computed the frequential form of the Love numbers, solving the viscoelasto-gravitational set of boundary conditions described in the Appendix, using the Maple software package.

Note that in the fluid limit (for a Maxwell model of rheology, the fluid limit $s = 0$ corresponds to an inviscid fluid mantle), an internal load in an homogeneous mantle is not stable (because there is no more viscous stress, see boundary condition in $y_2$ at the radius $r_o$) and that the linear theory of the viscoelasto-gravitational deformation becomes not valid. Nevertheless, we can investigate the viscous fluid limit, that is to say the deformation remaining after viscoelastic relaxation, which is equivalent to a Newtonian viscous limit.

We find that, in the Laplace domain, the form of the viscoelastic Love numbers associated with an internal load located at the

radius $r_o$ is

$$h'(s, r_o) = h^e(r_o) + \sum_{j=1}^{M} \frac{h_j(r_o) \tau_j}{1 + s \tau_j}$$

$$l'(s, r_o) = l^e(r_o) + \frac{l_0(r_o)}{s} + \sum_{j=1}^{M} \frac{l_j(r_o) \tau_j}{1 + s \tau_j}$$

$$k'(s, r_o) = k^e(r_o) + \sum_{j=1}^{M} \frac{k_j(r_o) \tau_j}{1 + s \tau_j}.$$

Note that there is an additional term in $1/s$ for the viscoelastic tangential surface displacement, $l_0(r_o)$, which vanishes for loads located at the core-mantle boundary (CMB) and at the surface. This term means that if the mantle is an inviscid fluid, the tangential displacement and stress induced by loads within the mantle are not defined. This term does not appear for the radial displacement: in the fluid limit, the radial problem is well defined and we obtain that the form is the equipotential within the mantle. The term in $1/s$ does not also appear in the $k'(s, r_o)$ Love number, because this last one is related to the mass redistribution potential and consequently for an incompressible Earth depends only on the radial displacements at interfaces with density jumps (in our model at the CMB and at the surface).

If the mantle is assumed to be rigid, the viscoelastic Love numbers will vanish excepted for the degree $n = 1$, for which we have

(i) $k'_1(s, r_o) = -(r_o/a)^2$ because of the conservation of the centre of mass ($[y_5(a, s) = 0]$)

(ii) $h'_1(s, r_o) = l'_1(s, r_o) = -(r_o/a)^2$, because the rigid translation is a solution of the viscoelasto-gravitational set of eqs (A2) for the degree $n = 1$ (see the term in $C_4^i$ in the propagators of $y_1(r)$ and $y_3(r)$ in the Appendix).

At the timescale of the mantle convection we have $s\tau_j \ll 1$ and consequently we can introduce a quasi-fluid approximation for the Love numbers

$$h'(s, r_o) \approx h^e(r_o) + \sum_{j=1}^{M} h_j(r_o) \tau_j$$

$$l'(s, r_o) \approx l^e(r_o) + \frac{l_0(r_o)}{s} + \sum_{j=1}^{M} l_j(r_o) \tau_j$$

$$k'(s, r_o) \approx k^e(r_o) + \sum_{j=1}^{M} k_j(r_o) \tau_j.$$

These frequential forms are valid, whatever the degree $n$ of the Love numbers.

#### Degree one and two Love numbers

In this part, we investigate in details the degree 1 and degree 2 viscoelastic internal loading Love numbers.

We plot on the Fig. 2 the three kernels for the surface radial displacement as a function of the radius of the load $r_o$, for the degree $n = 1$ (left-hand side) and the degree $n = 2$ (right-hand side), for our five-layers earth's model.

(i) Rigid Love number $h_1^R(r_o) = -(r_o/a)^2$ for $n = 1$ and $h_2^R(r_o) = 0$ for $n = 2$.

(ii) Elastic Love number $h^e(r_o)$.

(iii) Fluid Love number $h^e(r_o) + \sum_{j=1}^{M} h_j(r_o) \tau_j$.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Kernels for the surface radial displacement, for the degree $n = 1$ (left-hand panel) and the degree $n = 2$ (right-hand panel): elastic Love number (solid line) and fluid Love number (dashed line).</figcaption></figure>

Note that for the degree 1, the rigid and the elastic radial displacement are quasi identical. It can be explained by looking at the analytical form of the discrepancy between $h^t(r_o)$ and $h_1^R(r_o)$, which is proportional to the density contrast between mantle and fluid core and is very small with respect to $h_1^R(r_o)$. The rigid translation is the dominant term in the displacement of the planet for $n = 1$.

The fluid Love numbers differ from the elastic ones essentially when the load is located at the surface or at the CMB where it is isostatically compensated.

The elastic and fluid Love numbers are negative whatever the depth of the internal load. It means that a positive mass anomaly within the mantle will induce a negative surface topography.

We plot now, on Fig. 3, the four kernels for the surface tangential displacement, for the degree $n = 1$ (left-hand side) and the degree $n = 2$ (right-hand side): $l^R(r_o)$, $l^e(r_o)$, $l^e(r_o) + \sum_{j=1}^M l_j(r_o)\tau_j$ and $l_0(r_o)$ in Myr$^{-1}$, as a function of the radius of the load $r_o$, assuming a viscosity ratio $VR = 40$.

Similarly to radial displacements, rigid and elastic degree-one tangential displacements are quasi identical, because the global displacement of the planet is dominated by a rigid translation. The dimension of $l_0(r_o)$ is $[T]^{-1}$, where $T$ is time dimension. For a characteristic time of mantle dynamics $t \gg 25\,000$ yr, we have $l_0(r_o) t \gg 1$, that is to say very large in comparison with the rigid or elastic Love numbers. This is the limit of our approach. We cannot compute tangential displacement using our viscoelasto-gravitational approach

for timescale larger than 25 000 yr, because the shear modulus vanishes and consequently the tangential stress and displacement within the mantle are not defined. We have tested the influence of VR on this 'limit': for plausible value of VR ranging from 10 to 100, the maximum value of $l_0(r_o) t$ (reached at the upper-lower mantle discontinuity) is between 35 and 50. Consequently, our approach is valid for characteristic times smaller than 20 000 yr for low VR, or 29 000 yr for high VR.

### 2.3.3 Degree n viscous geoid and topography kernels

The Love number $k'(s, r_o)$ is relative to the mass redistribution potential. The surface degree $n$ geoid induced by a degree $n$ internal load may be written as a sum of the direct gravitational effect of the load + the viscoelastic mass redistribution potential

$$ \frac{y_5(a, s)}{g_o} = \left[ \left( \frac{r_o}{a} \right)^{n+1} + k'(s, r_o) \right] \frac{S^{\text{int}}(s, r_o)}{g_o} $$

We define the viscous geoid kernel as the quantity $[k^e(r_o) + \sum_{j=1}^M k'_j(r_o)\tau_j + (\frac{r_o}{a})^{n+1}]$, function of $r_o$ the radius where is located the load. This kernel is a good approximation when the time scale of internal loads greatly exceeds the viscoelastic relaxation times ($t \gg \tau_j$).

We plot on the Fig. 4 (left-hand panel) the geoid kernel for various degree $n = 2, 4, 8$ and $12$. For the degree $n = 2$, 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 its sign at the depth 1200 km. The depth in which the geoid vanishes depends on the degree $n$. The larger is $n$, the larger is the depth: for $n = 12$ this depth is 2340 km.

These kernels are close to those computed in studies based on the best fitting between the observed geoid and the one induced by mantle density heterogeneities derived from tomography and/or geodynamic models (e.g. Steinberger 2000; Marquart et al. 2205). These studies assumed a Newtonian viscous stratified mantle.

On the right-hand panel of the Fig. 4 is plotted the surface topography kernel $h^t(r_o) + \sum_{j=1}^M h_j(r_o)\tau_j$ for various degree $n = 2, 4, 8$ and $12$. This kernel is always negative such as a positive mass anomaly in the mantle will induce a negative topography. Note that the surface topography is essentially induced by mantle density in the upper mantle.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Kernels for the surface tangential displacement, for the degree $n = 1$ (left-hand panel) and the degree $n = 2$ (right-hand panel): elastic Love number (solid line), fluid Love number (dashed line) and $l_0(r_o)$ in Myr$^{-1}$ in dot-dashed line. The y-axis is $r_o/a$ the normalized radius of the load.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Geoid kernel and surface topography kernel for degree $n = 1$ (grey line), $n = 2$ (solid line), $n = 4$ (dotted line), $n = 8$ (dashed line) and $n = 12$ (dot-dashed line). The y-axis is $r_o$ the radius of the load.</figcaption></figure>

For a surface load, we can analytically show that the fluid Love number is

$$h'^e(a) + \sum_{j=1}^{M} h_j(a) \tau_j = - \frac{2n+1}{3} \frac{\tilde{\rho}}{\rho(a)}$$

which means that a surface mass density causes a large surface displacement such as the mass is isostatically compensated. $\tilde{\rho}$ is the mean density of the Earth as defined in the Appendix.

### 2.3.4 Viscoelastic Love numbers in temporal domain

To end this theoretical part, we write the viscoelastic Love numbers in the temporal domain

$$h'(t, r_o) = h'^e(r_o)\delta(t) + \sum_{j=1}^{M} h_j(r_o)e^{-t/\tau_j} H(t)$$

$$k'(t, r_o) = k'^e(r_o)\delta(t) + \sum_{j=1}^{M} k_j(r_o)e^{-t/\tau_j} H(t)$$

$$l'(t, r_o) = l'^e(r_o)\delta(t) + l_0(r_o)H(t) + \sum_{j=1}^{M} l_j(r_o)e^{-t/\tau_j} H(t).$$

## 3 GEOPHYSICAL APPLICATION

In this part, we compute the surface gravitational potential (geoid and $J_2$ coefficient), the surface topography and the geocentre motion for a simple model of time-dependent mass anomalies within the mantle.

### 3.1 Temporal evolution of the mantle mass anomalies since 120 Ma

We first present our simple model for the large-scale pattern of mantle dynamics. This model is described in details in Rouby *et al.* (2010). From mantle density heterogeneities, the model may explain the observed wander of the rotation axis since 120 Ma. In this section, we just point out its main features.

#### 3.1.1 Upwellings

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 with very large periodicities ranging from some hundred Myr up to 1 Gyr. Such domes could be responsible for the super-swells observed at the Earth's surface. 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 the subduction. Consequently, we have assumed that the mass anomalies in the lower mantle associated with these large-scale upwellings are constant with time since 120 Myr, and we have modelled their spatial distribution from the present-day tomography within the lower mantle, assuming a constant density contrast $\Delta\rho = -50 \text{ kg m}^{-3}$ between the 'hot' domes and the surrounding mantle.

#### 3.1.2 Downwellings

Subducted lithosphere is a major component of mantle density heterogeneity. We start with the model of mantle density heterogeneity derived by Ricard *et al.* (1993). This model uses plate-motion reconstruction under the assumption that subducted slabs sink vertically into the mantle (Fig. 5). When the slab is in the upper mantle, its diving velocity is its surface velocity. How such plates cross the 670 km discontinuity is not well known. In the model proposed by Ricard *et al.* (1993) that we used, the velocity of the plates is assumed to decrease by a factor $F$ whereas the size (the thickness) of the plates increases by the same factor $F$. They suggest that such a factor may be written as the log of the viscosity ratio between lower and upper mantle, that is to say for $VR$ ranging from 10 up to 100 may vary from 2.5 up to 4.5. For $VR = 40$, the factor is about $F = 4$. The density contrast between the cold plates and the surrounding mantle is assumed to remain constant with depth and equal to 80 kg m$^{-3}$. The plate-motion reconstruction and details about the slabs modelling are reviewed in (Lithgow-Bertelloni & Richards 1998).

Lithgow-Bertelloni & Richards (1998) modelled the temporal evolution of the mantle density heterogeneities associated with large-scale pattern of plate tectonic motions since 120 Ma, assuming that the plates sink until the CMB. However, the tomographic images, inferred from seismology and from mineral physics, show that
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Cartoon representing the construction of a density heterogeneity model for the Earth based on subduction history from Lithgow-Bertelloni & Richards (1998).</figcaption></figure>

the subducted plates do not reach systematically the CMB (Fukao et al. 2001; Ricard et al. 2005). We have thus also investigated the influence of the depth down to which the plates present a significant density contrast with respect to the surrounding mantle (Rouby et al. 2010).

### 3.1.3 Our simple model

Our 'preferred' model of mantle density heterogeneities is the model for which we obtain:

(1) A present-day heterogeneous structure of the mantle in good agreement with results obtained from seismic tomography,

(2) The best variance reduction between the present-day computed and observed geoid (for degree 2 up to 12).

For each test, the mantle mass anomalies have been expanded in spherical harmonics, at each radius r:

$$ \Delta\rho(r, \theta, \varphi, t) = \sum_{n=0}^{\infty} \sum_{m=0}^{n} [\Delta\rho_{n}^{cm}(r, t) \cos m\varphi + \Delta\rho_{n}^{sm}(r, t) \sin m\varphi] P_{n}^{m}[\cos(\theta)], $$

where $P_n^m$ are the Legendre functions, $\Delta\rho_n^cm$ and $\Delta\rho_n^sm$ are the coefficients of the spherical harmonics expansion, normalized to $4\pi$, in kg m⁻³.

Using the geoid kernel (Fig. 4), we have computed the associated present-day geoid (for degree 2 up to 12) and its variance reduction with respect to the observed geoid.

Our final model assumes time-independent large-scale upwellings within the lower mantle. All the plates reach the CMB with exception of North America sinking down to 2250 km and South America that stops its diving at a depth about 800 km. Topping its diving means that a plate does not display a significant density contrast with respect to the surrounding mantle.

The associated present-day geoid (for degree 2 up to 12) is compared to the observed one in Fig. 6. For this model, we obtain a variance reduction of $\Phi = 0.91$ for the degree 2 alone, and $\Phi = 0.79$ for degree 2–12. Cancelling in the deep mantle the density contrast of some plates with respect to the surrounding mantle allows to ameliorate significantly the present-day computed geoid in comparison with the case where all the plates sink until the CMB. We found that the degree 2 component of the geoid is now in good agreement with the observed one and the maximum principal inertia axis is close to the North Pole.

### 3.1.4 Surface dynamic topography

We then compute the present-day radial surface displacement using the surface topography kernel plotted in Fig. 4. We plot on Fig. 7 the dynamic topography associated with our model. For degree $n = 1, \dots, 12$, a slab sinking into the mantle will involve a negative topography whereas a hot dome within the lower mantle will involve a surface swell. The maximum negative anomaly of the topography is around the Pacific, below large subduction zones. Note that, in our

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Geoid in metre, from degree 2 to 12. Left-hand panel: observed geoid from EIGEN-GEODERMA (Joerste et al. 2008) in which the hydrostatic reference ellipsoid (Nakiboglu 1979, 1982) and the effect of postglacial rebound have been corrected. The gravitational contribution of the postglacial rebound has been corrected using the simple model (three ice-sheets with the same dependence) detailed in Greff-Lefftz (2000) and our surface loading viscoelastic Love numbers for VR = 40. Right-hand panel: present-day computed geoid induced by our mantle density heterogeneities model, for a mantle with a viscosity ratio VR = 40, using geoid kernels plotted in Fig. 4.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Present-day surface topography (in meter)</figcaption></figure>

**Figure 7.** Surface topography from degree 1 to 12. The contour interval is 200 m.

model, the upwelling plumes in the upper mantle are not taken into account. Because the surface topography is strongly dependent on mantle density heterogeneities in the upper mantle (see Fig. 4b) such plumes will involve large surface radial displacement. Nevertheless these plumes are narrow and therefore too small to significantly affect large wavelength geodetic observables.

The mean value of this surface topography over the continents is about −90 m, one order of magnitude smaller than the continental elevation due to the isostasy adjustment. This is due to the fact that the oceanic crust is globally denser than the continental crust (see Section 3.3.2).

These mantle density heterogeneities are time-dependent and their temporal evolutions are mostly due to the subducted plates sinking into the mantle.

From the degree 1 and the degree 2 components of these mass anomalies, we can compute using our Love number formalism the secular variation of the zonal degree 2 geopotential coefficient $J_2$ and of the geocentre motion.

## 3.2 Geological variation of $J_2$ since 120 Ma

Starting with this model of time-dependent internal loads, and taking into account the viscoelastic deformations, we compute, since 120 Ma, the perturbations of the degree 2 geopotential coefficient $J_2$ using the degree 2 geoid kernel introduced in Section 2.

$$ \delta J_2(t) = -\frac{3}{\sqrt{5}} \int_b^a \left[ \left(\frac{r}{a}\right)^3 \delta(t) + k_2(r, t) \right] \frac{r}{a} \frac{\Delta \rho_2^{c0}(r, t)}{\tilde{\rho}a} dr. \quad (11) $$

The temporal evolution of $\delta J_2(t)$ since 120 Ma is plotted on Fig. 8. Our present computed value is $\delta J_2(0) = 0.57 \times 10^{-5}$.

The difference between the observed $J_2^{\text{obs}} = 0.001082626360$ (Foerste et al. 2008) and the hydrostatic one $J_2^H = 0.001072701000$ (Nakiboglu 1982) is about $10^{-5}$. Consequently, we can conclude that our simple model of mantle density heterogeneities can explain most of the part of the non-hydrostatic oblateness of the Earth. This result was already proposed by Forte et al. (1995). This contribution of the mantle mass anomalies to the degree 2 component of the geoid has to be added to other contributions such as the stabilization of continent-ocean heterogeneity which reduces the vertical flow of mass and heat across the upper mantle (Forte et al. 1995) or to the stabilizing effect of the elastic upper part of the lithosphere on the equilibrium rotation form of the Earth (Mitrovica et al. 2005).

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Temporal evolution of the degree 2 geopotential coefficient $J_2$ induced by our time-dependent internal loads model, since 120 Ma.</figcaption></figure>

To compare this effect with the one induced by the postglacial rebound, we use the simple model (three ice-sheets with the same dependence) we developed in (Greff-Lefftz 2000). For the same earth model (homogeneous mantle with a viscosity ratio $VR = 40$), we will obtain a perturbation of the zonal degree two coefficient of the geopotential $\delta J_2^{\text{Pgr}} = 0.1 \times 10^{-5}$, one order of magnitude smaller than the contribution of mantle mass anomalies.

Let us now investigate the secular variation of the $J_2$ coefficient, that is to say the time-derivative $\dot{J}_2$. Note, in our Fig. 8, that since 15 Ma, there is a linear trend in the perturbation $\delta J_2$ induced by the dynamics of the mantle, with the present rate:

$$ \dot{J}_2 = -1.17 \times 10^{-13} \text{ yr}^{-1}. $$

According to space geodetic observations over the past 25 yr, the zonal degree 2 coefficient of the Earth's gravitational potential $J_2$ has undergone a decrease, with a secular term of about $\dot{J}_2 = \sim -2.6 \times 10^{-11} \text{ yr}^{-1}$ (Foerste et al. 2008). Our computed value is less than 1 per cent of the observed one which is usually assumed to result from the postglacial rebound excitation source, which perturbs the geoid at a timescale faster than the one of the mantle dynamics. For our simple model of deglaciation (Greff-Lefftz 2000), with $VR = 40$, we would obtain $\dot{J}_2 = -4 \times 10^{-11} \text{ yr}^{-1}$, value a little bit too large to explain the observations but with the correct order of magnitude.

To end this study about the secular variation of the $J_2$ coefficient, we can conclude that, contrary to the surface load induced by the last deglaciation which perturb significantly the $\dot{J}_2$ but not the static term $\delta J_2$, the mantle density heterogeneities within the mantle can explain the present-day non-hydrostatic flattening of the Earth but that they vary at a too slow timescale to significantly perturb the $\dot{J}_2$ coefficient.

## 3.3 Geological variation of the geocentre

In this part, we investigate the influence of the mantle dynamics and the plate tectonics on the degree 1 deformations at the geological timescale.

### 3.3.1 The mantle dynamics

We define the geocentre as the figure centre of the deformed surface

$$ [X, Y, Z] = \frac{1}{4\pi} \int_{\Sigma} \vec{u}_1 \sin\theta d\theta \, d\phi, \quad (12) $$
<figure><img src="image_10.png" /><figcaption>Figure 9. Temporal evolution of the coordinates of the geocentre (left-hand panel), and present-day degree 1 surface displacement with contour interval of 50 m (right-hand panel), in metre, induced by our time-dependent internal loads model, since 120 Ma.</figcaption></figure>

where $\vec{u}_1$ is the degree 1 displacement vector at the Earth's surface $\Sigma$. X, Y, Z are the coordinates of the geocentre with respect to the initial reference frame centred to the centre of mass of the Earth and defined in Section 2.1.

Using the degree one viscoelastic Love numbers introduced in Section 2 for the radial displacement $h'_1(t, r)$ and the tangential displacement $l'_1(t, r)$, we can compute the geocentre perturbations induced by the mantle mass anomalies

$$ \begin{bmatrix} X(t) \\ Y(t) \\ Z(t) \end{bmatrix} = \begin{bmatrix} \frac{1}{\bar{\rho}\sqrt{3}} \int_b^a \left(\frac{r}{a}\right) [h'_1(t, r) + 2l'_1(t, r)] * \Delta\rho_1^{c1}(r, t) dr \\ \frac{1}{\bar{\rho}\sqrt{3}} \int_b^a \left(\frac{r}{a}\right) [h'_1(t, r) + 2l'_1(t, r)] * \Delta\rho_1^{s1}(r, t) dr \\ \frac{1}{\bar{\rho}\sqrt{3}} \int_b^a \left(\frac{r}{a}\right) [h'_1(t, r) + 2l'_1(t, r)] * \Delta\rho_1^{c0}(r, t) dr \end{bmatrix}. \quad (13) $$

We have shown in the first section that our theoretical approach with the internal load viscoelastic Love number formalism is not valid to compute the tangential displacement at the timescale of the mantle convection. Because we want simply to quantify the order of magnitude of the secular geocentre variation, we will do a simple approximation: it seems, from the study about the radial and tangential displacement (see Fig. 2), that the most important contribution of an internal load to the degree one surface displacement is the rigid one (i.e. the direct effect of the mass on the degree one surface displacement). Consequently, we assume, in this preliminary study that the mantle is rigid. In this case, the Love numbers become

$$ h'_1(r, t) = - \left(\frac{r}{a}\right)^2 \delta(t) \quad \text{and} \quad l'_1(r, t) = - \left(\frac{r}{a}\right)^2 \delta(t) $$

and the geocentre perturbations induced by internal loads may be written

$$ \begin{bmatrix} X(t) \\ Y(t) \\ Z(t) \end{bmatrix} = \begin{bmatrix} -\frac{\sqrt{3}}{\bar{\rho}} \int_b^a \left(\frac{r}{a}\right)^3 \Delta\rho_1^{c1}(r, t) dr \\ -\frac{\sqrt{3}}{\bar{\rho}} \int_b^a \left(\frac{r}{a}\right)^3 \Delta\rho_1^{s1}(r, t) dr \\ -\frac{\sqrt{3}}{\bar{\rho}} \int_b^a \left(\frac{r}{a}\right)^3 \Delta\rho_1^{c0}(r, t) dr \end{bmatrix}. \quad (14) $$

We use our model of mantle density heterogeneities and first plot, in Fig. 9, the temporal evolution of the geocentre coordinates since 120 Ma (left-hand panel) and the present-day surface displacement (right-hand panel) induced by these mantle mass anomalies.

We find that the present-day difference between the centre of figure of the outer surface and the centre of mass of the Earth, that is, the coordinates of the geocentre, are, in metre:

$$ [X(0), Y(0), Z(0)] = [316, -368, -256] \text{ m} $$

that is to say a translation of about 550 m towards the direction $-28^\circ S, 311^\circ E$.

We can compare this translation with the one induced by the last deglaciation: for our simple model (Greff-Lefftz 2000), we obtained a translation of about 7 m towards the direction $-70^\circ S, 134^\circ E$ (i.e. towards the opposite direction of Greenland). It is negligible in comparison with present-day surface displacement induced by the mantle mass anomalies.

The physical interpretation of this result is difficult: it corresponds to the variation of the figure centre between an Earth with mantle density heterogeneity and our reference earth model which is radially stratified. However it could be interesting to develop a geodetical method allowing to measure the difference between the centre of mass (inferred from the studies of the satellites trajectories) and the centre of figure (inferred from the location of the geodetical stations at the Earth's surface).

Another result is the secular variation of the geocentre. We plot on Fig. 10 the time-derivatives of the geocentre coordinates since 120 Ma: the order of magnitude is a 10 m Myr$^{-1}$, that is, a 1 mm Cy$^{-1}$ (one millimetre per century).

<figure><img src="image_11.png" /><figcaption>Figure 10. Time-derivative of the coordinates of the geocentre, in metre per million years (or in millimetre per century), since 120 Ma.</figcaption></figure>
Note on Fig. 10, that since 10 Ma, there is a constant rate of secular variation of the geocentre coordinates

$$ \dot{X} = 0.17 \text{ mm Cy}^{-1}; \quad \dot{Y} = -0.85 \text{ m Cy}^{-1}; \quad \dot{Z} = 1.12 \text{ mm Cy}^{-1} $$

These values are negligible in comparison with the geocentre motions induced by postglacial rebound, which has been shown to be at the level of -0.4-0.2 mm yr⁻¹ (Greff-Lefftz 2000), or by recent climate changes (about a few tenths of a millimetre per year).

The rigid translation induced by this mantle dynamic depends only on the degree one component of the mass anomalies, that is to say on the distribution of the plates within the mantle at large

wavelengths. The mantle viscosity plays a role only via the factor F, related to how plates cross the 670 km discontinuity. We have investigated the influence of such a factor on the associated variation in the geocentre motion and have found that it does not change the order of magnitude of our results: the present-day geocentre has an amplitude varying from 400 m (for F = 2.5) up to 600 m (for F = 4.5) in a direction varying in latitude from -22°S up to -34°S and in longitude from 306°E up to 314°E. The rate of the secular variation of the geocentre remains closed to one millimetre per century.

From our results we can conclude that the observed secular variation of the geocentre is not induced by internal loading. Surface loading seems to be a better candidate to create a rate of about 1 mm yr⁻¹.

<figure><img src="image_12.png" /><figcaption>Figure 11. Simple model of continental and oceanic crust.</figcaption></figure>

### 3.3.2 The plate tectonics

To end this study, let us investigate another phenomenon related to mantle dynamics.

Since the oceanic crust is globally denser than the continental crust, the isostasy adjustment creates a higher surface topography over the continents than over the oceans. We use a simple model of continental and oceanic crust (e.g. Turcotte & Schubert 1982; Lowrie 1997) plotted on Fig. 11.

The study of the continental crust gives for the density profile

1. 0- > 20 km: $\rho = 2750 \text{ kg m}^{-3}$ (granitic rock and gneiss)

2. 20- > 35 km: $\rho = 2950 \text{ kg m}^{-3}$ (gabbros)

3. up to 35 km: $\rho = 3310 \text{ kg m}^{-3}$ (peridotite).

The study of the oceanic crust gives for the density profile

1. 0- > 4 km: $\rho = 1000 \text{ kg m}^{-3}$ (water)

2. 4- > 5 km: $\rho = 1500 \text{ kg m}^{-3}$ (water + sediments)

<figure><img src="image_13.png" /><figcaption>Figure 12. Continent distribution from 200 Ma to today.</figcaption></figure>
3. 5– > 11 km: $\rho$ = 2850 kg m⁻³ (basalts and gabbros)

4. 11– > 35 km: $\rho$ = 3310 kg m⁻³ (peridotite).

Because of the isostatic equilibrium at the depth of 35 km where the weight of a continental crust column and the one of an oceanic crust column are equal, there is an elevation $h$ = 843 m of the continental crust with respect to the oceanic one. Such isostatic topography does not modify the geoid in our first order perturbation theory of the viscoelasto-gravity and consequently the associated perturbation of the $J_2$ coefficient is zero. It will imply a discrepancy between the centre of figure of the surface and the centre of mass due to the associated degree one term in the continental topography.

The plate tectonic slow motions of continents over time (typically at a cm yr⁻¹ level, see Nuvell-NNR for example; Altamimi et al. 2007) should change the geocentre position, on a geological timescale.

Let us introduce the ocean–continent function $f(\theta, \varphi)$ which is equal to 1 over the ocean and to 0 over the continent. It may be

expanded in spherical harmonics

$$f(\theta, \varphi) = \sum_{n=0}^{\infty} \sum_{m=0}^{n} (a_{nm} \cos m\varphi + b_{nm} \sin m\varphi) P_{n}^{m}(\cos \theta), \quad (15)$$

where $a_{nm}$ and $b_{nm}$ are the coefficients of the spherical harmonics expansion, normalized to $4\pi$. The surface topography may be written as $h[1 - f(\theta, \varphi)]$ and consequently the geocentre motion is

$$\begin{bmatrix} X(t) \\ Y(t) \\ Z(t) \end{bmatrix} = \frac{h}{\sqrt{3}} \begin{bmatrix} a_{11}(t) \\ b_{11}(t) \\ a_{10}(t) \end{bmatrix}. \quad (16)$$

In order to quantify this phenomenon, we investigated the position of continents over the last 200 Myr, using information provided by hotspot reference frame, palaeomagnetism and kinematic data from the sea floor. Fig. 12 shows the ocean–continent repartition over time that we used in the present study. We compute the coefficients $a_{10}(t)$, $a_{11}(t)$ and $b_{11}(t)$ associated with such a motion.

<figure><img src="image_14.png" /><figcaption>Figure 13. Centre of mass wrt geocentre position and velocity induced by the continent motion from 200 Ma to today.</figcaption></figure>
Assuming as in the present-day situation that the continents present a mean topography of 843 m with respect to the mean sea level, we calculated the motion and the velocity of the geocentre over time. Fig. 13 presents the variations of geocentre position (X, Y, Z) on top and the velocity variations (Ẋ, Ẏ, Ż) on bottom. One can see that the presence of continents and their motion induce a departure between the geocentre and the centre of mass relatively large, about a few tens metres, with a maximum of 130 m on X component during the Mesozoic when there was only one continent on Earth (the Pangea). The present-day values are, in metre:

$$[X(0), Y(0), Z(0)] = [-53, -28, -58] \text{ m}$$

that is to say a translation of about 84 m towards the direction -62°S, 208°E.

The geocentre velocity is between -0.05 and 0.05 mm Cy⁻¹, which is smaller than the impact of mantle dynamics.

Let us end this section with comments about the mean elevation of the continents with respect to the oceans.

The mean topography over the continent induced by isostatic adjustment between oceanic and continental crust is about 843 m, although the one induced by mantle mass anomalies is negative and about of -90 m. The sum of this two effects leads to a mean elevation of about 753 m. The mean topography on the continent obtained from the observed surface topography is 793 m, based on ETOPO global relief model (Amante & Eakins 2008). Consequently, it remains a small positive part of the mean topography on the continent, which is not explained by our large-wavelength model. It may be induced by mantle density anomalies within the upper part of the mantle and especially by the narrow hot plumes existing in the upper mantle.

## 4 CONCLUSION

We quantified the J₂ gravitational potential coefficient and the geocentre motion induced by upwelling domes and subducted plates sinking into the mantle, using a simple model in which the slabs are modelled as blobs diving vertically through the mantle, and in which the domes are assumed to be stable over the last 120 Ma and modelled from tomography within the lower mantle.

We can conclude that the density heterogeneities within the mantle can explain an important part of the present-day non-hydrostatic flattening of the Earth but that they vary at a timescale which is too slow to significantly perturb the J₂ coefficient. Although there is a significant discrepancy of about a few hundred metres between the centre of figure and the centre of mass of the Earth, the secular variation of the geocentre motion is one order of magnitude smaller than the one induced by surface loads.

We also quantified the impact of continents motions on the geocentre motion: the present-day geocentre velocity is less than 0.05 mm Cy⁻¹ and consequently negligible.

The combination of the effects induced by the mantle density heterogeneities and those resulting from plates tectonics involves a discrepancy between the centre of figure and the centre of mass of the Earth about 570 m in the direction -33°S, 303E.

## ACKNOWLEDGMENTS

This study is IPGP contribution number 2590. We thank Hilaire Legros and Olivier de Viron for discussions and remarks on the

original manuscript, Hélène Rouby for the mantle density heterogeneities model and Y. Ricard for giving us his geodynamic model.

## REFERENCES

Altamimi, Z., Collilieux, X., Legrand, J., Garayt, B. & Boucher, C., 2007. ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth orientation parameters, *J. geophys. Res.*, **112**, B09401, doi:10.1029/2007JB004949.

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

Amante, C. & Eakins, B.W., 2008. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis, National Geophysical Data Center, NESDIS, NOAA, U.S. Department of Commerce, Boulder, CO, August 2008.

Besse, J. & Courtillot, V., 2002. Apparent and true polar wander and the geometry of the geomagnetic field over the last 200 Myr, *J. geophys. Res.*, **107**(B11), 2300, doi:10.1029/2000JB000050.

Bills, B.G. & James, T.S., 1997. Polar motion of a viscoelastic Earth due to glacial cycle mass loading, *J. geophys. Res.*, **102**(B4), 7579–7602.

Blewith, G., 2003. Self-consistency in reference frames, geocenter definition, and surface loading of the solid Earth, *J. geophys. Res.*, **108**(B2), 2103, doi:10.1029/2002JB002082.

Cazenave, A. & Nerem, S., 2002. Redistributing Earth’s mass, *Science*, **297**, 783–784.

Chen, J.L., Wilson, C.R., Eanes, R.J. & Nerem, R.S., 1999. Geophysical interpretation of observed geocenter variations, *J. geophys. Res.*, **104**, 2683–2690.

Chinnery, M.A., 1975. The static deformation of an Earth with a fluid core: a physical approach, *Geophys. J. R. astr. Soc.*, **42**, 461–475.

Collilieux, X., Altamimi, Z., Ray, J., van Dam, T. & Wu, X., 2009. Effect of the satellite laser ranging network distribution on geocenter motion estimation, *J. geophys. Res.*, **114**, B04402, doi:10.1029/2008JB005727.

Cox, C.M. & Chao, B.F., 2002. Detection of a large-scale mass redistribution in the terrestrial system since 1998, *Science*, **297**, 831–833.

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.

Dickey, J.O., Marcus, S.L., de Viron, O. & Fukumori, I., 2002. Recent Earth oblateness variations: unraveling climate and postglacial rebound effects, *Science*, **298**, 1975–1977.

Dong, D., Dickey, J.O., Chao, Y. & Cheng, M.K., 1997. Geocenter variations caused by atmosphere, ocean and surface ground water, *Geophys. Res. Letts.*, **24**, 1867–1870.

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

Farrell, W.E., 1972. Deformation of the Earth by surface loads, *Rev. Geophys. Space Phys.*, **10**(3), 761–797.

Foerste, C. et al., 2008. The GeoForschungsZentrum Potsdam/Groupe de Recherche de Geodesie Spatiale satellite-only and combined gravity field models: EIGEN-GL04S1 and EIGEN-GL04C. *J. Geod.*, **82**(6), 331–346, doi:10.1007/s00190-007-0183-8.

Forte, A.M. & Mitrovica, J.X., 1996. New inferences of mantle viscosity from joint inversion of long-wavelength mantle convection and postglacial rebound data, *Geophys. Res. Lett.*, **23**, 1147–1150.

Forte, A.M., Dziewonski, A.M. & O’Connell, R.J., 1995. Thermal and chemical heterogeneity in the mantle—a seismic and geodynamic study of continental roots, *Phys. Earth planet. Int.*, **92**(1-2), 45–55.

Fukao, Y., Widiyantoro, S. & Obayashi, M., 2001. Stagnant slabs in the upper and lower mantle transition region, *Rev. geophys. Res. Lett.*, **15**, 229–232.

Greff-Lefftz, M., 2000. Secular variation of the geocenter, *J. geophys. Res. Solid Earth*, **105**(B11), 25685–25692.

Greff-Lefftz, M. & Legros, H., 1997. Some remarks about the degree-one deformation of the Earth, *Geophys. J. Int.*, **131**, 699–723.
Greff-Lefftz, M. & Legros, H., 2007. Fluid core dynamics and degree-one deformations: Slichter mode and geocenter motions, *Phys. Earth planet. Int.*, **161**, 150–160.

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.

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

Longman, I.M., 1962. Greens function for determining deformation of Earth under surface mass loads. 1: theory, *J. geophys. Res.*, **67**(2), 845–850.

Lowrie, W., 1997. *Fundamentals of Geophysics*, Cambridge University Press, Cambridge, UK.

Marquart, G., Steinberger, B. & Niehuus, K., 2205. 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.

Métivier, L., Greff-Lefftz, M., Altamimi, Z. & Besse, J., 2009. Secular motions of the geocenter, *Geophys. Res. Abstracts*, **11**, EGU2009–9542.

Mitrovica, J.X. & Peltier, W.R., 1991. A complete formalism for the inversion of post-glacial rebound data: resolving power analysis, *Geophys. J. Int.*, **104**, 267–288.

Mitrovica, J.X. & Peltier, W.R., 1993. The inference of mantle viscosity from an inversion of the Fennoscandian relaxation spectrum, *Geophys. J. Int.*, **114**, 45–62.

Mitrovica, J.X., Wahr, J., Matsuyama, I. & Paulson, A., 2005. The rotational stability of an ice-age earth, *Geophys. J. Int.*, **161**(2), 491–506.

Nakiboglu, S.M., 1979. Hydrostatic figure and related properties of the Earth, *Geophys. J. R. astr. Soc.*, **57**, 639–648.

Nakiboglu, S.M., 1982. Hydrostatic theory of the Earth and its mechanical implications, *PEPI*, **28**, 302–311.

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

Peltier, W.R., 1985. The LAGEOS constraint on deep mantle viscosity— results from a new normal mode method for the inversion of viscoelastic relaxation spectra, *J. geophys. Res.*, **90**(B11), 9411–9421.

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

Ricard, Y., Mattern, E. & Matas, J., 2005. Tomographic Images of Slabs from Mineral Physics, *Geophys. Monogr.*, **160**, 283–300.

Richards, M.A., Ricard, Y., Lithgow-Bertelloni, C., Spada, G. & Sabadini, R., 1997. An Explanation for Earth's Long-Term Rotational Stability, *Science*, **275**, 372–375.

Rouby, H., Greff-Lefftz, M. & Besse, J., 2010. Mantle dynamics, geoid, inertia and TPW since 120 Ma, *E.P.S.L.*, in press.

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. & O'Connell, R., 1997. Changes of the Earth's rotation axis owing to advection of mantle density heterogeneities, *Nature*, **387**, 169–173.

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

Steinberger, B. & O'Connell, R., 2002. The convective mantle flows signal in rates of True Polar Wander, *Ice Sheets, Sea Level and the Dynamic Earth*, Geodynamics Series 29, Copyright 2002 by the American Geophysical Union.

Turcotte, D.L. & Schubert G., 1982. *Geodynamics. Applications of Continuum Physics to Geological Problems*, John Wiley & Sons, New York.

Tushingham, A.M. & Peltier, W.R., 1991. Ice-3g: a new global model of late Pleistocene deglaciation based upon geophysical predictions of post-glacial relative sea level change, *J. geophys. Res.*, **96**, 4497–4523.

Vermeersen, L.L.A., Sabadini, R., Devoti, R., Luceri, V., Rutigliano, P., Sciarretta, C. & Bianco, G., 1998. Mantle viscosity inferences from joint inversions of Pleistocene deglaciation-induced changes in geopotential with a new SLR analysis and polar wander, *Geophys. Res. Lett.*, **25**, 4261–4264.

## APPENDIX A

The elastogravitational system can be written, for each degree $n \ge 1$ as a first order differential equations in a form given by Alterman et al. (1959)

$$ \frac{dy_{in}(r)}{dr} = A_{ij}^{n} y_{jn}(r), \qquad (A1) $$

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

Our earth model consists of $N$ homogeneous incompressible layers: a lithosphere (layer 1), a stratified mantle (layer 2 to $N-2$), an inviscid fluid core (layer $N-1$) and an inner core (layer $N$). The density and rigidity stratification within the outer core and the mantle are obtained from the PREM earth model (Dziewonski & Anderson 1981).

In each homogeneous incompressible solid layer $i$ (that is in each layer except for the fluid core) with a density $\rho^i$ and a constant rigidity $\mu^i$, solving the elasto-gravitational system, we obtain

$$ y_{1n}^{i}(r) = \frac{C_{1n}^{i}}{r^{n}} + \frac{C_{2n}^{i}}{r^{n+2}} + C_{3n}^{i} r^{n+1} + C_{4n}^{i} r^{n-1} $$

$$ y_{5n}^{i}(r) = C_{5n}^{i} r^{n} + \frac{C_{6n}^{i}}{r^{n+1}} $$

$$ y_{3n}^{i}(r) = -\frac{n-2}{n(n+1)} \frac{C_{1n}^{i}}{r^{n}} - \frac{1}{n+1} \frac{C_{2n}^{i}}{r^{n+2}} + \frac{n+3}{n(n+1)} C_{3n}^{i} r^{n+1} + \frac{C_{4n}^{i}}{n} r^{n-1} $$

$$ y_{2n}^{i}(r) = 2\mu^{i} \left[ \frac{-n^2 - 3n + 1}{n+1} \frac{C_{1n}^{i}}{r^{n+1}} - (n+2) \frac{C_{2n}^{i}}{r^{n+3}} + \frac{n^2 - n - 3}{n} C_{3n}^{i} r^{n} + (n-1) C_{4n}^{i} r^{n-2} \right] + \rho^{i} g(r) y_{1n}^{i}(r) - \rho^{i} y_{5n}^{i}(r) $$

$$ y_{4n}^{i}(r) = 2\mu^{i} \left[ \frac{n-1}{n} \frac{C_{1n}^{i}}{r^{n+1}} + \frac{n+2}{n+1} \frac{C_{2n}^{i}}{r^{n+3}} + \frac{n+2}{n+1} C_{3n}^{i} r^{n} + \frac{n-1}{n} C_{4n}^{i} r^{n-2} \right] $$

(A2)

$$ y_{6n}^{i} = n C_{5n}^{i} r^{n-1} - (n+1) \frac{C_{6n}^{i}}{r^{n+2}} - 4\pi G \rho^{i} y_{1n}^{i}(r). $$
In each layer i there are 6 constants ($C_{1n}^i, \dots, C_{6n}^i$) which will be determined from the boundary conditions and consequently will be dependent on the excitation sources deforming the planet.

Within the inviscid fluid core, the differential system of 6 equations is not valid. We use the theoretical approach introduced by Chinnery (1975): the outer core is assumed to be in hydrostatic equilibrium, that is to say the form is the equipotential [$y_{5n}^{N-1}(r) = g(r) y_{1n}^{N-1}(r)$ and $y_{2n}^{N-1}(r) = 0$] and the tangential displacement $y_{3n}^{N-1}(r)$ is undetermined. For an homogeneous incompressible fluid core, we obtain:

$$y_{5n}^{N-1}(r) = C_{5n}^{N-1} r^n + \frac{C_{6n}^{N-1}}{r^{n+1}}$$

$$y_{1n}^{N-1}(r) = y_{5n}^{N-1}(r)/g(r)$$

$$y_{6n}^{N-1} = n C_{5n}^{N-1} r^{n-1} - (n+1) \frac{C_{6n}^{N-1}}{r^{n+2}} - 4\pi G \rho^i y_{1n}^{N-1}(r).$$

There are two constants ($C_{5n}^{N-1}$ and $C_{6n}^{N-1}$) which will be determined from the boundary conditions.

## A1 Boundary conditions

We are interested in the deformations induced by a surface load as well as by internal loading.

A degree n surface mass distribution $\sigma_n$, located at the radius $r_o$ produces two effects.

(i) a degree n pressure effect $-g(r_o) \sigma_n$ acting at the interface $r_o$

(ii) a degree n potential effect $\frac{4\pi G}{2n+1} \sigma_n r_o$ 
$$ \begin{cases} \left(\frac{r}{r_o}\right)^n & \text{if } r \le r_o \\ \left(\frac{r_o}{r}\right)^{n+1} & \text{if } r \ge r_o. \end{cases} $$

We note $S_n^{\text{int}} = \frac{4\pi G}{2n+1} \sigma_n r_o$ the internal load, and $S_n^e = \frac{4\pi G}{2n+1} \sigma_n a$ the surface load, where a is the surface radius.

In the above section the subscript n will be suppressed for simplicity of notation.

When the Earth is submitted to these different excitation sources, the boundary conditions may be written.

1. At the surface: the tangential stress is equal to the acting tangential traction, that-is-to-say zero and the radial stress is equal to the acting pressure; the discontinuity in the gravitational potential leads to a relation between $y_6(a)$, $y_5(a)$ and the acting potential (Longman 1962; Farrell 1972):

$$ \begin{cases} y_4^1(a) = 0 \\ y_2^1(a) = -\frac{2n+1}{3} \tilde{\rho} S^e \\ y_6^1(a) + \frac{n+1}{a} y_5^1(a) = \frac{2n+1}{a} S^e \end{cases} $$

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

2. At an interface $a_i$ within the mantle: the radial and tangential displacements, the tangential stress and the potential are continuous; taking into account an internal loading potential $S^{\text{int}}$ acting at the radius $r_o$ the continuity of the radial stress and that of the gravity at an interface $a_i$ will depend on this excitation source if $a_i = r_o$,

$$ \begin{cases} y_1^i(a_i) = y_1^{i+1}(a_i) \\ y_2^i(a_i) = y_2^{i+1}(a_i) + \frac{2n+1}{3} \frac{g(a_i)}{g(a)} \frac{a}{a_i} \tilde{\rho} S^{\text{int}} \delta(a_i - r_o) \\ y_5^i(a_i) = y_5^{i+1}(a_i) \\ y_3^i(a_i) = y_3^{i+1}(a_i) \\ y_4^i(a_i) = y_4^{i+1}(a_i) \\ y_6^i(a_i) = y_6^{i+1}(a_i) - \frac{2n+1}{a_i} S^{\text{int}} \delta(a_i - r_o), \end{cases} $$

where $\delta(a_i - r_o)$ is the Dirac function which is equal to 1 when $a_i = r_o$ and zero when $a_i \neq r_o$.

3. At the CMB: the radial displacement is continuous but we have to introduce an unknown parameter $K_1$ in order to take into account the difference between the form of the solid mantle and the form of the fluid core which is an equipotential (Chinnery 1975); the tangential displacement is not determined; the tangential stress is equal to zero; the gravitational potential is continuous; the continuity of the radial stress depends on fluid pressure acting within the core and the gravitational attraction of the mantle is related to the one of the fluid core,
taking into account the discrepancy $K_1$

$$
\left\{
\begin{aligned}
y_1^{N-2}(a_{N-2}) &= y_1^{N-1}(a_{N-2}) + a_{N-2}K_1 \\
y_3^{N-2}(a_{N-2}) &= a_{N-2}K_3 \\
y_4^{N-2}(a_{N-2}) &= 0 \\
y_5^{N-2}(a_{N-2}) &= y_5^{N-1}(a_{N-2}) \\
y_2^{N-2}(a_{N-2}) &= \rho^{N-1}g(a_{N-2})a_{N-2}K_1 \\
y_6^{N-2}(a_{N-2}) &= y_6^{N-1}(a_{N-2}) - 4\pi G\rho^{N-1}a_{N-2}K_1
\end{aligned}
\right.
$$

4. At the ICB: the radial displacement is continuous but we have to introduce an unknown parameter $K_2$ in order to take into account the difference between the form of the solid inner core and the form of the fluid core which is an equipotential; the tangential displacement is not determined; the tangential stress is equal to zero at the ICB; the gravitational potential is continuous; the continuity of the radial stress depends on fluid pressure within the core acting at the ICB and the gravitational attraction of the inner core is related to the one of the fluid core, taking into account the discrepancy $K_2$

$$
\left\{
\begin{aligned}
y_1^{N-1}(a_{N-1}) &= y_1^N(a_{N-1}) - a_{N-1}K_2 \\
y_3^{N-1}(a_{N-1}) &= a_{N-1}K_4 \\
y_4^{N-1}(a_{N-1}) &= 0 \\
y_5^{N-1}(a_{N-1}) &= y_5^N(a_{N-1}) \\
y_2^{N-1}(a_{N-1}) &= -\rho^{N-1}g(a_{N-1})a_{N-1}K_2 \\
y_6^{N-1}(a_{N-1}) &= y_6^N(a_{N-1}) + 4\pi G\rho^{N-1}a_{N-1}K_2.
\end{aligned}
\right.
$$

For the degree $n = 1$, Farrell (1972) has shown that there is a Consistency Relation, that is, a special condition in which the degree one valid solutions have to obey, which is written, in the layer $i$, when the Earth is submitted to surface or internal loads

$$y_2^i(r) + 2y_4^i(r) + \frac{g(r)}{4\pi G} \left[ y_6^i(r) + \frac{2}{r} y_5^i(r) \right] = 0.$$

Because of that, only two of the three surface conditions are needed and the Consistency Relation ensures that the third boundary condition is met automatically. To solve this problem, we have to add a new boundary condition: the conservation of the centre of mass of the Earth plus internal load. It will simply require that the degree one surface potential is equal to zero, that is to say, we just have to add a surface boundary condition:

$$y_5^i(a) = 0.$$

These 6 $N-3$ boundary conditions for surface and internal loads allow us to determine the 6($N-2$) constants ($C_{1n}^i, \dots, C_{6n}^i$) for $i=1, \dots, N-2$ of the mantle, the six constants of the fluid core $C_{5n}^{N-1}$, $C_{6n}^{N-1}$, $K_1, K_2, K_3, K_4$, and the three constants $C_{3n}^N$, $C_{4n}^N$ and $C_{5n}^N$ of the inner core. The constants known, we can compute the surface displacement and the gravitational potential and deduce the internal loading Love numbers.