# Upwelling plumes, superswells and true polar wander

Marianne Greff-Lefftz

► To cite this version:

Marianne Greff-Lefftz. Upwelling plumes, superswells and true polar wander. Geophysical Journal International, 2004, 159, pp.1125-1137. 10.1111/j.1365-246X.2004.02440.x. insu-03600273

HAL Id: insu-03600273

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

Submitted on 7 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
# Upwelling plumes, superswells and true polar wander

Marianne Greff-Lefftz

Institut de Physique du Globe de Paris, Département de Géomagnétisme et Paléomagnétisme, 4 place Jussieu, 75252 Paris 05, France.

E-mail: greff@ipgp.jussieu.fr

Accepted 2004 July 23. Received 2004 July 22; in original form 2003 July 23

## SUMMARY

The geological evolution of the rotational axis of the Earth is most likely controlled by internal mass redistribution within the mantle. Palaeomagnetic observations suggest that it is episodic in nature, with periods of quasi-standstill alternating with periods of faster wander. Here, we investigate two models for the influence of mantle plumes that vary at different spatial wavelengths on the time variations of the rotational axis (true polar wander; TPW). In the first model, we represent an upwelling plume as a sphere whose radius varies as a function of the flux of material in the conduit and that traverses the mantle at the Stokes velocity. Such a plume produces very little wander of the rotational axis. We then study the effects of two superswells that mimic the ones observed with seismic tomography and conclude that a doming regime within the mantle involves significant polar wander. Some of the features of this TPW that are directly linked to the periodicity of doming are reminiscent of observed phases of slow and fast TPW, with similar peak velocities.

**Key words:** geodynamics, plumes, polar motion.

## 1 INTRODUCTION

True polar wander (TPW) is the coherent motion of Earth with respect to its rotation axis. In view of the angular momentum theorem, it may be interpreted physically as the motion of the rotational axis with respect to the fixed terrestrial frame described by hotspots. On geological timescales, this motion is most likely controlled by internal mass redistribution within the mantle.

Palaeomagnetic evidence from recent Earth history suggests that this motion has been small. Palaeomagnetic studies combine apparent polar wander based on palaeomagnetic data (the palaeomagnetic reference frame) with motions of hotspots with respect to lithospheric plates (the hotspot reference frame). Given the relative motions between the plates, the position of these plates with respect to the rotation axis (provided by palaeomagnetic results), and the relative motions between plates and a reference frame attached to the mantle, TPW can be defined as the motion of the mantle with respect to the rotation axis. A recent study (Besse & Courtillot 2002), hereafter referred to as BC02, confirms earlier findings that TPW appears to be episodic in nature, with periods of quasi-standstill alternating with periods of faster TPW (Fig. 1). Specifically, the path back to 200 Ma shows a standstill at 160–130 Ma, then a circular track from 130 to 70 Ma and finally a standstill at 50–10 Ma. As a result of the uncertainties in models of hotspot kinematics prior to 130 Ma, the behaviour of the TPW prior to that time is doubtful. The major event since 130 Ma is the end of the 130–60 Ma period of relatively fast polar wander (where the TPW rate averages 30 km Ma⁻¹), with a standstill (i.e. no or little TPW) from 50 Ma (possibly as early as 80 Ma because of larger 95 per cent confidence circles) to 10 Ma. It seems that TPW may have been negligible for approximately 50 Ma, but accelerated a few millions years ago, with a velocity of the order of 100 km Ma⁻¹; this recent acceleration is probably a result of the glacial effect.

The most plausible explanation for the slow TPW rate of the Earth is probably related to the great stability of subduction zones and the large-scale heterogeneity of the mantle. The effects of these on TPW have been explored in models in which subduction history or advected mass anomalies inferred from seismic tomography and global plate motions are used to estimate the time evolution of mantle heterogeneity for the Cenozoic and Mesozoic (Spada et al. 1992; Richards et al. 1997; Steinberger & O’Connell 1997; Richards et al. 1999). A recent study (Steinberger & O’Connell 2002) showed that large-scale upwellings (superswells) implied by tomography may play an important role in TPW. In this paper, we also study the effects of superswells on TPW, but we use an analytical approach to understand how the location of the domes close to the equator and their periodicity influences the episodic behaviour of TPW.

There have been suggestions for superfast episodes of TPW or excursions (see references in BC02). Most of these are found only by one author, on one subset of the data and do not seem to be robust or reproducible. BC02 conclude that none of these events can at present be documented with any certainty. However, because of the way the data are averaged in 10 to 20 Myr long time windows, this may lead to smoothing that could hide brief, 1 Myr long, excursions. Because of this uncertainty and because of the brevity of appearance of hotspots as
<figure><img src="image_2.png" /><figcaption>Observed TPW
(Besse and Courtillot, 2002)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Observed palaeomagnetic TPW (Besse & Courtillot 2002). (a) 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. (b) Temporal evolution of the latitude, longitude and velocity of the shift of the rotational axis at the surface of the Earth.</figcaption></figure>

traps at the surface of the Earth (e.g. for a review, Courtillot & Renne 2003; Courtillot 2004), it is interesting to investigate if mantle plumes could be able to generate TPW excursions.

