# Global plate motion and the geoid: a physical model

Y. Ricard, C. Froidevaux, L. Fleitout

► To cite this version:

Y. Ricard, C. Froidevaux, L. Fleitout. Global plate motion and the geoid: a physical model. Geophysical Journal International, 1988, 93 (3), pp.477-484. 10.1111/j.1365-246X.1988.tb03875.x. hal-02046609

HAL Id: hal-02046609

https://hal.science/hal-02046609

Submitted on 20 Aug 2021

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
# Global plate motion and the geoid: a physical model

Y. Ricard, C. Froidevaux and L. Fleitout

Département de Géologie, Ecole Normale Supérieure, 24 rue Lhomond, 75231 Paris, France

Accepted 1987 November 12. Received 1987 November 9; in original form 1987 July 31

## SUMMARY

The effect of lateral variations in the sub-lithospheric viscosity, i.e. in the lithosphere-asthenosphere coupling, is incorporated in dynamic earth models driven by plate velocities. If the coupling is stronger under fast-moving continents, it is shown that a geoid low will develop in their wake and a high in front of them. Thus the well-known low in the Indian ocean could be explained, at least in part, in terms of induced upper-mantle dynamics. Similarly, if the thickness of oceanic plates increases with age the models show that trenches should be associated with marked geoid highs, whereas ridges could correspond to much weaker geoid lows. This also seems to agree with some features of the observed geoid at very long wavelengths. The mathematical framework of such dynamic earth models is developed extensively here. These models are characterized by lateral viscosity variations inside an outer shell and a purely radial viscosity structure at greater depth. Their internal flow patterns are only driven by imposed surface velocities, not by internal loads.

**Key words:** geoid, plate motion, mantle dynamics

## INTRODUCTION

The long wavelength geoid ($λ > 4000$ km) has been well known over both continents and oceans for more than two decades. Its interpretation remains controversial. No obvious correlation is found with the present distribution of continents and oceans or with the geometric pattern of plate tectonics. The now emerging seismic tomography structure of the Earth's mantle (Dziewonski 1984; Nataf et al. 1986) has revealed the existence of deep lateral velocity heterogeneities, resolved for the first six to eight degrees of spherical harmonics. The corresponding density variations can create sizeable geoid undulations. Indeed it has been shown that these deep masses will induce internal flows and deflect the various density interfaces within the Earth, including of course its outer surface (Hager, Ricard et al. 1984). The resulting dynamic configuration generates a geoid which correlates significantly with observations, in particular for degrees 2 or 3 (Ricard et al., Hager et al. 1985). The oceanic lithospheric slabs dipping into the Earth's interior around the Pacific Ocean give rise to a well defined belt of dense masses, not easy to resolve by present-day tomographic techniques. This deep pattern correlates strikingly with the geoid of degree $l \ge 4$ (Hager 1984).

Another meaningful correlation is found between the observed geoid and the divergence of the velocity field defined by plate tectonics, in particular for degree 4 (Forte & Peltier 1987). This observation was also predicted by dynamic earth models (Ricard et al. 1984). Indeed the surface kinematics requires a return flow driven by horizontal pressure gradients. The latter support an uncompensated topography which produces a sizeable

contribution to the geoid. One thus expects geoid highs across the zones of convergence and lows across the zones of opening. The maxima are well marked around the Pacific, whereas the minima corresponding to ocean ridges are less conspicuous (Jacoby & Seidler 1981).