The paper is organized as follows. In Section 1, we recall the equations governing the rotational dynamics of a viscous planet. In Section 2, we consider the influence on the polar motion of a plume crossing the mantle in order to explain fast excursions of TPW. Using a simple model in which the plume is represented as a sphere whose radius varies as a function of the flux of material in the conduit and that traverses the mantle at the Stokes velocity, we have investigated how the density jump and viscosity contrast at 670 km depth and the viscosity of the lithosphere influence the perturbations of the rotational axis. In Section 3, we extend these results to the case of two superswells, one beneath Africa and the other beneath Polynesia, modelled as domes oscillating vertically within the mantle with a period of 100 Myr. The conclusions are presented in Section 4.

## 2 THEORY OF THE ROTATION OF THE EARTH ON GEOLOGICAL TIMESCALE

We recall now the equations governing the rotational dynamics of a viscous planet. The temporal evolution of the rotation axis of the Earth is governed by the conservation of angular momentum **H** in the absence of external torques. In a reference frame corotating with the mantle,

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

with **H** = **J** · ω. **J** is the inertia tensor of the Earth and ω is the angular velocity. 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 δ$_{ij}$ the Kronecker symbol.

The first term $I_o$ is the inertia tensor of a spherical non-rotating Earth ($I_o = 0.33 Ma^2$ for the actual Earth). The second term $C_{ij}$ represents the centrifugal potential: the difference between the associated principal moments of inertia represents the rotational bulge. The last term $I_{ij}$ describes changes in inertia as a result of 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 gravity constant and * denotes time convolution. From the fluid limit of this eq. (3), we have $\frac{a^5}{3G} \simeq \frac{\alpha I_o}{\Omega^2 k^T}$, where k^T is the tidal fluid Love number, a is the radius of the Earth, Ω is the sidereal rotation rate and α is the hydrostatic flattening.

The tidal Love number k^T may be written in the time domain as (e.g. Peltier 1974; Spada et al. 1990)

$$ k^T(t) = k^e \delta(t) + \sum_{j=1}^{M} r_j e^{-t/\tau_j} H(t), \qquad (4) $$
<table><caption>**Table 1.** Geometrical and physical parameters of the five-layer reference Earth model (PREM-averaged values). For the model without a density jump at the 670-km-depth discontinuity, the density of the mantle and of the D'' layer is equal to 4558 kg m⁻³.</caption><thead><tr><th>Radius in km</th><th>Density in kg m<sup>-3</sup></th><th>Rigidity in Pa</th><th>Viscosity in Pa s</th></tr></thead><tbody><tr><td>a<sub>l</sub> = 6221 &lt; r &lt; a</td><td>3232</td><td>6.114 &times; 10<sup>10</sup></td><td>10<sup>26</sup></td></tr><tr><td>a<sub>trans</sub> = 5701 &lt; r &lt; a<sub>l</sub></td><td>3666</td><td>9.169 &times; 10<sup>10</sup></td><td>10<sup>21</sup></td></tr><tr><td>a<sub>D''</sub> = 3771 &lt; r &lt; a<sub>trans</sub></td><td>4904</td><td>2.225 &times; 10<sup>11</sup></td><td>3 &times; 10<sup>22</sup></td></tr><tr><td>b = 3480 &lt; r &lt; a<sub>D''</sub></td><td>4904</td><td>2.225 &times; 10<sup>11</sup></td><td>10<sup>21</sup></td></tr><tr><td>0 &lt; r &lt; b</td><td>10901</td><td>0</td><td>0</td></tr></tbody></table>

<table><caption>**Table 2.** Relaxation times (units of 10³ yr) for the Earth model described in Table 1.</caption><thead><tr><th></th><th>T<sub>r1</sub></th><th>T<sub>r2</sub></th><th>T<sub>r3</sub></th><th>T<sub>r3</sub></th><th>M<sub>o</sub></th><th>C</th><th>L</th><th>M<sub>1</sub></th></tr></thead><tbody><tr><td>τ (kyr)</td><td>0.151</td><td>0.199</td><td>0.358</td><td>0.403</td><td>9.034</td><td>10.800</td><td>57.454</td><td>1822.537</td></tr></tbody></table>

where $k^e$ is the tidal elastic Love number, $τ_j$ are the M relaxation modes associated with the Earth model, $r_j$ are the viscous tidal coefficients and $H(t)$ is the Heaviside step function. Our Earth model consists of five incompressible homogeneous layers: a viscoelastic lithosphere (viscosity $ν^L$), viscoelastic upper and lower mantles (viscosity $ν_2$ and $ν_3$), a viscoelastic D'' layer and an inviscid fluid core. This model has $M = 9$ or 10 relaxation modes (Peltier 1974): $M_o$ and $C$ resulting from the density contrast at the surface of the Earth and at the core–mantle boundary (CMB); $L$, $L_1$ and $L_2$ resulting from the viscosity and density contrast between the lithosphere and the upper mantle (for an elastic lithosphere, only the mode $L$ resulting from the viscosity contrast between the lithosphere and the viscoelastic upper mantle exists); $M_1$ resulting from the possible non-adiabatic density jump at the 670-km-depth discontinuity (the existence of such a mode is still controversial because it requires that the 670-km-depth discontinuity is a chemical boundary); transition modes $T_{r1}$, $T_{r2}$, $T_{r3}$ and $T_{r4}$ associated with the jump in the Maxwell times between upper and lower mantle and lower mantle and D'' layer. In this paper, we systematically investigate two mantle models: one with a density jump at 670 km depth (i.e. with a $M_1$ mode) and one without density jump. The geometrical and physical parameters of our five-layers Earth model and the M relaxation times (in 10³ yr) are listed in Tables 1 and 2, respectively.

For the long timescales associated with mantle convection, the Love number $k^T(t)$ is approximately equal to the quasi-fluid Love number. We have shown (Lefftz 1991) that the viscoelastic relaxation times resulting from a density jump within the mantle ($M_1$ mode) and the viscoelasticity of the lithosphere ($L_1$ and $L_2$ modes) are not significantly excited by the rotational potential, so that their viscous amplitudes are negligible. Assuming therefore that $t/τ_j \gg 1$, we obtain (Munk & MacDonald 1960; Lefftz 1991; Ricard et al. 1993a)

$$k^T(t) = k^f[\delta(t) - T_1\dot{\delta}(t)] \quad \text{with} \quad T_1 = \frac{1}{k^f} \sum_{j=1}^{M} r_j \tau_j^2, \qquad (5)$$

where $k^f = k^e + \sum_{j=1}^M r_j τ_j$ is the tidal fluid Love number. If the Earth is an inviscid fluid, the time $T_1$ is identically zero and the rotational bulge will be instantaneously perpendicular to the rotational axis. $T_1$ is the characteristic time of the delayed readjustment of the equatorial bulge of a rotating planet toward this fluid position. Its value, for the Earth model described before, depends on the viscosity profile within the mantle: $T_1$ increases with the viscosity ratio between lower and upper mantle, and is of the order 2–40 kyr for ratios in the range 5–100. Because the $M_1$ mode is not significantly excited by a tidal potential, the time $T_1$ does not depend on the density stratification at 670 km depth (e.g. Ricard et al. 1993a). With this quasi-fluid approximation, the inertia perturbation resulting from the centrifugal potential becomes

$$C_{ij}(t) = \frac{\alpha I_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]. \qquad (6)$$

On geological timescales, $dH/dt \ll \omega \wedge 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 this system, we can show that neglecting $dH/dt$ forbids any temporal variation of the amplitude of the rotation vector, implying $\omega^2 = \Omega^2$: the length of day cannot vary with this approach. The equations $\omega \wedge H = 0$ governing polar wander may now be written, using indicial notation:

$$\epsilon_{ijk} \omega_j \dot{\omega}_k = \epsilon_{ijk} \frac{I_{kl}}{\alpha I_o T_1} \omega_l \omega_j. \qquad (7)$$

The excitation function of TPW depends on two terms: the relaxation of the rotational bulge ($\alpha I_o T_1$), and the amplitude and the temporal evolution of the inertia tensor perturbations induced by the internal load ($I_{kl}$).

The mantle mass anomaly at each radius $r$ may be expressed as a load with surface density $\sigma(r, \theta, \varphi)$ that can be expanded in spherical harmonics:

$$\sigma(r, \theta, \varphi) = \sum_{n=0}^{\infty} \sum_{m=0}^{n} [\sigma_n^{cm}(r) \cos m\varphi + \sigma_n^{sm}(r) \sin m\varphi] P_n^m(\cos\theta),$$

where $P_n^m$ are the Legendre functions, $\theta$ is colatitude and $\varphi$ is longitude.
By MacCullagh's 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 the 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)^3 k'(r, t) \right] * \left[ \frac{1}{3}\sigma_2^{c0}(r, t) - \sigma_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)^3 k'(r, t) \right] * \left[ \frac{1}{3}\sigma_2^{c0}(r, t) + \sigma_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)^3 k'(r, t) \right] * \sigma_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)^3 k'(r, t) \right] * \sigma_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)^3 k'(r, t) \right] * \sigma_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)^3 k'(r, t) \right] * \sigma_2^{c1}(r, t) dr.
\end{aligned}
\right.
\quad (8) $$

Here a and b are the surface and core radii, respectively, and $\delta I_o(t)$ is the inertia perturbation induced by $\sigma_2^{c0}$, the degree-0 mantle mass anomaly. This term disappears in eq. (7) and consequently does not create TPW; its only effect is to perturb the norm of the rotation vector (through the term $\frac{dH}{dt}$, which has been neglected in our approximation), i.e. the geological length of day.

Inertia changes resulting from mantle density heterogeneities act immediately on the planet through the $\delta(t)$ Dirac function. However, they also have a delayed effect as a result of viscoelastic deformations of the interfaces of our incompressible model (CMB, 670-km-depth discontinuity, surface...). These effects are taken into account through the time convolution with the degree-2 internal loading Love number $k'(r, t)$. This Love number expresses the contribution of viscoelastic deformations induced by internal loads within the Earth to the surface gravitational potential. In the time domain, it may be written

$$ k'(r, t) = k^{\prime E}(r)\delta(t) + \sum_{j=1}^{M} k'_{j}(r)e^{-t/\tau_{j}}H(t), \quad (9) $$

where $k^{\prime E}(r)$ is the elastic internal loading Love number and $k'_{j}(r)$ are the viscous internal loading coefficients.