As far as physical mechanism are concerned, two approaches have thus been tried. The first starts from the observed density heterogeneities, which are taken as sources driving the internal dynamics, and hence the surface flow. Up to now this first class of models was built on the assumption of a purely radial variation of the mantle viscosity. Therefore, the predicted surface velocities are unlikely to match real plate velocities which have equal strengths in poloidal and toroidal components (Hager & O'Connell 1978). These models only yield a poloidal flow, i.e. with zones of divergence and convergence, but no shear. Having no lithospheric plates, they prevent small-size masses, like slabs or compensated ridge topographies, from driving large-scale surface displacements. Nevertheless this first approach partly explains the existing geoid.

The second approach takes the observed plate velocities as its starting point. It imposes these boundary conditions to the computed internal flow field, and may or may not include known internal density variations. Up to now this second class of models has also assumed a purely radial dependence for the viscosity. Although the driving forces are *ad hoc*, these models can simulate large-scale return flows and their corresponding geoid signatures.

Ideally, one would like to construct a model predicting both the plate velocities and the total geoid. It would necessarily include lateral viscosity variations, in particular, at shallow depths. This paper will address this general
problem quantitatively. A mathematical formulation will be worked out for an earth model with a simple geometry. Two idealized applications with dynamically induced geoid will then be presented. One will show the resulting geoid near a continent-ocean passive margin. The other will deal with ridges and trenches.

## THE MATHEMATICAL FORMULATION

A lithospheric plate is subjected to both driving and resistive forces. The coupling between the plate and the underlying asthenosphere is generally believed to belong to the second group (Richardson *et al.* 1979). Another view is that the plates are moving passively, driven by the underlying mantle convection (Jarvis & Peltier 1982). In this case the uniform velocity imposed by a rigid plate could lead to a situation where the mantle drives the lithosphere beneath the plate interior and is driven by that plate near its edges. If so, only this last portion would indeed be subjected to a resistive drag. In any case the coupling is likely to be larger under old continents than under young oceans (Forsyth & Uyeda 1975; Froidevaux *et al.* 1977).

Our model earth starts from the hypothesis that the plates are mainly driven by mass heterogeneities near their boundaries giving rise to slab pull and ridge push. Thus our model lithosphere drags the asthenosphere. Our earth geometry comprises one central sphere where the mechanical properties vary only radially, and an outer shell of thickness L in which the viscosity varies laterally. This outer shell includes a lithosphere of viscosity η₀ and of thickness H(θ, φ) where θ is the colatitude and φ the longitude. The remaining portion of thickness L - H(θ, φ) has an asthenospheric viscosity η₁. Fig. 1 depicts this configuration.

What is new here is the existence of lateral viscosity variations within the outer shell. For the inner sphere with a purely radial viscosity profile, the dynamics is governed by known equations which are most suitably expressed in terms of spherical harmonics $Y^{lm}(\theta, \varphi)$ (Kaula 1975; Hager & O'Connell 1978; Ricard *et al.* 1984). These are defined as follows:

$$Y^{lm}(\theta, \varphi) = P^{lm|l}(\cos \theta) \exp(im\varphi) \quad (1)$$

where $P^{lm|l}$ is a Legendre function. It is more efficient and elegant to introduce generalized spherical harmonics $Y^{lm+}(\theta, \varphi)$ (Edmonds 1960; Phinney & Burridge 1973) for describing vectorial quantities such as the horizontal velocity and the horizontal drag, and keep $Y^{lm}(\theta, \varphi)$ for scalar quantities such as the radial velocity, the radial stress or the geoid height. We have the following definition:

$$Y^{lm+}(\theta, \varphi) = \frac{i}{\sqrt{l(l+1)}} \left( \frac{1}{\sin \theta} \frac{\partial Y^{lm}(\theta, \varphi)}{\partial \varphi} + i \frac{\partial Y^{lm}(\theta, \varphi)}{\partial \theta} \right). \quad (2)$$

If horizontal and vertical velocities are imposed on the surface of the inner sphere of our earth model, which has no internal loads, the induced flow creates stresses and deflects the various interfaces. These deflections generate geoid undulations of amplitude N if the various layers have different densities. Rewriting the linear relationships derived

<figure><img src="image_2.png" /><figcaption>Figure 1. Schematic diagram of the model earth comprising an outer shell of thickness L and an inner sphere made up of mantle and core material. The inner sphere has a viscosity which varies only radially. In numerical applications presented in the next section, this radial variation will simply consist of a jump at 650 km. The viscosity in the outer shell is η₀ down to a depth H(θ, φ) which depends upon geographical position. In the lower part of this shell the viscosity decreases to a value η₁. One can thus think of the depth H(θ, φ) as being the interface between lithosphere and asthenosphere.</figcaption></figure>

in the papers we just quoted gives:

$$(\tau_{+}^{L}) = M_{1}(V_{+}^{L}) + M_{2}(V_{r}) \quad (3)$$

$$(N) = M_{3}(V_{+}^{L}) + M_{4}(V_{r}). \quad (4)$$

All of the above quantities are clearly defined in the mathematical appendix. $(\tau_{+}^{L})$ and $(V_{+}^{L})$ are vectors representing all the components of the drag and horizontal velocity at the surface of our inner sphere. They are expressed on the $Y^{lm+}$ basis. $(N)$ and $(V_{r}^{L})$ are also vectors made up of the components of the geoid and of the radial velocity expressed on the basis of ordinary spherical harmonics. The index L refers to the depth L of our earth model in Fig. 1, i.e. to the surface of the inner sphere. The four matrices in (3) and (4) depend both upon the radial viscosity profile of the mantle, and upon the degree l of the applied surface velocity field.

Equation (3) relates three unknown physical quantities which are continuous at the interface of depth L between the inner sphere and the outer shell of our model earth. The dynamics of this outer shell will yield two more equations relating the same quantities to the imposed velocities at zero depth. Thus all quantities in (3) and (4), in particular the geoid, will be determined once the surface kinematics has been given. We now construct these additional mathematical relationships.

For the inhomogeneous outer shell, the imposed horizontal surface velocity $V^O$ is a function of θ and φ, as well as the unknown horizontal velocity $V^L$ and shear $\tau^L$ at the interface with the inner sphere. A classical expression relates the horizontal shear stress $\tau$ and velocity $V$:

$$\frac{\partial V}{\partial r} + \nabla_2 v_r - \frac{V}{r} = \frac{\tau}{\eta}, \quad (5)$$
where $v_r$ is the radial component of the total velocity, $\nabla_2$ the horizontal component of the gradient and $r$ the distance to the Earth's centre. In the thin-shell approximation where the shell thickness $L$ is negligible compared to both the Earth radius $R$ and the horizontal wavelength, the first term dominates the left-hand side of the above equation, which yields after integration:

$$ \nabla^O - \nabla^L = L \left\langle \frac{1}{\eta} \right\rangle \tau^L. \qquad (6) $$

The equation of mass conservation in the outer shell can also be integrated through the thickness $L$. Assuming a vanishing radial velocity at zero depth, this gives:

$$ v_r^L = L \nabla_2 \langle \nabla \rangle. \qquad (7) $$

The quantities between brackets in (6) and (7) represent vertical averages. Considering Fig. 1, we can formulate them in terms of the functionals $H(\theta, \varphi)$, $\nabla^O$ and $\nabla^L$ and of the two viscosities $\eta_0$ and $\eta_1$:

$$ \left\langle \frac{1}{\eta} \right\rangle = \frac{1}{\eta_0 \eta_1 L} (\eta_1 H + \eta_0 (L - H)) \qquad (8) $$

$$ \langle \nabla \rangle = c \nabla^O + (1-c) \nabla^L \qquad (9) $$

with

$$ c = \frac{1}{2} + \frac{H(L-H)}{L^2} \frac{\eta_0 - \eta_1}{2\eta_0\eta_1 \langle 1/\eta \rangle}. \qquad (10) $$

On the basis of ordinary and generalized spherical harmonics the local equations (6) and (7) become linear spectral relationships which are derived in the Appendix. They read:

$$ (\nabla_L^O) - (\nabla_+^L) = M_5(\tau_+^L) \qquad (11) $$

$$ (v_r^L) = M_6(V_+^O) + M_7(V_+^L). \qquad (12) $$

The symbols have the same meaning as in the equations (3) and (4) derived for the inner sphere. Here the index $O$ refers to the outer surface of the model earth and the matrices $M_5$, $M_6$ and $M_7$ reflect the coupling properties of the thin shell.

In conclusion we use equations (3), (11) and (12) to determine $(v_r^L)$, $(V_+^L)$ and $(\tau_+^L)$, i.e. the velocity and shear stress pattern at the interface between the outer shell and the inner sphere of our earth model. The geoid can then be computed on the basis of equation (4). The velocity and stress fields below this interface can also be determined, if necessary, following standard procedures (see e.g. Ricard et al. 1984).

One disadvantage of the proposed scheme lies in the rapidly increasing memory space required in numerical calculations. Indeed the number of elements in the matrices $M_1$ to $M_7$ is proportional to $l_{\max}^4$. This leads us to choose a low cut-off at $l = 14$ for the applications presented below.

One should realize that the presence of lateral variations in the viscosity leads to non-linear equations (see Appendix). For example the expression of (2) on the basis of spherical functions has shown that a stress component of degree $l_3$ interacting with a lateral viscosity distribution of degree $l_2$ will generate a velocity pattern of degrees between $|l_3 - l_2|$ and $(l_3 + l_2)$. In other words, short-wavelength excitations may drive long-wavelength responses, unless

the power spectra of the coupling functions $\langle 1/\eta \rangle$ and $C$ in (2) and (5) decrease rapidly enough with $l$. Because of these non-linear effects, in the examples given in the next section the cut-off at $l = 14$ seems to allow numerically stable predictions only up to degree $l = 12$.

## TWO HEURISTIC EXAMPLES

A moving lithospheric plate will drag the underlying mantle material more efficiently beneath its continental than beneath its oceanic domain. A model where the continent has simply a circular shape will first help us to deduce in which way the geoid may be affected near leading and trailing continent-ocean edges.

As a first step we consider a flat channel of thickness $d$ for a fluid with viscosity $\eta$. If a velocity $v_c$ is imposed over a circular area of radius $R_c$ a bulge of fluid will form in the front and a depression in the wake of the moving circle. These dynamic effects can easily be expressed analytically, without using the rather involved mathematics presented in the first section. The Navier-Stokes equation for an incompressible fluid without lateral density variations implies that $\nabla^2 p = 0$. In polar coordinates $r_c$ and $\theta_c$ with origin at the centre of the moving continental circle, the pressure reads therefore either

$$ p_{\text{in}} = A r_c \cos \theta_c \qquad (13) $$

or

$$ p_{\text{out}} = A \frac{R_c^3}{r_c^2} \cos \theta_c \qquad (14) $$

depending on whether the fluid volume beneath or outside the circle is considered. Here $\theta_c$ is the angle with respect to the orientation of the velocity. Thus under the continent the pressure $p_{\text{in}}$ increases linearly with $x = r_c \cos \theta_c$, the distance along the velocity direction. For a shallow fluid channel, $d \ll R_c$, the corresponding Couette flow defines the value of the above arbitrary constant:

$$ A = 6\eta v_c / d^2. \qquad (15) $$

The corresponding vertical deflection of the fluid surface can therefore be written as

$$ h_{\text{in}} = \frac{6\eta v_c}{\rho g d^2} r_c \cos \theta_c \qquad (16) $$

and

$$ h_{\text{out}} = \frac{6\eta v_c R_c^3}{\rho g d^2 r_c^2} \cos \theta_c. \qquad (17) $$

It represents a dipolar shape parallel to the velocity with a maximum at the front edge and a minimum at the wake edge. The associated pressure drives a return flow within the fluid channel. An uncompensated topography $h$ of characteristic wavelength $\lambda$ generates a geoid anomaly of amplitude

$$ N = (G \rho \lambda / g) h, \qquad (18) $$

where $G$ is the universal gravitational constant. Using this crude estimate and taking $\lambda = 4R_c$ with $R_c = 2000$ km, $\eta = 5.10^{19}$ Pa s, $d = 650$ km and $v = 5$ cm yr$^{-1}$, the above approximate solutions yield a topographic amplitude $h_c(R_c, 0) = 82$ m and a geoid undulation of amplitude
12.5 m. The topography effect could hardly be detectable, but the geoid signature would certainly stand out.

The same circular continental configuration, but on a spherical earth, has been examined in the framework of the mathematical theory developed in the previous section. Now a velocity $v = v_0 \sin \theta$ which depends only on the colatitude $\theta$ is imposed everywhere on the sphere surface. This poloidal velocity field is of degree one. For a radially symmetrical earth model, it does not generate a geoid undulation; it merely displaces the centre of gravity. This is not true in the presence of lateral viscosity variation. Let us take an outer shell of thickness $L = 200$ km, with a continent of radius $R_c = 2000$ km in equatorial position. The lithosphere of viscosity $\eta_0 = 5.10^{22}$ Pa s has a thickness $H(\theta, \varphi)$ of 200 km under the continent, but vanishes under the remaining oceanic Earth surface. Everywhere else the upper mantle has a viscosity $\eta_1 = 5.10^{20}$ Pa s down to the

interface with the lower mantle at 650 km depth. There, the viscosity increases by a factor of 30. The Earth's core is defined as inviscid.

Imposing a northward equatorial velocity $v_0 = 5$ cm yr$^{-1}$, we computed the numerical solutions for the velocities and stresses on the basis of equations (3), (4), (11) and (12) up to degrees and order 14. Fig. 2 depicts the predicted geoid by means of isolines and also in 3-D view. We can recognize the expected dipolar pattern. The amplitude amounts to 20.6 m. This value is lower than the simple analytical estimate, as the flow is now less confined in depth, despite the larger viscosity. A radial viscosity structure which would tend to confine the return flow at shallower depths reduces the lateral extent of this perturbation. It will also enhance the amplitude near both extrema. If all viscosity values are multiplied by one given factor, the geoid will of course be enhanced by the same factor. As in the simple Couette flow

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Computed dipolar geoid pattern induced by a northward-moving continent. The shaded circle represents an ideal continental area located on the equator with a radius of 2000 km and a 200 km thick lithospheric root. The whole surface moves northward, with a velocity of $5$ cm yr$^{-1}$ at the equator. The mantle viscosity increases by a factor 30 at 650 km depth. The upper-mantle viscosity beneath the lithosphere is $5.10^{20}$ Pa s. There is no lithosphere outside the circular continental area. The spherical surface is mapped through a linear representation of both longitude and latitude. On the top map the geoid isolines are drawn every 5 m. On the N edge of the continent the geoid maximum reaches 20.6 m, whereas the S edge is marked by a minimum of similar amplitude. The bottom diagram is a 3-D representation of the same geoid pattern.</figcaption></figure>
model described above, the computed geoid is also proportional to the imposed velocity.

This first example is instructive in the sense that the striking geoid minima observed on the Indian plate are precisely located in the wake of the fast-moving continental masses, India and Australia. The same is true in the wake of North and South America, i.e. along their eastern margins (Forsyth *et al.* 1980). These puzzling features are broadly explained here. Is the predicted maximum in front of the moving continent also observed? The answer is positive for Australia, although the very large amplitude of the existing maximum must also derive from other sources, such as lower-mantle density anomalies (Richards & Hager 1986). The situation is equally favourable north of India, but the thick crust of Tibet can also at least in part explain the existence of a local maximum (Froidevaux & Ricard 1987).

Our second idealized example will address the problem of

the zones of opening, the oceanic ridges and the zones of convergence or subduction. It is worth remembering that the flow within our model earth is merely excited by the imposed velocities at the surface, not by internal loads. Thus, there is no dense subducting slab nor light ridge mantle material contributing to the computed geoid. In other words our dynamic models do not include the associated ridge-push effect, nor are they capable of predicting the well-known 15 m of geoid high over the Atlantic ridge (Sandwell & Schubert 1980).

Let us consider an earth model with two sectorial rigid plates of equal dimension separated by one ridge and one trench, running from N to S. On the equator, they diverge or converge with relative velocities amounting to 10 cm yr⁻¹. Two different viscosity structures will be tested: One with a constant lithospheric thickness H(θ, φ) = 100 km, the other where this value varies as the square root

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

**Figure 3.** Computed geoid induced by the diverging and converging velocity pattern of the earth model with two rigid hemispheric plates. The lithospheric thickness is uniform, 100 km. The diverging and converging velocities are 10 cm yr⁻¹ at the equator. The viscosity structure beneath the lithosphere is identical to that of the model shown in Fig. 2. The map representation also remains the same. The 90° and 270° meridians correspond to a ridge (diverging arrows) and to a subduction zone (converging arrows). The top representation depicts isolines 10 m apart. Thus the geoid maximum and minimum located at the converging and diverging plate boundaries amount to ±63.6 m. The pattern is fully antisymmetrical, as can also be seen in the 3-D block diagram.
of age, reaching the thickness of 200 km at the trenches. In both situations, the lithospheric viscosity amounts to $5.10^{22}$ Pa s, and the mantle viscosity beneath that lithosphere takes the value $5.10^{20}$ Pa s down to 650 km and again increases by a factor 30 in the lower mantle.

Figure 3 illustrates the first situation, i.e. the computed geoid associated with the dynamically supported topography when the two lithospheric plates are of uniform thickness. One finds complete antisymmetry between converging and diverging zones. The extrema lie on the equator with absolute values reaching 63.6 m. The sensitivity of this amplitude to the relative viscosity distribution, the absolute value of this viscosity or the velocity is similar to what was already stated for the case of the moving round continent.

For comparison, Fig. 4 shows the computed geoid when the model lithosphere thickens as the square root of age. In this second situation the mechanical decoupling due to the asthenospheric wedge beneath the ridge clearly damps the corresponding minimum which drops to 36 m. On the other hand the maximum over the trench is enhanced to 108 m.

This last model is physically meaningful for speculative comparisons with the real Earth. Indeed, it has been noticed by various authors that a marked broad geoid high exists along the circum-Pacific belt and is particularly conspicuous when the degrees 2 and 3 are filtered out (Hager 1984). This geoid high is also present, along the Alpine-Himalayan belt. Although we find that plate velocities can in principle explain these remarkable features, we should note that the very existence of the dense slabs has also been invoked for this purpose (Hager 1984). This last proposal requires a rather stiffer lower-mantle structure. Turning now to the minima of the observed geoid, we notice that they do not clearly overlap with ridges. The reduced mechanical coupling imposed under young oceans can explain this. The aforementioned positive geoid signature linked with the compensated ridge topography (Haxby & Turcotte 1978; Fleitout & Froidevaux 1983) could of course partly hide a broader but weak geoid minimum.

In conclusion we emphasize that the first aim of this paper was to present a mathematical framework to deal with a

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Computed geoid induced by the same plate velocity pattern as in Fig. 3. The difference from the previous model is that the hemispherical rigid plates now have a square root of age dependence for their thickness, which vanishes at the ridge (90°) and reaches 200 km at the converging boundary (270°). Otherwise the viscosity structure is identical. Notice the enhancement of the maximum which reaches 108 m and the damping of the minimum to 36 m.</figcaption></figure>
dynamic Earth with lateral variations of the viscosity structure within its outer shell. Potentially the proposed formalism is capable of incorporating the existence of lithospheric plates with their weak boundaries and their varying thicknesses. The broader general goal for the future would be to build a self-consistent model predicting both surface velocities and geoid undulations, given the distribution of internal loads defined by seismic tomography data. In this first study we ignored the density sources, but imposed surface velocities relying on the fact that for the real Earth the present plate velocity pattern is well known. This type of approach has enabled us, by means of two simplified geometrical configurations, to propose new explanations for several important features of the long wavelength geoid.

## ACKNOWLEDGMENT

This work was supported by the Institut National des Sciences de l'Univers (CNRS) and this publication is the contribution number 1 of its programme 'Dynamique et Bilans de la Terre: couplage global'.

## REFERENCES

* Dziewonski, A. M., 1984. Mapping the lower mantle: determination of lateral heterogeneity in P velocity up to degree and order 6, *J. geophys. Res.*, **89**, 5929–5952.

* Edmonds, A. R., 1960. *Angular Momentum in Quantum Mechanics*, Princeton University Press, New Jersey.

* Forsyth, D. W. & Uyeda, S., 1975. On the relative importance of the driving forces of plate motion, *Geophys. J. R. astr. Soc.*, **43**, 163–200.

* Forsyth, D. W., Parmentier, E. M. & Bowin, 1980. Continental wakes, *EOS*, **61**, 382.

* Fleitout, L. & Froidevaux, C., 1983. Tectonic stresses in the lithosphere, *Tectonics*, **2**, 315–324.

* Forte, A. M. & Peltier, W. R., 1987. Surface plate kinematics and mantle convection, in *The Composition, Structure, and Dynamics of the Lithosphere-Asthenosphere System*, Geodyn. Ser.

* Froidevaux, C. & Ricard, Y., 1987. Tectonic evolution of high plateaus, *Tectonophys.*, **134**, 227–238.

* Froidevaux, C., Schubert, G. & Yuen, D. A., 1977. Thermal and mechanical structure of the upper mantle: a comparison between continental and oceanic models, *Tectonophys.*, **37**, 233–246.

* Hager, B. H., 1984. Subducted slabs and the geoid: constraints on mantle rheology and flow, *J. geophys. Res.*, **89**, 6003–6015.

* Hager, B. H. & O'Connell, R. J., 1978. Kinematic models of large-scale flow in the Earth's mantle, *J. geophys. Res.*, **84**, 1031–1048.

* Hager, B. H., Clayton, R. W., Richards, M. A., Comer, R. P. & Dziewonski, A. M., 1985. Lower mantle heterogeneity, dynamic topography, and the geoid, *Nature*, **313**, 541–545.

* Haxby, W. F. & Turcotte, D. L., 1978. On isostatic geoid anomalies, *J. geophys. Res.*, **83**, 5473–5478.

* Jacoby, W. L. & Seidler, E., 1981. Plate kinematics and the gravity field, *Tectonophys.*, **74**, 155–167.

* Jarvis, G. T. & Peltier, W. R., 1982. Mantle convection as a boundary layer phenomenon, *Geophys. J. R. astr. Soc.*, **68**, 389–427.

* Kaula, W. M., 1975. Product-sum conversion of spherical harmonics with applications to thermal convection, *J. geophys. Res.*, **80**, 225–231.

* Nataf, H. C., Nakanishi, I. & Anderson, D. L., 1987, Measurements of mantle wave velocities and inversion for lateral heterogeneities and anisotropy. 3: inversion, *J. geophys. Res.*, **91**, 7261–7307.

* Phinney, R. A. & Burridge, R., 1973. Representation of the elastic-gravitational excitation of a spherical Earth model by generalized spherical harmonics, *Geophys. J. R. astr. Soc.*, **74**, 451–487.

* Ricard, Y., Fleitout, L. & Froidevaux, C., 1984. Geoid heights and lithospheric stresses for a dynamic earth, *Ann. Geophysicae*, **2**, 267–286.

* Ricard, Y., Hannetelle, B. & Fleitout, L., 1985. Earth seismic tomography, geoid and surface topography, *E.U.G., Terra Cognita*, **5**, 143.

* Richards, M. A. & Hager, B. H., 1986. The Earth's geoid and the large-scale structure of mantle convection, *Proc. NATO Advanced Study Inst., The Physics of Planets*.

* Richardson, R. M., Solomon, C. S. & Sleep, N. H., 1979. Tectonic stresses in the plates, *Rev. Geophys. Space. Phys.*, **17**, 981–1019.

* Sandwell, D. & Schubert, G., 1980. Geoid height versus age for symmetric spreading ridges, *J. geophys. Res.*, **85**, 7235–7240.

## APPENDIX

At any depth within the Earth one can define the components of the flow field $v_r$, $v_\theta$ and $v_\varphi$ and the stress acting on a horizontal plane $\tau_{rr}$, $\tau_{r\theta}$ and $\tau_{r\varphi}$. The horizontal components can be joined to form complex functions of the spherical coordinates $r, \theta, \varphi$

$$V_+ = \frac{1}{\sqrt{2}} (v_\theta + i v_\varphi) \qquad (A1)$$

$$\tau_+ = \frac{i}{\sqrt{2}} (\tau_{r\theta} + i \tau_{r\varphi}). \qquad (A2)$$

With the help of ordinary and generalized spherical harmonics, defined by (1) and (2) in the main text, we have:

$$v_r = \sum_{lm} u_{1lm}(r) Y^{lm}(\theta, \varphi) \qquad (A3)$$

$$V_+ = \sum_{lm} [u_{2lm}(r) + i v_{1lm}(r)] \sqrt{\frac{l(l+1)}{2}} Y^{lm+}(\theta, \varphi) \qquad (A4)$$

$$\tau_{rr} = \sum_{lm} u_{3lm}(r) Y^{lm}(\theta, \varphi) \qquad (A5)$$

$$\tau_+ = \sum_{lm} [u_{4lm}(r) + i v_{2lm}(r)] \sqrt{\frac{l(l+1)}{2}} Y^{lm+}(\theta, \varphi). \qquad (A6)$$

In these expressions we find the poloidal scalar functions called $u_{1lm}(r)$, $u_{2lm}(r)$, $u_{3lm}(r)$ and $u_{4lm}(r)$ and the toroidal scalar functions $v_{1lm}(r)$ and $v_{2lm}(r)$. This notation is in agreement with other authors (Hager & O'Connell 1978; Ricard et al. 1984).

Let us use the symbols ($v_r$), ($V_+$), ($\tau_{rr}$), and ($\tau_+$) to describe vectors made of the $r$-dependent coefficients of $v_r$, $V_+$, $\tau_{rr}$, and $\tau_+$ expressed either by ordinary or generalized spherical harmonics. We can then introduce an operator $F(k_1, k_2)$ which multiplies the real component of degree $l$ of such a vector by $k_1(l)$, and the imaginary part of degree $l$ by $k_2(l)$.

For the inner sphere of our earth model shown in Fig. 1 the above formalism can readily be applied. Indeed the case of a purely radial viscosity structure has been solved on the basis of the Navier-Stokes equation for an incompressible fluid with infinite Prandtl number (Kaula 1975; Hager & O'Connell 1978; Ricard et al. 1984). In the absence of internal loads we deduce that a horizontal velocity pattern with poloidal and toroidal components $u_{2lm}$ and $v_{1lm}$
imposed at the surface of such a sphere creates stresses. The shear stresses at the surface have poloidal and toroidal components $u_{4lm}$ and $v_{2lm}$ associated with the above $u_{2lm}$ and $v_{1lm}$ through proportionality factors $k_H^p(l)$ and $k_H^t(l)$. The vertical stress $\tau_{rr}$ with components $u_{3lm}$ deflects the outer surface of the sphere as well as possible deeper interfaces. This creates a geoid proportional to the poloidal $u_{2lm}$ only, with a proportionality factor $n_H^p(l)$. On the other hand an imposed vertical surface velocity $v_r$ with components $u_{1lm}$ creates a poloidal horizontal shear stress and a geoid, with proportionality factors $k_V^p(l)$ and $n_p V(l)$. These linear relationships can be written in a compact form by means of the operator $F$ defined above:

$$ (\tau_+^L) = F(k_H^p, k_H^t)(V_+^L) + F(k_V^p \sqrt{\frac{l(l+1)}{2}}, 0)(v_r) \quad (A7) $$

$$ (N) = F(n_H^p \sqrt{\frac{2}{l(l+1)}}, 0)(V_+^L) + F(n_V^p, 0)(v_r). \quad (A8) $$

Notice that the square-root term appears only in the proportionality between quantities expressed on different sets of spherical harmonics, i.e. generalized and ordinary. By inspection of these last equations and of equation (3) and (4) in the main text, one can identify the matrices $M_1, M_2, M_3$ and $M_4$.

We now turn our attention to the outer shell with a laterally varying viscosity structure (see Fig. 1). The coupling coefficients $\langle 1/\eta \rangle$ and $c$ found in (6) and (9) can be expanded in spherical harmonics:

$$ \left\langle \frac{1}{\eta} \right\rangle = \sum_{lm} \left\langle \frac{1}{\eta} \right\rangle_{lm} Y^{lm}(\theta, \varphi) \quad (A9) $$

$$ c = \sum_{lm} c_{lm} Y^{lm}(\theta, \varphi). \quad (A10) $$

Using the definitions (A1), (A2) and (A9), we can now write the dynamic relationship (6):

$$ \sum_{l_1 m_1} (V_{l_1 m_1}^0 + -V_{l_1 m_1}^L) Y^{l_1 m_1} + L \sum_{l_2 m_2} \sum_{l_3 m_3} \left\langle \frac{1}{\eta} \right\rangle_{l_2 m_2} \tau_{l_3 m_3}^L Y^{l_2 m_2} Y^{l_3 m_3} + \quad (A11) $$

On the right-hand side, the products of two spherical harmonic functions can be represented by means of $Y^{l_1 m_1 +}$:

$$ Y^{l_2 m_2} Y^{l_3 m_3 +} = \sum_{l_1 m_1} \alpha_{l_1 m_1}^{l_2 m_2 l_3 m_3} Y^{l_1 m_1 +}, \quad (A12) $$

where the coefficients of this expansion are

$$ \alpha_{l_1 m_1}^{l_2 m_2 l_3 m_3} = (-1)^{l_1 + m_1} \sqrt{\frac{(2l_1 + 1)(2l_2 + 1)(2l_3 + 1)}{4\pi}} \times \begin{pmatrix} l_1 & l_2 & l_3 \\ -1 & 0 & 1 \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l_3 \\ -m_1 & m_2 & m_3 \end{pmatrix}. \quad (A13) $$

Here the two brackets stand for the Wigner-3j symbols

(Edmonds 1960). The coefficients $\alpha_{l_1 m_1}^{l_2 m_2 l_3 m_3}$ vanish unless $m_2 = m_1 + m_3$, $|m_2| \le l_2$, and $|l_1 - l_3| \le l_2 \le l_1 + l_3$. This yields

$$ v_{l_1 m_1 +}^0 - v_{l_1 m_1 +}^L = L \sum_{l_3 m_3} \left( \sum_{l_2 m_2} \alpha_{l_1 m_1}^{l_2 m_2 l_3 m_3} \left\langle \frac{1}{\eta} \right\rangle_{l_2 m_2} \right) \tau_{l_3 m_3 +}^L \quad (A14) $$

which, in matrix form, exhibits complete similarity with (6), and corresponds to equation (11) in the main text:

$$ (V_+^0) - (V_+^L) = M_5(\tau_+^L). \quad (A15) $$

The matrix $M_5$ contains the following elements:

$$ m_{l_1 m_1}^{l_3 m_3} = \sum_{l_2 m_2} \alpha_{l_1 m_1}^{l_2 m_2 l_3 m_3} \left\langle \frac{1}{\eta} \right\rangle_{l_2 m_2}. \quad (A16) $$

In addition to this last dynamic equation, we also have equation (7) which deals with the vertical velocity at the base of the outer shell of our earth model. To deal with this last kinematic relationship we first must formulate the vertically-averaged horizontal velocity given by (9) with the help of the same spherical functions. This readily leads to

$$ \langle (V_+) \rangle = C(V_+^0) + (\mathrm{Id} - C)(V_+^L), \quad (A17) $$

where $\mathrm{Id}$ is the identity matrix and the new matrix $C$ is similar to $M_5$ in (A15). It contains the elements

$$ c_{l_1 m_1}^{l_3 m_3} = \sum_{l_2 m_2} \alpha_{l_1 m_1}^{l_2 m_2 l_3 m_3} c_{l_2 m_2}. \quad (A18) $$

To express (7) we have to calculate the horizontal divergence of the averaged horizontal velocity, whose components have just been given. We find that this divergence depends only on the poloidal scalar $u_{2lm}$. Thus,

$$ \nabla_2 \langle \mathbf{V} \rangle = - \sum_{lm} \frac{l(l+1)}{R} \langle u_{2lm} \rangle Y^{lm}(\theta, \varphi). \quad (A19) $$

We therefore need to express $\langle u_{2lm} \rangle$ in terms of $\langle V_+ \rangle$. By inspection of equation (A4) we find that

$$ \langle (u_2) \rangle = F \left( \sqrt{\frac{2}{l(l+1)}}, 0 \right) \langle (V_+) \rangle. \quad (A20) $$

Combining the last two equations we can now write the kinematic relationship (7) in the following form:

$$ (v_r^L) = \frac{L}{R} F(-\sqrt{2l(l+1)}, 0) \langle (V_+) \rangle. \quad (A21) $$

Finally this vertical velocity at the base of the outer shell can be given in terms of the horizontal velocities at the base and surface of the same shell. Substituting (A17) into (A21) we find:

$$ (v_r^L) = \frac{L}{R} F(-\sqrt{2l(l+1)}, 0) \{C(V_+^0) + (\mathrm{Id} - C)(V_+^L)\}. \quad (A22) $$

This is identical to equation (12) and therefore defines the matrices $M_6$ and $M_7$.