We define the viscous geoid kernel as the quantity $[k^{\prime E}(r) + \sum_{j=1}^{M} k'_{j}(r)\tau_{j} + (\tau_{a})^{3}]$. This kernel is a good approximation when the timescale of internal loads greatly exceeds the viscoelastic relaxation times ($t \gg \tau_{j}$). Although this is not the case in our computation, this kernel gives useful information about the viscous behaviour of the Earth in terms of geoid and internal loads.

The behaviour of the 670-km-depth discontinuity plays an important role because the buoyancy effect associated with a density jump influences the inertia tensor perturbations. This is evident in the viscous geoid kernels of harmonic degree 2, shown in Fig. 2 for models

<figure><img src="image_4.png" /><figcaption>Figure 2. Viscous geoid kernel for the Earth model described in Table 1, with (left panel) or without (right panel) a density jump at 670 km depth. Kernels are shown for three values of the viscosity ratio between the lower and upper mantle: $\frac{1}{2} = 5$ (dotted line), $\frac{1}{2} = 30$ (solid line) and $\frac{1}{2} = 50$ (dashed line).</figcaption></figure>
with (left) or without (right) density jump at the 670-km-depth discontinuity and for different viscosity ratio $v_3/v_2$ from 5 to 50 between the lower and the upper mantle. Note that the viscous geoid kernel vanishes at 670 km depth for the model with a density jump. Note especially that the most important contributions to the inertia perturbations come from deep mantle density anomalies, an effect that is amplified by the presence of the less viscous D'' layer. In the computation of Fig. 2, we have assumed that the upper mantle and the D'' layer have a viscosity of approximately $10^{21}$ Pa s and that the viscosity of the lithosphere $v^L = 10^{26}$ Pa s. Different values of $v^L$ in the range $10^{22}-10^{26}$ Pa s do not change significantly the viscous geoid kernel curve for mantle density heterogeneity. Hereafter, our five-layer Earth model with $v^L = 10^{26}$ Pa s and $v_3/v_2 = 30$ will be called the reference model.

In the next sections, we will compute the shift of the rotation axis induced by inertia tensor perturbations for degree-2 mantle mass anomalies related to large-scale upwellings.

The most plausible cause for the slow TPW rate of the Earth is probably related to the great stability of subduction zones and the large-scale heterogeneity structure of the mantle. Richards et al. (1997) investigated the influence of the subducted plates on the TPW, using a geodynamic model in which mantle mass anomalies were estimated using plate motion reconstructions together with the assumption that subducted slabs sink vertically into the mantle (Ricard et al. 1993b). They studied the motion of the rotation axis of the Earth relative to hotspots for isoviscous and layered viscosity models, and obtained TPW rates in good agreement with palaeomagnetic observations. However, their present-day computed geoid had non-zero tesseral degree-2 coefficients and consequently their computed present-day maximum principal inertia axis (PIA) did not coincide with the geographic North Pole. Their model did not take into account plumes or domes coming from the CMB; but it is now known that such features (especially the domes responsible for superswells) may create large wavelength anomalies within the mantle and induce TPW (Steinberger & O'Connell 2002). In the next two sections, we investigate in more detail the influence of an ascending mantle plume on the polar motion and then extend these results to study the effects of superswells on TPW.

# 3 INFLUENCE OF A PLUME CROSSING THE MANTLE ON THE ROTATION AXIS: THEORETICAL STUDY

Our aim here is to understand if an upwelling plume crossing the mantle or its appearance at the surface of the Earth may be correlated with a fast excursion of the pole.

## 3.1 Modelling of the plume crossing the mantle

The influence of cold blob (or slab) sinking in the mantle (Ricard et al. 1993a) or hot blob rising in the mantle (Spada et al. 1996) with constant velocity has been investigated using a simple geometrical model. Here, we model a plume crossing the mantle from the D'' layer to the bottom of the lithosphere using the analytical approach proposed by (Bercovici & Mahoney 1994). The plume comprises a spherical head with radius $R$ and a trailing conduit. The time $t = t_o$ is that at which the sphere detaches from the D'' layer (Fig. 3b).

The volumetric flux of material in the conduit $Q$ is assumed to be constant until the plume reaches the bottom of the lithosphere and zero thereafter. The initial radius $R_o$ of the plume head at $t = t_o$ is (Bercovici & Mahoney 1994)

$$R_o = \left( \frac{3\eta_{lm}Q}{4\pi\Delta\rho g} \right)^{1/4}, \qquad (10)$$

where $\eta_{lm}$ is the lower-mantle viscosity, $g$ is the gravitational acceleration and $\Delta\rho$ is the density contrast between the plume fluid and the surrounding mantle.

If the volume of the conduit is negligible, then the radius $R$ of the head increases with time as

$$\frac{4\pi}{3} R^3 = \frac{4\pi}{3} R_o^3 + Q(t - t_o) \Rightarrow R = R_o \left[ 1 + \frac{3Q}{4\pi R_o^3} (t - t_o) \right]^{1/3}. \qquad (11)$$

In the following, we assume that the ascent speed of the hot plume head is the Stokes speed $v = -\frac{2}{9} \frac{\Delta\rho R^2 g}{\eta}$, where $\eta = \eta_{lm}$ if the plume head is in the lower mantle or $\eta = \eta_{um}$ if it is in the upper mantle (the depth where the centre of the spherical plume head is located determines the velocity when the head is crossing the boundary). In models with a chemical interface that prevents passage of the plume head (Liu et al. 1991; Kumagai & Kurita 2000), we assume for simplicity that thermal coupling across the interface induces a mass anomaly that then rises through the upper mantle with the appropriate Stokes velocity.

Let $t_1$ be the time when the top of the plume head reaches the bottom of the lithosphere (Fig. 3c). Because the lithosphere is more viscous than the upper mantle, the plume head stops rising. We assume that $(1-x)$ per cent of the mass anomaly of the plume head spreads under the bottom of the lithosphere to form a spherical cap subtending an angle $\alpha_L$. The remaining $x$ per cent crosses the lithosphere and arrives at the surface of the Earth as another spherical cap subtending an angle $\alpha_S$. The thicknesses $h_L$ and $h_S$ of these two caps increase with time, assuming that the spherical plume head disappears like a setting sun. For $t \ge t_M$ there is no more plume head (Fig. 3d) and the spherical caps have a constant thickness:

$$2\pi \Delta\rho h_L(t_M)(1 - \cos \alpha_L) a_L^2 = \frac{4\pi}{3} \Delta\rho R(t_1)^3 (1 - x),$$

$$2\pi \Delta\rho h_S(t_M)(1 - \cos \alpha_S) a_S^2 = \frac{4\pi}{3} \Delta\rho R(t_1)^3 x.$$
<figure><img src="image_5.png" /><figcaption>Figure 3. Modelled motion of a plume crossing the mantle.</figcaption></figure>

Within the mantle, the density contrast between the fluid plume and the colder surrounding mantle is negative, implying positive (ascending) Stokes velocities. In order to conserve the mass of the Earth, we make the rather arbitrary hypothesis that the rising hot material is compensated by a positive mass anomaly at the base of the D'' layer. We assume that this internal load too is a spherical cap with an angle $\alpha_{CMB}$, a density anomaly $-\Delta\rho$ and a thickness $h_{CMB}$. This thickness increases until $t = t_o$ while the initial sphere forms (Fig. 3a) and then continues to grow as the plume head rises through the mantle according to

$$+ 2\pi \Delta\rho h_{CMB}(t)(1 - \cos \alpha_{CMB})b^2 = \frac{4\pi}{3} \Delta\rho R(t)^3.$$

Because this mass anomaly is located close to the density contrast of the CMB, it will be perfectly compensated on long timescales and therefore will not perturb the inertia tensor (the geoid kernels plotted in Fig. 2 vanish for internal loads located at this interface).

## 3.2 Perturbations of the inertia tensor and true polar wander

The mass anomaly within the mantle comprises four parts: the sphere crossing the mantle, the cap at the CMB, and the two spherical caps at the bottom of the lithosphere and at the surface of the Earth. The inertia tensor perturbations $I_{ij}$ may be computed for our five-layer reference model taking into account the viscoelastic deformations from eq. (6).

As an example, we assume that the hot plume crosses the mantle vertically at latitude $\theta_o = 60^\circ$ and longitude $\varphi_o = 25^\circ$. We compute (Fig. 4) the perturbations in the inertia tensors for $\Delta\rho = -50 \text{ kg m}^{-3}$, $R_o = 150 \text{ km}$ (Bercovici & Mahoney 1994) and $\alpha_S = 20^\circ$ (Courtillot et al. 1999), assuming $\alpha_L = \alpha_S$ and $x = 0.25$. We put $\alpha_{CMB} = 20^\circ$: this value is quite arbitrary but it does not influence the inertia perturbation because this load is located at the CMB and is isostatically compensated. For $\frac{v_3}{v_2} = 30$, we find $t_o = 77.6 \text{ Myr}$; the arrival time of the plume at the base of the lithosphere is $t_1 = 216 \text{ Myr}$ and the final time is $t_M = 217 \text{ Myr}$.

The amplitude of the perturbations is approximately $5 \times 10^{-7} \times I_o$, where $I_o$ is the inertia moment of the non-rotating spherical Earth. The inertia perturbations are even smaller for a mantle model with a density jump at 670 km depth as a result of the buoyancy effect at this interface.

We have tested two extreme values for the lithospheric viscosity ($v^L = 10^{26} \text{ Pa s}$, i.e. a quasi-elastic lithosphere and $v^L = 10^{22} \text{ Pa s}$), but this parameter does not change the results significantly (Fig. 4).

If the Earth is in hydrostatic equilibrium, the maximum PIA and the rotational axis will coincide; consequently, if inertia products appear, the rotation axis will shift towards the new maximum PIA. It is thus interesting to investigate the PIA associated with the mantle mass anomalies induced by the temporal evolution of the upwelling plume.

For this system, the minimum principal axis can easily be shown to coincide with the axis of the plume (here $60^\circ \text{N}$ and $25^\circ \text{E}$) and be associated with a large negative eigenvalue. The negative sign is a consequence of the viscoelastic deformations of the interfaces; the eigenvalue would be positive if these deformations were neglected by setting $k'(r, t) = 0$, because the plume has a negative density jump relative to the surrounding mantle. The two other PIA are in the plane perpendicular to this axis and are associated with a positive eigenvalue with a multiplicity of 2: all axes located in this plane are maximum PIA. Consequently, the pole of rotation drifts from its initial position towards this plane, following a great circle on the sphere (i.e. taking the shortest way). The trajectory of the pole will therefore depend strongly on its initial position.

The inertia tensor perturbation induced by slabs sinking in the mantle is approximately $10^{-5} \times I_o$ (Richards et al. 1997), or two orders of magnitude larger than that induced by a plume. Consequently, the amplitude of the term $\frac{I_{kl}}{\alpha I_o T_1}$ in eq. (7) is large and then the rotational bulge is quasi-instantaneously perpendicular to the maximum PIA. In our study the amplitude of $I_{kl}$ is weaker (approximately $2 \times 10^{-7} \times I_o$) and consequently the rotational axis will shift more slowly towards the maximum PIA.

We compute the perturbation of the rotation vector induced by the inertia tensor perturbations, using eq. (7). The latitude, longitude and velocity of the rotational axis are plotted in Fig. 5 as a function of time. Note that the pole tends to the direction ($30^\circ \text{N}, 205^\circ \text{E}$), i.e. to
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Inertia tensor perturbations [normalized to the spherical moment ($I_o = 0.33 \text{ Ma}^2$)] induced by the model plume described in the text, with (left) or without (right) a density jump at 670 km depth. Results are shown for two values of the lithosphere viscosity: $10^{26}$ Pa s (top) and $10^{22}$ Pa s (bottom); (b) eigenvalue associated with the axis of the plume.</figcaption></figure>

the maximum PIA, which is perpendicular to the axis of the plume, and that the velocity of the polar wander is higher for the mantle model without a density jump because of the larger amplitude of the inertia tensor perturbations. In any case, the amplitude of wander is small.

In our computation, the lithosphere is assumed to be viscoelastic, so that the loads beneath the lithosphere and at the surface of the Earth are isostatically compensated on scale of millions of years. If we assume instead a global purely elastic shell at the top of the lithosphere of thickness 10–40 km (Mound *et al.* 2003), the density anomalies beneath the lithosphere and at the surface of the Earth will be only partly
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Perturbation of the rotational axis induced by the viscoelastic inertia tensor perturbations plotted in Fig. 4(a), for lithosphere viscosity $10^{26}$ Pa s (solid line) and $10^{22}$ Pa s (dashed-dotted line). Shown are the temporal evolutions of the longitude $\phi$ (top), latitude $\theta$ (centre) and the velocity $\sqrt{\dot{\theta}^2 + \cos^2 \theta \dot{\phi}^2}$ (bottom), where $\dot{\phi}$ denotes the time derivative.</figcaption></figure>

compensated: after the viscoelastic relaxation, a constant term will remain in the inertia tensor perturbations for $t > 217$ Myr. Nevertheless, the amplitude of this inertia perturbation will be weak and the associated shift of the rotational axis that tends to coincide with the final PIA will be very slow.

To conclude this simple example, we have investigated the influence of different parameters on the polar wander, assuming a homogeneous mantle. Fig. 6(a) shows the velocity of the polar wander for different values of the lower-mantle viscosity from $5 \times 10^{21}$ up to $5 \times 10^{22}$ Pa s: note that an increase in the lower-mantle viscosity strongly decreases the rate of polar wander. We also varied $x$, but because of the quasi-isostatic compensation of the loads located at the bottom and at the top of the lithosphere, this parameter does not perturb significantly the TPW.

Finally we have investigated, for our reference model, the influence of the radius of the plume head. In our previous computation, $R$ increased because of the constant flux $Q$ in the conduit and varied from 150 km at the top of the D'' layer up to more than 300 km at the top of the upper mantle. Assuming that $R$ remains constant and equal to $R_o$ produces two effects on the inertia tensor perturbations (Fig. 6b): it decreases the Stokes velocity, thereby increasing the transit time $t_M$ across the mantle by approximately a factor of 2, and it decreases the amplitude of the inertia perturbations (because the degree-2 mass anomaly associated with the plume head is weaker). The resulting polar wander velocity is plotted in Fig. 6(c), for $v_3 = 3 \times 10^{22}$ Pa s.

A more realistic modellization would be to treat the effect of the plume in the presence of other, and much larger, inertia tensor perturbations (resulting from slabs, as explained at the end of the next section), which would cause plume-related TPW to be even much slower than computed here.

We consequently conclude that the amplitude of the polar wander induced by an upwelling plume is weak. It would therefore seem difficult to correlate a fast excursion of the TPW with a plume crossing the mantle or with the appearance of a hotspot at the surface of the Earth.

## 4 SUPERSWELLS AND TPW

We now extend the results obtained above for a plume to investigate the effects of two superplumes located beneath Africa and Polynesia. A recent 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 periodicities ranging from 100 Myr up to 1 Gyr. Such domes could be responsible for the superswells observed at the surface of the Earth.

In this preliminary computation, we model the domes as hemispheres (Fig. 7) and use the results of the previous section. From tomography studies (e.g. Romanowicz & Gung 2002), we assume that the initial radius $R_o$ is 1600 km. An angle $\alpha_{CMB} = 65^\circ$ of the spherical cap at the CMB is necessary to compensate for the negative mass of the dome within a cap whose thickness does not exceed that of the D'' layer. We
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Velocity of the shift of the rotational axis at the surface of the Earth induced by a plume crossing the homogeneous mantle for different values of the viscosity ratio between the lower and upper mantle: $\frac{v_3}{v_2} = 5$ (dotted line), $\frac{v_3}{v_2} = 30$ (solid line) and $\frac{v_3}{v_2} = 50$ (dashed line). (b) Inertia tensor perturbations assuming a plume head with a constant radius, for our reference model ($\frac{v_3}{v_2} = 30$) without the density jump at 670 km depth. (c) Polar wander induced by the inertia tensor perturbations plotted in (b).</figcaption></figure>

assume a density contrast between the plume and the surrounding mantle of –50 kg m⁻³, based on the conclusion of Davaille (1999) that the doming regime should occur for density contrasts less than approximately 1 per cent. We stop the rise of both domes at the top of the lower mantle after 50 Myr. The first dome is located beneath Africa at latitude and longitude (θ₁, φ₁) and the second beneath Polynesia at (θ₂, φ₂). The dome beneath Africa starts to oscillate 30 Myr before the other.

The computation of the PIA of this system is not trivial. It is easy to determine the minimum principal axis, which is associated with a large negative eigenvalue: it corresponds approximately to the axis where the domes, which are quasi(but not exactly)-antipodal, are located. The maximum principal axis will be in a plane perpendicular to this axis. The two others eigenvalues are positive and very close; consequently the maximum PIA is strongly dependent on the position of the domes.

In our test, we first assume for simplicity that the two domes are on the equator, i.e. θ₁ = θ₂ = 0° and φ₁ = 30°; φ₂ = 230°. In this case, the maximum PIA is close to the north axis. Fig. 8 shows the temporal evolution of the inertia tensor perturbations and the minimum principal axis for the mantle model introduced in Section 3. The perturbations of the inertia tensor of the Earth oscillate with a period approximately 100 Myr, the period of the domes themselves. Note that the minimum PIA is in the equatorial plane and oscillates in longitude between 30°E and 50°E.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Geometrical model of two oscillating domes. The domes are modelled as two hemispheres oscillating vertically within the mantle with a periodicity of 100 Myr. The first dome is located beneath Africa and starts to oscillate 30 Ma before the other, located beneath Polynesia. The density contrast between the plume and the surrounding mantle is −50 kg m⁻³.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Inertia tensor perturbations induced by the simple model of equatorial domes described in the text, for mantle models with (left) or without (right) a density jump at 670 km depth: (a) components of the tensor; (b) minimal principal inertia axis (PIA).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Temporal evolution of the direction of the rotation axis using an orthographic projection (top) and velocity of this polar wander (bottom), for mantle models with (left) or without (right) a density jump at 670 km depth and different initial conditions for the rotational pole: $\theta_o = 50°$, $\phi_o = 150°$ (filled circle and solid line) and $\theta_o = 25°$, $\phi_o = 250°$ (star and dashed line). The time step is 10 Myr. The two grey semi-spheres represent the location of the centre of the domes: one beneath Africa and one beneath Polynesia.</figcaption></figure>

Fig. 9 shows the associated temporal evolution of the direction of the rotation axis and its velocity, for two different initial conditions ($\theta_o, \phi_o$) for the position of the rotational pole at $t = 0$ and for mantle models with and without a density jump at 670 km depth. Note that the pole first shifts along a great circle at increasing speed toward the plane, which bisects the starting point locations of the two superswells. The pole next heads to its final position more or less along the meridian, which contains the maximum PIA, with oscillations both in velocity and in trajectory (slight oscillations). The trajectory of the pole is consequently dependent on the initial position of the rotational axis. The fast original polar motion shown here appears to be caused by this arbitrary starting position. These initial values have to be compared to the PIA induced by subducted plates.

The temporal evolution of the velocity of this polar wander shows that there are alternating periods of slow and faster TPW.

Comparing the left and the right column of Figs 8 and 9, we show that the amplitude of the inertia perturbations is smaller for the mantle models with a density jump at 670 km depth, because of the associated buoyancy effect at this interface. The rotational axis will consequently tend more rapidly (in less than 200 Myr) toward the maximum PIA (here the North Pole) for the model without a density jump. This explains the higher velocity of the polar wander during the first 200 Myr for this model.

To investigate a more realistic model, we position the domes at $\theta_1 = -15°$, $\phi_1 = 30°$ and $\theta_2 = -5°$, $\phi_2 = 230°$. In this case, the minimum PIA oscillates around the positions (–15°N, 30°E) and (5°N, 50°E) and the maximum PIA is in a plane perpendicular to this direction. The final position of the rotational axis is 45°N, 315°E.

Comparing these results with the ones of Figs 8 and 9, we see that both the maximum PIA and the TPW are strongly dependent on the initial position of the domes, always adjoining the equator. To decrease this strong dependence, it is necessary to consider at the same time the effect of subducted plates sinking into the mantle. The latter will break the symmetry of our system and thereby strengthen the PIA; they will also influence the trajectory of the rotational axis through the initial values of the pole of rotation.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** (a) PIA induced by the inertia tensor perturbation resulting from the equatorial domes shown in Fig. 8 superimposed on the additional inertia perturbation related to the subducted plates, for mantle models without a density jump at 670 km depth. (b) Temporal evolution of the direction of the rotation axis, using an orthographic projection assuming that the initial rotational axis is aligned with the PIA related to subductions, induced by both subductions and domes (filled circle and solid line) and by domes only (stars and dashed line). The time step is 10 Myr. The two grey semi-spheres represent the location of the centre of the domes: one beneath Africa and one beneath Polynesia. (c) Velocity of the TPW described in (b).</figcaption></figure>

This effect can be quantified here already very easily. We do not model anything about subduction, we just assume that there is a constant additional inertia tensor contribution. As an example, we assume an inertia tensor perturbation $I_{ij}^S$ related to subduction:

$$ \frac{I^S}{I_o} = \begin{pmatrix} -0.158604 \times 10^{-4} & 0.771024 \times 10^{-6} & -0.698933 \times 10^{-5} \\ 0.771024 \times 10^{-6} & -0.158604 \times 10^{-4} & -0.524940 \times 10^{-5} \\ -0.698933 \times 10^{-5} & -0.524940 \times 10^{-5} & 0.531581 \times 10^{-5} \end{pmatrix}. \qquad (12) $$

This value is quite arbitrary but plausible: it corresponds to the inertia perturbation induced, 95 Ma, by the mantle density heterogeneities associated with the subducted plates sinking into the mantle derived by Ricard *et al.* (1993b). The PIA related to this mass distribution is (59°N, 247°E).

We then investigate the PIA induced by both equatorial domes (plotted in Fig. 8) and subductions: there is now a larger positive eigenvalue associated with the maximum PIA axis, which oscillates around the position (80°N, 230°E; see Fig. 10a).

We next compute the TPW induced by both domes and subductions, starting with the rotation axis aligned with the PIA related to subductions. We plot the resulting temporal evolution of the rotational axis (Fig. 10b) and its velocity (Fig. 10c), assuming a homogeneous mantle. The comparison between the curve corresponding to the TPW when subductions and domes are simultaneously taken into account (filled circle) and the curve where only the domes are taken into account (start) is very interesting: first, we see that the amplitude of the initial fast TPW velocity is decreased and, second, we see that the rotational axis, similarly to the PIA, oscillates around the position (80°N, 230°E);
the velocity of the TPW is periodic with a period of 100 Myr and an amplitude approximately 0.6° Myr⁻¹. The damped effect appearing in Fig. 9(b) disappears here because the maximum PIA has gained in strength.

Our conclusion from this preliminary study is that a doming regime within the mantle involves significant TPW and that some of the features of this TPW that are directly linked to the periodicity of doming are reminiscent of observed phases of slow and fast TPW, with similar peak velocities.

## 5 CONCLUSION

In this paper, we have investigated the influence of mantle plumes with different spatial wavelengths on the temporal variation of the rotational axis over geological timescales. TPW depends on two terms: the relaxation of the rotational bulge and the amplitude and temporal evolution of the inertia tensor perturbations induced by internal loads. We found that an upwelling plume whose head has a radius of 300 km induces very slow wander of the rotational axis. Finally, we showed that a doming regime within the mantle involves significant TPW and that it is necessary to add this effect to that induced by subducted plates before comparing the computed TPW with palaeomagnetic observations.

## ACKNOWLEDGMENTS

This study was supported by a CNRS-INSU (Dyeti) grant and is IPGP contribution number 2002. The author thanks Jean Besse, Vincent Courtillot, Anne Davaille and Neil Ribe for discussions and remarks on the original manuscript. The author gratefully acknowledges B. Steinberger and G. Spada for their helpful reviews.

## REFERENCES

* Bercovici, D. & Mahoney, J., 1994. Double Flood Basalts and Plume Head Separation at the 660-kilometer Discontinuity, *Science*, **266**, 1367–1369.

* Besse, J. & Courtillot, V., 1991. Revised and synthetic apparent polar wander paths of the African, Eurasian, North American and Indian plates, and true polar wander since 200 Ma, *J. geophys. Res.*, **96**, 4029–4050.

* 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 (BC02), 2300, doi: 10.1029/2000JB000050.

* Courtillot, V., 2004. True Polar Wander, in *Kluwer Encyclopedia of Geomagnetism and Paleomagnetism*, eds Gubbins, D. & Herrero-Bervera, E., submitted.

* Courtillot, V. & Renne, P.R., 2003. On the ages of flood basalt events, *Comptes Rendus Geoscience*, **335**(1), 113–140.

* Courtillot, V., Manighetti, I., Jaupart, C. & Besse, J., 1999. On causal links between flood basalts and continental break-up, *Earth planet. Sci. Lett.*, **166**, 177–195.

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

* Kumagai, I. & Kurita, K., 2000. On the fate of mantle plumes at density interfaces, *Earth planet. Sci. Lett.*, **179**, 63–71.

* Lefftz, M., 1991. Aspects théoriques de la rotation de la Terre et de son noyau: influence de la viscosité, *PhD thesis*, Université de Strasbourg, p. 261.

* Lefftz, M., Legros, M. & Hinderer, J., 1991. Non-linear equations for the rotation of a viscoelastic planet taking into account the influence of a liquid core, *Celestials Mechanics and Dynamical Astronomy*, **52**, 14–43.

* Liu, M., Yuen, D.A., Zhao, W. & Honda, S., 1991. Development of diapiric structures in the upper mantle due to phase transition, *Science*, **252**, 1836–1839.

* Mound, J.E., Mitrovica, J.X. & Forte, A., 2003. The equilibrium form of a rotating Earth with an elastic shell, *Geophys. J. Int.*, **152**, 237–241.

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

* Peltier, W.R., 1974. Impulse response of a Maxwell Earth, *Rev. Geophys. Space Physics*, **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. geophys. Res.*, **B98**(12), 21895–21909.

* 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.

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

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

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

* Spada, G., Ricard, Y. & Sabadini, R., 1992. Excitation of true polar wander by subduction, *Nature*, **360**, 452–454.

* Spada, G., Sabadini, R. & Boschi, E., 1996. The spin and inertia of Venus, *Geophys. Res. Lett.*, **23**(15), 1997–2000.

* 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. & O'Connell, R., 2002. The convective mantle flow signal in rates of True Polar Wander, in *Ice Sheets, Sea Level and the Dynamic Earth*, Geodynamics Series 29, pp. 233–256, eds Mitrovica, J.X. & Vermeersen, B.L.A., American Geophysical Union, Washington, DC.