Thin elastic shells with variable thickness
for lithospheric flexure of one-plate planets
Mikael Beuthe
Royal Observatory of Belgium, Brussels, Belgium. E-mail: mbeuthe@oma.be
Abstract
Planetary topography can either be modeled as a load supported by the lithosphere, or as a
dynamical effect due to lithospheric flexure caused by mantle convection. In both cases the response
of the lithosphere to external forces can be calculated with the theory of thin elastic plates or shells.
On one-plate planets the spherical geometry of the lithospheric shell plays an important role in the
flexure mechanism. So far the equations governing the deformations and stresses of a spherical shell
have only been derived under the assumption of a shell of constant thickness. However local studies
of gravity and topography data suggest large variations in the thickness of the lithosphere. In this
article we obtain the scalar flexure equations governing the deformations of a thin spherical shell with
variable thickness or variable Young’s modulus. The resulting equations can be solved in succession,
except for a system of two simultaneous equations, the solutions of which are the transverse deflection
and an associated stress function. In order to include bottom loading generated by mantle convection,
we extend the method of stress functions to include loads with a toroidal tangential component. We
further show that toroidal tangential displacement always occurs if the shell thickness varies, even
in the absence of toroidal loads. We finally prove that the degree-one harmonic components of
the transverse deflection and of the toroidal tangential displacement are independent of the elastic
properties of the shell and are associated with translational and rotational freedom. While being
constrained by the static assumption, degree-one loads can deform the shell and generate stresses.
The flexure equations for a shell of variable thickness are useful not only for the prediction of the
gravity signal in local admittance studies, but also for the construction of stress maps in tectonic
analysis.
1 Introduction
Terrestrial planets are flattened spheres only at first sight: close-up views reveal a rich topography with
unique characteristics for each planet. Unity in diversity is found by studying the support mechanism
for topographic deviations from the hydrostatic planetary shape. A simple mechanism, called isostasy,
postulates mountains floating with iceberg-like roots in the higher density mantle. Another simple mech-
anism assumes that mountains stand on a rigid membrane called the mechanical lithosphere. The two
simple models predict very different gravity signals, from very weak in the first to very strong in the
second. The truth lies in-between: an encompassing model views topographic structures as loads on an
elastic shell with finite rigidity, called the elastic lithosphere (the elastic lithosphere is a subset of the
mechanical lithosphere). The rigidity depends on the elastic properties of the rocks and on the apparent
elastic thickness of the lithosphere. The latter parameter is the main objective of many studies, since its
value can be related to the lithospheric composition and temperature. Another important application of
the model of lithospheric flexure is the determination of stress maps, which can then be compared with
the observed distribution of tectonic features.
http://arxiv.org/abs/0704.1627v2
Besides the important assumption of elasticity, the model of lithospheric flexure is often simplified
by two approximations. The first one states that the area to be analyzed is sufficiently small so that the
curved lithospheric shell can be modeled as a flat plate. The second approximation states that the shell
(or plate) is thin, which means that deformations are small and that the elastic thickness is small with
respect to the wavelength of the load.
On Earth, the model of lithospheric flexure has been very successful for the understanding of the
topography of the oceanic plates, whereas the analysis of continental plates is fraught with difficulties
due to the very old and complex structure of the continents. As far as we know, plate tectonics do not
occur at the present time on other terrestrial planets, which are deemed one-plate planets [Solomon ,
1978], although the term single-shell would fit better because of the curvature. A single shell can support
loads of much larger extent and greater weight, the best example of which is the huge Tharsis volcanic
formation covering a large portion of Mars. Indeed such a load cannot be supported by the bending
moments present in a thin flat plate, whereas it can be supported by stresses tangent to the shell: the
shell acts as a membrane.
The first application of thin shell theory to a planet was done by Brotchie and Silvester [1969] for
the lithosphere of the Earth (see also Brotchie [1971]). However the approximation of flat plate theory
was seen to be sufficient when it was understood that the lithosphere of the Earth is broken into several
plates [Tanimoto, 1998]. Brotchie and Silvester [1969] do not consider tangential loads and their flexure
equation only includes dominant terms in derivatives; their equation is thus a special case of the equations
of Kraus [1967] and Vlasov [1964] discussed below. The articles reviewed hereafter use the thin shell
theory of Kraus or Vlasov unless mentioned otherwise.
On the Moon, Solomon and Head [1979] used Brotchie’s equation to study displacement and stress
in mare basins. Turcotte et al. [1981] estimated gravity-topography ratios for the mascons and discussed
the type of stress supporting topographic loads. Arkani-Hamed [1998] modeled the support of mascons
with Brotchie’s equation. Sugano and Heki [2004] and Crosby and McKenzie [2005] estimated the elastic
thickness of the lithosphere from Lunar Prospector data.
On Mars, the dominance of the Tharsis rise in the topography led to numerous applications of
the theory of thin elastic shells to lithospheric flexure. Thurber and Toksöz [1978], Comer et al. [1985],
Hall et al. [1986] and Janle and Jannsen [1986] used Brotchie’s equation to estimate the lithospheric
thickness under Martian volcanoes. Turcotte et al. [1981] studied the transition between bending and
membrane regimes. Willemann and Turcotte [1982] analyzed the lithospheric support of Tharsis. Sleep and Phillips
[1985] analyzed the membrane stress distribution on the whole surface. Banerdt et al. [1992], Banerdt and Golombek
[2000] and Phillips et al. [2001] used a model of lithospheric flexure including membrane support, bending
stresses and tangent loads [Banerdt , 1986] in order to study the global stress distribution. Arkani-Hamed
[2000] determined the elastic thickness beneath large volcanoes with Brotchie’s equation whereas Johnson et al.
[2000] estimated the elastic thickness beneath the North Polar Cap. Using local admittance analysis with
spatiospectral methods, McGovern et al. [2002] determined the elastic thickness at various locations [see
also McGovern et al., 2004]. McKenzie et al. [2002] made local admittance analyses of line-of-sight grav-
ity data both with flat plates and with spherical shell models. Turcotte et al. [2002] used the spherical
shell formula for a one-dimensional wavelet analysis of the admittance in order to determine the average
elastic thickness of the lithosphere. Zhong and Roberts [2003] and Lowry and Zhong [2003] studied the
support of the Tharsis rise with an hybrid model including the flexure of a thin elastic shell as well as the
internal loading of a thermal plume in the mantle. Belleguic et al. [2005] determined the elastic thickness
and the density beneath large volcanoes. Searls et al. [2006] investigated the elastic thickness and the
density beneath the Utopia and Hellas basins.
Venus is considered as a one-plate planet but does not have giant volcanic or tectonic structures
comparable to Tharsis. The spherical shell model has thus not been used as often for Venus as for Mars.
Banerdt [1986] studied the global stress distribution. Janle and Jannsen [1988] and Johnson and Sandwell
[1994] used Brotchie’s equation to estimate the lithospheric thickness in various locations. Sandwell et al.
[1997] computed global strain trajectories for comparison with observed tectonics. Lawrence and Phillips
[2003] inverted the admittance in order to estimate the elastic thickness and mantle density anomalies
over two lowland regions and one volcanic rise. Anderson and Smrekar [2006] established a global map
of the elastic thickness based on local admittance analysis. Mercury’s topography and gravity fields are
not yet known well enough to warrant the application of a thin elastic shell model. We refer to Wieczorek
[2007] for a review of recent results regarding the lithosphere of terrestrial planets.
The elastic thickness of the lithosphere is not at all homogeneous over the surface of a planet. For
example, McGovern et al. [2004] (for Mars) and Anderson and Smrekar [2006] (for Venus) find litho-
spheric thickness variations of more than 100 km. The former study explains the variation in lithospheric
thickness in terms of different epochs of loading, as the lithosphere is thickening with time. Other studies
however inferred that spatial variations in lithospheric thickness on Mars are as important as temporal
variations [Solomon and Head , 1982, 1990; Comer et al., 1985]. Thin spherical shell models have always
been applied with a constant elastic thickness for the whole lithosphere. Local studies are done by win-
dowing both the data (gravity and topography) and the model predictions for gravity [Simons et al.,
1997; Wieczorek and Simons , 2005]. The assumptions behind these methods are that the elastic thick-
ness is constant within the window and that the area outside the window can be neglected. The first
assumption is of course true for a small enough window, but the size of the window is limited from below
by the resolution of the data [Wieczorek and Simons , 2005]. Even if the first assumption were true, the
second assumption is violated in two ways (unless the elastic thickness is spatially constant). First, the
deformation of the shell within the window as well as the associated stress field are both modified if the
elastic thickness is changed in the area outside the window. Second, the value of the predicted gravity
field within the window depends on the shell deflection outside the window.
These reasons make it interesting to develop a model of the lithospheric flexure for a spherical shell
of variable thickness. Although a full inversion of the gravity and topography data is impractical with
such a model because of the huge size of the parameter space, other applications are of high interest. For
example a two-stage inversion can be considered: in the first stage a constant elastic thickness is assumed,
and the resulting values are used in the second stage as a starting point for an inversion with variable
elastic thickness (the parameter space can also be constrained with an a priori). Moreover this model
can be used to produce synthetic data and thus allows us to check the validity of inversions assuming a
constant elastic thickness. Finally, stress and strain fields can be computed for given variations of the
elastic thickness, with the aim of comparing stress and strain maps with tectonic features.
General equations governing the deformations of a thin elastic shell have been given by various
authors [e.g. Love, 1944; Vlasov , 1964; Kraus, 1967]. However the possibility of variable shell thickness
is only considered at the early stage where the strain-displacement relationships, Hooke’s law and the
equilibrium equations are separately derived for the thin shell. The combination of these three sets of
equations into a unique equation for the transverse deflection is made under the restriction of constant
elastic thickness, ‘owing to the analytical complications which would otherwise arise’ [see Kraus, 1967,
p. 199].
This article is dedicated to the derivation of the minimum set of equations governing the deformations
of a thin elastic spherical shell with variable thickness. Using Kraus’ method of stress functions, we find
that the transverse deflection is the solution of a simultaneous system of two differential equations of the
fourth order. Contrary to the case of constant thickness, these equations cannot be combined due to
the presence of products of derivatives of the thickness and derivatives of the deflection or of the stress
function. In order to include bottom loading generated by mantle convection, we extend the method of
stress functions to include toroidal tangential loads, which were not considered by Kraus [1967]. Non-
toroidal loading is for example generated by tangential lithostatic forces whereas toroidal loading could
be due to mantle flow generating drag at the base of the lithosphere (mantle flow also produces non-
toroidal loading). With applications to tectonics in mind, we derive the equations relating the tangential
displacements and the stresses to the transverse deflection and the stress functions. We further show
that toroidal tangential displacement occurs even if there is no toroidal loading (unless the shell thickness
is constant). Finally we prove three properties specific to the degree-one harmonic components: (1)
the degree-one transverse deflection and the degree-one toroidal tangential displacement drop from the
elasticity equations because they represent rigid displacements of the whole shell, (2) the transverse and
tangential components of degree-one loads are related so that the shell does not accelerate, (3) degree-
one loads can deform the shell and generate stresses. Though our aim is to introduce a variable shell
thickness, the final equations are also valid for a variable Young’s modulus (Poisson’s ratio must be kept
constant).
Another way to take into account variations of the lithospheric thickness consists in treating the
lithosphere as a three-dimensional spherical solid that is elastic [Métivier et al., 2006] or viscoelastic [e.g.
Zhong et al., 2003; Latychev et al., 2005; Wang and Wu, 2006]. The resulting equations are exact (no
thin shell approximation) and can be solved with finite element methods. The thin shell assumption
is probably satisfied for known planetary lithospheres; in case of doubt, it is advisable to compare the
results of thin shell theory with three-dimensional models assuming constant elastic thickness: static
thick shell models [e.g. Banerdt et al., 1982; Janes and Melosh , 1990; Arkani-Hamed and Riendler , 2002]
or time-dependent viscoelastic models [e.g. Zhong and Zuber , 2000; Zhong , 2002]. The advantage of thin
shell equations is their two-dimensional character, making them much easier to program and quicker to
solve on a computer. Solving faster either gives access to finer two-dimensional grids or allows to examine
a larger parameter space.
We choose to work within the formalism of differential calculus on curved surfaces, without which
the final equations would be cumbersome. Actually the only tool used in this formalism is the covariant
derivative, which can be seen by geophysicists as just a way of combining several terms into one ‘deriva-
tive’. All necessary formulas are given in the Appendix. The possibility of writing a differential equation
in terms of covariant derivatives (in a tensorial form) also provides a consistency check. The presence of
derivatives that cannot be included into covariant derivatives is simply forbidden. This is not without
meaning for differential equations of the fourth order including products of derivatives.
In section 2, we show how to obtain the strain-displacement relationships, Hooke’s law and the
equilibrium equations for a thin spherical shell. These equations are available in the literature for the
general case of a thin shell [e.g. Kraus, 1967], but we derive them for the spherical case in a simpler way,
starting directly with the metric for the spherical shell. We examine in detail the various approximations
made to obtain the thin shell theory of flexure, refraining until the end from taking the ‘thin shell’
quantitative limit in order to ascertain its influence on the final equations. In section 3, we use the
method of stress functions to obtain the flexure equations governing the displacements and the stresses.
In section 4, we give the final form of the flexure equations in the thin shell approximation. We also study
the covariance and the degree-one projection of the flexure equations. In section 5, we examine various
limit cases in which the flexure equations take a simpler form: the membrane limit, the Euclidean limit
and the limit of constant thickness.
2 Fundamental equations of elasticity
2.1 Three-dimensional elasticity theory
Linear elasticity theory is based on three sets of equations. We directly state them in tensorial form for
an isotropic material, since they are derived in Cartesian coordinates in many books [e.g. Ranalli , 1987;
Synge and Schild , 1978]. Recall that, in tensorial notation, there is an implicit summation on indices
that are repeated on the same side of an equation. The first set of equations includes strain-displacement
relationships:
ǫij =
(ui,j + uj,i) , (1)
where ǫij is the infinitesimal strain tensor and ui are the finite displacements. The ‘comma’ notation
denotes the spatial derivative (see Appendix 7.2).
The second set includes the constitutive equations of elasticity, or Hooke’s law, relating the strain
tensor and the stress tensor σij :
σij = λ ǫ δij + 2Gǫij , (2)
where ǫ = ǫ11 + ǫ22 + ǫ33 and δij = 1 if i= j, otherwise it equals zero. The parameter λ is known as the
first Lamé constant. The parameter G is known as the second Lamé constant, or the shear modulus, or
the modulus of rigidity. Boundary conditions are given by
σij nj = Ti , (3)
with nj being the normal unit vector of the surface element and Ti being the surface force per unit area.
The third set includes equations of motion which reduce to equilibrium equations for stresses if the
problem is static:
σij,j = 0 . (4)
Body forces, such as gravity, are assumed to be absent. Both strain and stress tensors are symmetric:
ǫij = ǫji and σij = σji.
These three-dimensional equations do not yet have the right form for the description of the defor-
mations of a two-dimensional spherical shell. Various methods have been used to generate appropriate
equations. Love [1944] and Timoshenko and Woinowsky-Krieger [1964] derive strain-displacements and
equilibrium equations directly on the surface of the sphere (the latter only for the special cases of no
bending or axisymmetrical loading). Sokolnikoff [1956] derives strain-displacement equations in three-
dimensional curvilinear coordinates using an arbitrary diagonal metric but states without proof the equi-
librium equations in curvilinear coordinates. Kraus [1967] uses Sokolnikoff’s form of strain-displacement
equations and derives equilibrium equations for an arbitrary two-dimensional surface using Hamilton’s
principle (i.e. virtual displacements).
Instead of directly deriving equations on the two-dimensional surface of the sphere, we will first
obtain their form in three-dimensional curvilinear coordinates and then restrict them to the surface of
the sphere. The first step can elegantly be done through the use of tensors [Synge and Schild , 1978].
Equations (1)-(4) are tensorial with respect to orthogonal transformations, but not with respect to other
coordinate transforms (one reason being the presence of usual derivatives). In other words, they are
only valid in Cartesian coordinates. In a three-dimensional Euclidean space, tensorial equations have
a simplified form in Cartesian coordinates because supplementary terms that make them tensorial with
respect to arbitrary coordinate transformations are zero. The missing terms can be reconstructed by using
a set of rules, such as the replacement of usual derivatives by covariant derivatives and the substitution
of tensorial contraction to sum on components. Correspondence rules lead to the following three sets of
equations:
ǫij =
ui|j + uj|i
, (5)
σij = λ ǫ gij + 2Gǫij , (6)
gjk σij|k = 0 , (7)
where ǫ = gkl ǫkl. The notation ui|j denotes the covariant derivative of ui (see Appendix 7.2). The
metric and its inverse are noted gij and g
ij , respectively. Tensorial components cannot be expressed in
a normalized basis (except for Cartesian coordinates) which is more common for physical interpretation
(see Appendix 7.1). Covariant components in equations (5)-(7) are related to components defined in a
normalized basis (written with a hat) by:
gii ûi ,
ǫij =
giigjj ǫ̂ij ,
σij =
giigjj σ̂ij ,
where there is no implicit summation on repeated indices.
In the next section we will introduce additional assumptions in order to restrict the equations to the
two-dimensional surface of a spherical shell.
2.2 Spherical shell
2.2.1 Assumptions of the thin shell theory
Suppose that the two first coordinates are the colatitude θ and longitude ϕ on the surface of the sphere,
whereas the third coordinate ζ is radial. R is the shell radius.
Assumptions of the thin shell theory are [see Kraus, 1967, chap. 2.2]:
1. The shell is thin (say less than one tenth of the radius of the sphere).
2. The deflections of the shell are small.
3. The transverse normal stress is negligible: σζζ = 0.
4. Normals to the reference surface of the shell remain normal to it and undergo no change of length
during deformation: ǫθζ = ǫϕζ = ǫζζ = 0.
The second assumption allows us to use linear equations to describe the deflections. The third and fourth
assumptions are not fully consistent: we refer to Kraus [1967] for more details. We will relax them in
the derivation of the equations for the deflection of a spherical shell. The crucial assumption is σζζ = 0
which is essential for the restriction of Hooke’s law to the two-dimensional shell. As we will see later, σζζ
cannot be zero since it is related to the non-zero transverse load (besides the fact that it is incompatible
with a vanishing transverse strain). What is absolutely necessary is that σζζ ≪ σii for i = (θ, ϕ). In
section 5.3.3, we will show that this condition is satisfied if the wavelength of the load is much larger
than the thickness of the shell.
The reference surface is the middle surface of the shell. With the aim of integrating out the third
coordinate, a coordinate system is chosen so that the radial coordinate ζ is zero on the reference surface.
The metric is given by
ds2 = (R+ ζ)
dθ2 + sin2 θ dϕ2
+ dζ2 . (8)
Christoffel symbols necessary for the computation of the covariant derivatives are given in Appendix 7.3.
2.2.2 Strain-displacement equations
With the metric (8), the strain-displacement equations (5) become
ǫ̂θθ =
(ûθ,θ + ûζ) ,
ǫ̂ϕϕ =
(csc θ ûϕ,ϕ + cot θ ûθ + ûζ) ,
ǫ̂θϕ =
(csc θ ûθ,ϕ − cot θ ûϕ + ûϕ,θ) , (9)
ǫ̂ζζ = ûζ,ζ ,
ǫ̂θζ =
((R + ζ) ûθ,ζ − ûθ + ûζ,θ) ,
ǫ̂ϕζ =
((R + ζ) ûϕ,ζ − ûϕ + csc θ ûζ,ϕ) ,
where all quantities are given in a normalized basis. The fourth assumption of the thin shell theory
implies that the displacements are linearly distributed across the thickness of the shell, with the transverse
displacement being constant. Displacements can thus be expanded to first order in ζ:
(ûθ, ûϕ, ûζ) = (vθ + ζβθ, vϕ + ζβϕ, w) . (10)
The coefficients (vθ , vϕ, w, βθ, βϕ) are independent of ζ: (vθ, vϕ, w) represent the components of the
displacement vector of a point on the reference surface, whereas (βθ, βϕ) represent the rotations of tangents
to the reference surface oriented along the tangent axes. We determine βθ and βϕ by applying the fourth
assumption of the thin shell theory and the expansion (10) to the last two strain-displacement equations
(vθ − w,θ) ,
(vϕ − csc θ w,ϕ) . (11)
The substitution of the expansion (10) into the fourth strain-displacement equation (ǫ̂ζζ = ûζ,ζ) leads to
the condition ǫ̂ζζ = 0, as postulated by the thin shell theory.
After the insertion of the expansion formulas (10) and (11), the first three strain-displacement equa-
tions (9) become
ǫ̂θθ = ǫ
1 + ζ/R
κ0θ ,
ǫ̂ϕϕ = ǫ
1 + ζ/R
κ0ϕ , (12)
2 ǫ̂θϕ = γ
1 + ζ/R
The extensional strains ǫ0θ, ǫ
ϕ and γ
θϕ are defined by
ǫ0θ =
(vθ,θ + w) ,
ǫ0ϕ =
(csc θ vϕ,ϕ + cot θ vθ + w) , (13)
γ0θϕ =
(vϕ,θ − cot θ vϕ + csc θ vθ,ϕ) .
They represent the normal and shearing strains of the reference surface. The flexural strains κ0θ, κ
ϕ and
τ0 are given by
κ0θ = −
O1 w ,
κ0ϕ = −
O2 w , (14)
τ0 = −
O3 w .
They represent the changes in curvature and the torsion of the reference surface during deformation
[Kraus, 1967]. The differential operators O1,2,3 are defined by
+ 1 ,
O2 = csc2 θ
+ cot θ
+ 1 , (15)
O3 = csc θ
− cot θ
The zero upper index in (ǫ0θ, ǫ
θϕ, κ
0) refers to the reference surface and appears in order to
follow Kraus’ notation.
In Love’s theory, one makes the approximation 1
∼= 1R in equations (12). However this does not
simplify the calculations when the shell is a sphere because it has only one radius of curvature (for a
shell with two radii of curvature, this approximation is a great simplification). Our choice to keep the
factor 1
leads to the same results as the theory of Flügge-Lur’e-Byrne, explained in Kraus [1967,
chap. 3.3a] or Novozhilov [1964, p. 53], in which this factor is expanded up to second order. In any
case the choice between the approximation or the expansion of this factor does not affect the equations
(derived in the next section) relating the stress and moment resultants to the strains: they are the same
is approximated to zeroth order, expanded to second order or fully kept.
2.2.3 Hooke’s law
When the metric is diagonal and the basis is normalized, Hooke’s law (6) becomes
σ̂ii = λ
ǫ̂kk + 2G ǫ̂ii (i = j) ,
σ̂ij = 2G ǫ̂ij (i 6= j) .
There is no implicit summation on repeated indices.
The third assumption of the thin shell theory, σ̂ζζ = 0, can be used to eliminate ǫ̂ζζ from Hooke’s
σ̂θθ =
1− ν2
(ǫ̂θθ + ν ǫ̂ϕϕ) ,
σ̂ϕϕ =
1− ν2
(ǫ̂ϕϕ + ν ǫ̂θθ) .
Young’s modulus E and Poisson’s ratio ν are related to Lamé parameters by
(1 + ν)(1 − 2ν)
2(1 + ν)
In principle, the fourth assumption of the thin shell theory, ǫ̂θζ = ǫ̂ϕζ = 0, leads to σ̂θζ = σ̂ϕζ = 0 but
non-vanishing values must be retained for purposes of equilibrium.
Figure 1: Stress resultants and stress couples acting on a small element of the shell. The directions of the
stress resultants (simple arrows) and the rotation sense of the stress couples (double arrows) correspond
to positive components (tensile stress is positive). Loads (qθ, qϕ, q) act on the reference surface.
The substitution of the expansion (12) into the thin shell approximation of Hooke’s law gives
σ̂θθ =
1− ν2
ǫ0θ + νǫ
1 + ζ/R
κ0θ + ν κ
σ̂ϕϕ =
1− ν2
ǫ0ϕ + νǫ
1 + ζ/R
κ0ϕ + ν κ
, (16)
σ̂θϕ =
2(1 + ν)
γ0θϕ +
1 + ζ/R
with ǫ0θ, ǫ
θϕ, κ
ϕ and τ
0 defined by equations (13) and (14). The expressions for σ̂θζ and σ̂ϕζ will
not be needed.
We now integrate the stress distributions across the thickness h of the shell (see Figure 1). The stress
resultants and couples obtained in this way are defined per unit of arc length on the reference surface:
∫ h/2
σ̂ii (1 + ζ/R) dζ (i = θ, ϕ) ,
Nθϕ = Nϕθ =
∫ h/2
σ̂θϕ (1 + ζ/R) dζ ,
∫ h/2
σ̂iζ (1 + ζ/R) dζ (i = θ, ϕ) , (17)
∫ h/2
σ̂ii (1 + ζ/R) ζ dζ (i = θ, ϕ) ,
Mθϕ = Mϕθ =
∫ h/2
σ̂θϕ (1 + ζ/R) ζ dζ .
We evaluate these integrals with the expansion (16). The tangential stress resultants are
Nθ = K
ǫ0θ + ν ǫ
Nϕ = K
ǫ0ϕ + ν ǫ
, (18)
Nθϕ = K
γ0θϕ .
Explicit expressions for the transverse shearing stress resultants Qi are not needed since these quantities
will be determined from the equilibrium equations. The moment resultants are
Mθ = D
κ0θ + ν κ
ǫ0θ + ν ǫ
Mϕ = D
κ0ϕ + ν κ
ǫ0ϕ + ν ǫ
, (19)
Mθϕ = D
The extensional rigidity K and the bending rigidity D are defined by
1− ν2
, (20)
12(1− ν2)
. (21)
Their dimensionless ratio ξ is a large number,
ξ = R2
, (22)
the inverse of which will serve as an expansion parameter for thin shell theory.
2.2.4 Equilibrium equations
With the metric (8), the components θ, ϕ and ζ of the equilibrium equations (7) respectively become
(R+ ζ)
(sin θ σ̂θθ),θ + σ̂θϕ,ϕ − cos θ σ̂ϕϕ + sin θ σ̂ζθ
+ sin θ
(R+ ζ)
= 0 ,
(R + ζ)
(sin θ σ̂θϕ),θ + σ̂ϕϕ,ϕ + cos θ σ̂θϕ + sin θ σ̂ζϕ
+ sin θ
(R+ ζ)
= 0 ,
(R + ζ)
(sin θ σ̂θζ),θ + σ̂ϕζ,ϕ − sin θ (σ̂θθ + σ̂ϕϕ)
+ sin θ
(R+ ζ)
= 0 ,
where the equations have been multiplied by sin θ (R + ζ), (R + ζ) and sin θ (R + ζ)2, respectively. The
stress components are given in a normalized basis.
The integration on ζ of these three equations in the range [−h/2, h/2] yields the equilibrium equations
for the forces:
(sin θ Nθ),θ +Nθϕ,ϕ − cos θ Nϕ + sin θ Qθ +R sin θ qθ = 0 , (23)
(sin θNθϕ),θ +Nϕ,ϕ + cos θ Nθϕ + sin θ Qϕ +R sin θ qϕ = 0 , (24)
(sin θ Qθ),θ +Qϕ,ϕ − sin θ (Nθ +Nϕ)−R sin θ q = 0 , (25)
where qθ and qϕ are the components of the tangential load vector per unit area of the reference surface:
(R+ ζ)
= R2 qi (i = θ, ϕ) .
We choose the convention that tensile stresses are positive (see Figure 1). The transverse load per unit
area of the reference surface is noted q and is taken to be positive toward the center of the sphere:
(R+ ζ)
= −R2 q . (26)
The first two equilibrium equations for the stresses can also be multiplied by ζ before the integration
to yield the equilibrium equations for the moments:
(sin θMθ),θ +Mθϕ,ϕ − cos θMϕ −R sin θ Qθ = 0 , (27)
(sin θMθϕ),θ +Mϕ,ϕ + cos θMθϕ −R sin θ Qϕ = 0 . (28)
We have neglected small terms in
ζ (R+ ζ)
where i = (θ, ϕ). A third equilibrium equation for
the moments exists but has the form of an identity: Mθϕ =Mϕθ.
3 Resolution
3.1 Available methods
At this stage the elastic theory for a thin spherical shell involves 17 equations: six strain-displacement
relationships (13)-14), six stress-strain relations (18)-(19) making Hooke’s law, and five equilibrium equa-
tions (23)-(23) and (27)-(28). There are 17 unknowns: six strain components (ǫ0θ, ǫ
θϕ, κ
three displacements (w, vθ , vϕ), three tangential stress resultants (Nθ, Nϕ, Nθϕ), two transverse shearing
stress resultants (Qθ, Qϕ), and three moment resultants (Mθ,Mϕ,Mθϕ). The three equations (16) are
also needed if the tangential stresses (σ̂θθ, σ̂θϕ, σ̂ϕϕ) are required. The quantities of primary interest to
us are the transverse deflection and the tangential stresses (sometimes tangential strain is preferred, as
in Sandwell et al. [1997] or Banerdt and Golombek [2000]). We thus want to find the minimum set of
equations that must be solved to determine these quantities.
We are aware of two methods of resolution [Novozhilov , 1964, p. 66]. In the first one, we insert the
strain-displacement relationships into Hooke’s law, and substitute in turn Hooke’s law into the equilibrium
equations. This method yields three simultaneous differential equations for the displacements. Once the
displacements are known, it is possible to compute the strains and the stresses. The second method
supplements the equilibrium equations with the equations of compatibility [Novozhilov , 1964, p. 27] that
relate the partial derivatives of the strain components. It is then convenient to introduce the so-called
stress functions [Kraus , 1967, p. 243], without direct physical interpretation, which serve to define the
stress resultants without introducing the tangential displacements. Equations relating the transverse
displacement and the stress functions are then found by applying the third equation of equilibrium and
the third equation of compatibility (it is also possible to use all three equations of compatibility in order
to directly solve for the stress and moment resultants). Once the transverse displacement and the stress
functions are known, stresses can be computed.
If the shell thickness is constant, the deformations of a thin spherical shell can be completely calcu-
lated with both methods. If the shell thickness is variable, the three equations governing displacements,
obtained with the first method, cannot be decoupled and are not easy to solve. Kraus’ method with
stress functions leads to a system of three equations (relating the transverse displacement and the two
stress functions), in which the first equation is decoupled and solved before the other two. This method
thus provides a system of equations much easier to solve and will be chosen in this article.
When solving the equations, one usually assumes from the beginning the large ξ limit, i.e. 1+ ξ ∼= ξ
where ξ is defined by equation (22). We will only take this limit at the end of the resolution. This
procedure will not complicate the computations, since we have to compute anyway many new terms
because of the variable shell thickness.
3.2 Differential operators
We will repeatedly encounter the operators Oi which intervene in the expressions (14) for the flexural
strains (κ0θ, κ
0). Since we are looking for scalar equations, we need to find out how the operators Oi
can be combined in order to yield scalar expressions, i.e. expressions that are invariant with respect to
changes of coordinates on the sphere.
The first thing is to relate the Oi to tensorial operators. Starting from the covariant derivatives on the
sphere ∇i, we construct the following tensorial differential operators of the second degree in derivatives:
Dij = ∇i∇j + gij , (29)
where ∇i denotes the covariant derivative (see Appendix 7.2). These operators give zero when applied
on spherical harmonics of degree one (considered as scalars):
Dij Y1m = 0 (m = −1, 0, 1) . (30)
This property can be explicitly checked on the spherical harmonics (109) with the metric and the formulas
for the double covariant derivatives given in Appendix 7.4.
In two-dimensional spherical coordinates (θ, ϕ), the three operators Oi defined by equations (15) are
related to the operators Dij acting on a scalar function f through
O1 f = Dθθ f ,
O2 f = csc2 θDϕϕ f , (31)
O3 f = csc θDθϕ f .
The operators O1,2,3f actually correspond to normalized Dijf , i.e. Dijf/(
giigjj).
The usual derivatives of the operators Oi satisfy the useful identities (112)-(113) which are proven
in Appendix 7.7. These identities are the consequence of the path dependence of the parallel transport
of vectors on the curved surface of the sphere.
Invariant expressions are built by contracting all indices of the differential operators in their tensorial
form. The indices can be contracted with the inverse metric gij or with the antisymmetric tensor εij (see
Appendix 7.4), which should not be confused with the strain tensor ǫij . In the following, a and b are
scalar functions on the sphere.
At degree 2, the only non-zero contraction of the Dij is related to the Laplacian (104):
∆′a ≡ gij Dij a
= (∆+ 2) a (32)
= (O1 +O2) a
= a,θ,θ + cot θ a,θ + csc
2 θ a,ϕ,ϕ + 2 a .
At degree 4, a scalar expression symmetric in (a, b) is given by
A(a ; b) ≡ [∆′a][∆′b]− [Dij a][Dij b] ,
= [∆ a][∆ b]− [∇i∇j a][∇i∇j b] + [∆ a] b+ a [∆ b] + 2 a b (33)
= [O1 a][O2 b] + [O2 a][O1 b]− 2 [O3 a][O3 b]
= (a,θ,θ + a)
csc2 θ b,ϕ,ϕ + cot θ b,θ + b
csc2 θ a,ϕ,ϕ + cot θ a,θ + a
(b,θ,θ + b)
−2 csc2 θ (a,θ,ϕ − cot θ a,ϕ) (b,θ,ϕ − cot θ b,ϕ) .
where upper indices are raised with the inverse metric: Dij = gikgjlDkl. The action of an operator does
not extend beyond the brackets enclosing it.
If a is constant, A(a ; b) = a∆′b. It is useful to define an associated operator A0 that gives zero if
its first argument is constant:
A0(a ; b) = A(a ; b)− a [∆′b] . (34)
A scalar expression of degree 4 antisymmetric in (a, b) is given by
B1(a ; b) ≡ gij εkl [Dik a] [Djl b]
= gij εkl [∇i∇k a] [∇j∇l b] (35)
= [(O1 −O2) a] [O3b]− [O3a] [(O1 −O2) b]
= csc θ
a,θ,θ − csc2 θ a,ϕ,ϕ − cot θ a,θ
(b,θ,ϕ − cot θ b,ϕ)
− csc θ (a,θ,ϕ − cot θ a,ϕ)
b,θ,θ − csc2 θ b,ϕ,ϕ − cot θ b,θ
We will also need another operator of degree 4:
B2(a ; b) ≡ εij [∇ia ] [∇j ∆′b]
= csc θ
a,θ [∆
′b],ϕ − a,ϕ [∆
′b],θ
. (36)
The sum of the operators B1 and B2 is noted B:
B(a ; b) = B1(a ; b) + B2(a ; b) . (37)
If either a or b is constant, B(a ; b) = 0.
The operators A and B have an interesting property proven in Appendix 7.8: for arbitrary scalar
functions a and b, A(a ; b) and B(a ; b) do not have a degree-one term in their spherical harmonic expansion.
This is not true of A0, B1 and B2.
3.3 Transverse displacement
3.3.1 Resolution of the equations of equilibrium
The first step consists in finding expressions for the moment resultants (Mθ,Mϕ,Mθϕ) in terms of the
transverse displacement and the stress resultants. The extensional strains (ǫ0θ, ǫ
θϕ) can be eliminated
from the equations for stress and moment resultants (18)-(19). The flexural strains (κ0θ, κ
0) depend
on the transverse displacement w through equations (14). We thus obtain
Mθ = −
(O1 + νO2)w +
Mϕ = −
(O2 + νO1)w +
Nϕ , (38)
Mθϕ = −
(1− ν)O3 w +
Nθϕ ,
where the parameter ξ is defined by equation (22).
The second step consists in solving the equilibrium equations for moments in order to find the
transverse shearing stress resultants (Qθ, Qϕ). We substitute expressions (38) into equations (27)-(28).
Knowing that the stress resultants satisfy the equilibrium equations (23)-(24), we obtain new expressions
for Qθ and Qϕ (identities (112)-(113) are helpful):
Qθ = −
(D∆′w),θ − (1− η)Rqθ
(1− ν) (D,θ O2 − csc θD,ϕO3) w −
(η,θNθ + csc θ η,ϕNθϕ) , (39)
Qϕ = −
csc θ (D∆′w),ϕ − (1− η)Rqϕ
(1− ν) (−D,θ O3 + csc θD,ϕ O1) w −
(η,θ Nθϕ + csc θ η,ϕNϕ) . (40)
The operator ∆′ is defined by equation (32) and η is a parameter close to 1:
1 + ξ
. (41)
The third step consists in finding expressions for the stress resultants (Nθ, Nϕ, Nθϕ) in terms of
stress functions by solving the first two equilibrium equations (23)-(24). Let us define the following linear
combinations of the stress and moment resultants:
(Pθ, Pϕ, Pθϕ) =
Mθ, Nϕ +
Mϕ, Nθϕ +
. (42)
We observe that these linear combinations satisfy simplified equations of equilibrium:
(sin θ Pθ),θ + Pθϕ,ϕ − cos θ Pϕ +R sin θ qθ = 0 ,
(sin θ Pθϕ),θ + Pϕ,ϕ + cos θ Pθϕ +R sin θ qϕ = 0 . (43)
Comparing these equilibrium equations with the identities (112)-(113), we see that the homogeneous
equations (i.e. equations (43) with a zero tangential load qθ = qϕ = 0) are always satisfied if
(Pθ, Pϕ, Pθϕ) = (O2,O1,−O3)F , (44)
where F is an auxiliary function called stress function. For the moment, this function is completely
arbitrary apart from being scalar and differentiable.
Particular solutions of the full equations (43) can be found if we express the tangential load qT =
qθθ̂ + qϕϕ̂ in terms of the surface gradient of a scalar potential Ω (consoidal or poloidal component) and
the surface curl of a vector potential V r̂ (toroidal component):
qT = −
∇̄Ω + 1
∇̄ × (V r̂) . (45)
Surface operators are defined in Appendix 7.5, where the terms consoidal/poloidal are also discussed.
The covariant components of qT are (qθ, sin θqϕ) and can be expressed as − 1RΩ,i +
gjkεikV,j , which
gives
qθ = −
Ω,θ +
R sin θ
V,ϕ ,
sin θ qϕ = −
Ω,ϕ −
sin θ
V,θ . (46)
If the tangential load is consoidal (V = 0), a particular solution of equations (43) is given by
(Pθ, Pϕ, Pθϕ) = (1, 1, 0)Ω . (47)
If the tangential load is toroidal (Ω = 0), a particular solution of equations (43) is given by
(Pθ, Pϕ, Pθϕ) = (2O3,−2O3,O2 −O1)H , (48)
where we have introduced a second stress function H which satisfies the constraint
∆′H = −V + V0 , (49)
where V0 is a constant (identities (112)-(113) are useful). This equation allows us to determine the stress
function H if the toroidal source V is known.
The general solution of the equations (43) is given by the sum of the general solution (44) of the
homogeneous equations and the two particular solutions (47)-(48) of the full equations:
(Pθ, Pϕ, Pθϕ) = (O2F +Ω + 2O3H,O1F +Ω− 2O3H,−O3F + (O2 −O1)H) . (50)
The stress resultants (Nθ, Nϕ, Nθϕ) can now be obtained from (Pθ, Pϕ, Pθϕ) by using equations (42) and
(38):
(Nθ, Nϕ, Nθϕ) = η (Pθ, Pϕ, Pθϕ) + η
(O1 + νO2,O2 + νO1, (1− ν)O3)w , (51)
which finally give
Nθ = η
O2 F +
(∆′ − (1− ν)O2)w +Ω+ 2O3H
Nϕ = η
O1 F +
(∆′ − (1− ν)O1)w +Ω− 2O3H
, (52)
Nθϕ = η
−O3 F + (1 − ν)
O3 w + (O2 −O1)H
The fourth step consists in expressing the third equation of equilibrium (25) in terms of the transverse
displacement w and the stress functions F and H . For this purpose, it is handy to express the transverse
shearing stress resultants (Qθ, Qϕ) in terms of (w,F,H). We thus substitute Nϕ, Nθ, Nθϕ, given by
equations (52), into the expressions for Qθ and Qϕ, given by equations (39)-(40):
Qθ = −
(ηD∆′w),θ +
(ηD),θ O2 − csc θ (ηD),ϕ O3
− (η,θ O2 − csc θ η,ϕO3)F + Ω,θ − (ηΩ),θ − csc θ V,ϕ
− 2 (η,θ O3 + csc θ η,ϕ O2)H + csc θ (η,ϕ ∆′H − η(∆′H),ϕ) , (53)
Qϕ = −
csc θ (ηD∆′w),ϕ +
− (ηD),θ O3 + csc θ (ηD),ϕ O1
− (−η,θ O3 + csc θ η,ϕ O1)F + csc θ (Ω,ϕ − (ηΩ),ϕ) + V,θ
+2 (η,θ O1 + csc θ η,ϕ O3)H − (η,θ ∆′H − η(∆′H),θ) . (54)
The various terms present in Qθ and Qϕ can be classified into generic types according to their
differential structure. Each generic type will contribute in a characteristic way to the third equation
of equilibrium. It is worthwhile to compute the generic contribution of each type before using the full
expressions of the stress resultants with their multiple terms. Identities (112)-(113) are helpful in this
calculation. We thus write the third equation of equilibrium (25) as follows:
I (Qθ ;Qϕ ; 0)− (Nθ +Nϕ)−Rq = 0 , (55)
where the differential operator I is defined for arbitrary expressions (X,Y, Z) by
I(X ;Y ;Z) = csc θ (sin θX),θ + csc θ Y,ϕ − cot θ Z . (56)
This operator must be evaluated for the following generic types present in (Qθ, Qϕ):
I (a,θ ; csc θ a,ϕ ; 0) = ∆a ,
I (− csc θ a,ϕ ; a,θ ; 0) = 0 ,
I (a,θ O2b− csc θ a,ϕ O3b ;−a,θ O3b+ csc θ a,ϕO1b ; 0) = A0(a ; b) , (57)
I (a,θ O3b+ csc θ a,ϕ O2b ;−a,θ O1b− csc θ a,ϕO3b ; 0) = B1(a ; b) ,
I (− csc θ (a,ϕ∆′b − a[∆′b],ϕ) ; a,θ ∆′b− a[∆′b],θ ; 0) = 2B2(a ; b) ,
where a and b are scalar functions. The operators A0, B1 and B2 are defined in section 3.2.
We now substitute (Qθ, Qϕ) and (Nϕ, Nθ) into the third equation of equilibrium (55) and use formulas
(57). We thus obtain the first of the differential equations that relate w and the stress functions F and
∆′ (ηD∆′w)− (1 − ν)A(ηD ;w) +R3 A(η ;F ) + 2R3 B(η ;H) = −R4 q +R3 (∆Ω−∆′(ηΩ)) . (58)
The operators ∆′, A and B are defined in section 3.2.
3.3.2 Compatibility relation
A second equation relating the transverse displacement and the stress functions comes from the compat-
ibility relation which is derived by eliminating (vθ, vϕ) from the strain-displacement equations (13):
sin θ γ0θϕ,ϕ
sin2 θ ǫ0ϕ,θ
+ ǫ0θ,ϕ,ϕ − sin θ cos θ ǫ0θ,θ + 2 sin2 θ ǫ0θ −
sin2 θ
∆′w . (59)
The strain components are related to the stress resultants through Hooke’s law (18). The substitution of
equations (18) into the compatibility equation (59) gives
∆′ (α (Nθ +Nϕ))−
∆′w − (1 + ν)J (αNθ ;αNϕ ;αNθϕ) = 0 . (60)
where α is the reciprocal of the extensional rigidity:
α ≡ 1
K(1− ν2)
. (61)
For arbitrary expressions (X,Y, Z), J (X ;Y ;Z) is defined by
J (X ;Y ;Z) = csc2 θ
sin2 θX,θ
+ Y,ϕ,ϕ + 2 (sin θ Z,ϕ),θ
− cot θ Y,θ + 2 Y . (62)
As in the case of the third equation of equilibrium, it is practical to classify the terms present in
(Nθ, Nϕ, Nθϕ) into generic types and evaluate separately their contribution to the equation of compati-
bility. There are three types of terms for which we must evaluate the operator J (identities (112)-(115)
are helpful in this calculation):
J (a ; a ; 0) = ∆′a ,
J (aO2b ; aO1b ;−aO3b) = A(a ; b) , (63)
J (2aO3b ;−2aO3b ; a (O2 −O1)b) = 2B(a ; b) ,
where a and b are scalar functions. The operators ∆′, A and B are defined in section 3.2.
We now evaluate J in equation (60) with expressions (52) and formulas (63):
J (αNθ ;αNϕ ;αNθϕ) = A(ηα ;F ) +
∆′ (ηαD∆′w)− 1
(1− ν)A(ηαD ;w)
+ ∆′ (ηαΩ) + 2B(ηα ;H) . (64)
The term A(ηαD ;w) can be rewritten with the help of the following equality:
1− ν2
A(ηαD ;w) = ∆′w −A(η ;w) , (65)
since (1 − ν2)ηαD/R2 = η/ξ and (η/ξ),i = −η,i.
We finally substitute Nθ + Nϕ, given by equations (52), into the compatibility equation (60) and
use expressions (64)-(65). We thus obtain the second of the differential equations that relate w and the
stress functions F and H :
∆′ (ηα∆′F )− (1 + ν)A(ηα ;F )−
A(η ;w) − 2(1 + ν)B(ηα ;H) = −(1− ν)∆′(ηαΩ) . (66)
3.4 Tangential displacements
Assuming that the flexure equations (58) and (66) for the transverse displacement w and the stress
functions (F,H) have been solved, we now show how to calculate the tangential displacements.
In analogy with the decomposition of the tangential load in equations (45)-(46), the tangential
displacement can be separated into consoidal and toroidal components:
v = ∇̄S + ∇̄ × (T r̂) , (67)
where S and T are the consoidal and toroidal scalars, respectively. The covariant components of v are
(vθ, sin θ vϕ) and can be expressed as S,i + g
jkεikT,j (see Appendix 7.5), which gives
vθ = S,θ + csc θ T,ϕ ,
sin θ vϕ = S,ϕ − sin θ T,θ . (68)
The strain-displacement equations (13) become
ǫ0θ =
((O1 − 1)S +O3 T + w) ,
ǫ0ϕ =
((O2 − 1)S −O3 T + w) , (69)
γ0θϕ =
(2O3S + (O2 −O1)T ) .
The stress resultants (Nθ, Nϕ, Nθϕ) given by equations (18) become
(∆S + (1 + ν)w − (1− ν) ((O2 − 1)S −O3 T )) ,
(∆S + (1 + ν)w − (1− ν) ((O1 − 1)S +O3 T )) , (70)
Nθϕ =
(2O3 S + (O2 −O1)T ) .
The toroidal potential T cancels in the sum Nθ +Nϕ:
Nθ +Nϕ =
(1 + ν) (∆S + 2w) .
The consoidal displacement potential S can thus be related to (w,F,Ω) by using expressions (52) for the
stress resultants:
∆S = Rηα (1− ν) (∆′F + 2Ω) + η
∆′w − 2w , (71)
where α is defined by equation (61).
It is more difficult to extract the toroidal displacement potential T . When the shell thickness is
constant, decoupled equations for the displacements can be found by suitable differentiation and combi-
nation of the three equilibrium equations (23)-(25) for the stress resultants. This method does not work
if the shell thickness is variable because the resulting equations are coupled. The trick consists in relating
the tangential displacements to (w,F,H) by the way of equations similar to the homogeneous part of the
first two equilibrium equations, but with (Nθ, Nϕ, Nθϕ) replaced by
(Nθ, Nϕ, Nθϕ), so that derivatives
do not mix with derivatives of T . We will thus calculate the following expression in two different
ways (from equations (52) and (70)):
= csc2 θ (−X,ϕ + sin θ Y,θ) ,
where X and Y are defined by
X = sin θ I
Nθϕ ;
Y = sin θ I
Nθϕ ;
Nϕ ;−
The operator I is defined by equation (56).
As before, it is easier to begin with the evaluation of the operator Z for generic contributions:
Z(a ; a ; 0) = 0 ,
Z(aO2b ; aO1b ;−aO3b) = −B(a ; b) ,
Z(2aO3b ;−2aO3b ; a(O2 −O1)b) = 2A(a ; b)−∆′ (a∆′b) ,
On the one hand, the evaluation Z for (Nθ, Nϕ, Nθϕ) given by equations (70) gives
−1− ν
∆∆′ T .
On the other hand, the evaluation of Z for (Nθ, Nϕ, Nθϕ) given by equations (52) gives
+ (1− ν)B
+ 2RA
The equality of the two previous formulas yields the sought equation for T :
∆∆′ T = 2R (1 + ν) (B (ηα ;F )− 2A (ηα ;H) + ∆′ (ηα∆′H)) + 2B (η ;w) . (72)
We have used the relation (η/ξ),i = −η,i.
The equation for T shows that toroidal displacement always occurs when the shell thickness varies.
The right-hand side of equation (72) only vanishes when two conditions are met: (1) there is no toroidal
source (so that the stress function H vanishes) and (2) the shell thickness is constant (so that the terms
in B vanish).
3.5 Stresses
Stresses can be computed from (w,F,H) and Ω by substituting equations (14), (18) and (52) into equa-
tions (16):
σ̂θθ =
(O2 F +Ω+ 2O3H) +
R(1− ν2)
(∆′ − (1− ν)O2)w ,
σ̂ϕϕ =
(O1 F +Ω− 2O3H) +
R(1− ν2)
(∆′ − (1− ν)O1)w , (73)
σ̂θϕ =
(−O3 F + (O2 −O1)H) +
R(1 + ν)
O3w .
Stresses at the surface are obtained by setting ζ = h/2.
4 Flexure equations and their properties
4.1 Thin shell approximation
The flexure equations derived in section 3 already include several assumptions of the thin shell theory,
but not yet the first one stating that the shell is thin. Of course, the three other assumptions can be
seen to be consequences of the first one [see Kraus, 1967, chap. 2.2], but we have not imposed in a
quantitative way the thinness condition on the equations. We thus impose the limit of small h/R or,
equivalently, the limit of large ξ (defined by equation (22)) on the flexure equations for (H,F,w, S, T ).
This procedure amounts to expand η ≈ 1 − 1/ξ (neglecting terms in 1/ξ wherever appropriate) and to
neglect the derivatives of η in equations (49), (58), (66), (71) and (72). We thus obtain the final flexure
equations for the displacements (w, S, T ) and for the stress functions (F,H):
∆′H = −V + V0 , (74)
∆′ (D∆′w)− (1 − ν)A(D ;w) +R3∆′F = −R4 q − 2R3Ω+ R
∆Ω , (75)
∆′ (α∆′F )− (1 + ν)A(α ;F )− 1
∆′w = −(1− ν)∆′(αΩ) + 2(1 + ν)B(α ;H) , (76)
∆S = Rα (1 − ν) (∆′F + 2Ω) + 1
∆′w − 2w , (77)
∆∆′ T = 2R (1 + ν) (B (α ;F )− 2A (α ;H) + ∆′ (α∆′H)) . (78)
Recall that the differential operators ∆′, A and B are defined by equations (32), (33) and (37). The
potential pairs (Ω, V ) and (S, T ) are related to the tangential load and displacement by equations (45)
and (67), respectively. In the second equation, the term R3∆Ω/ξ has been kept since it could be large if
Ω has a short wavelength. For the same reason, the term ∆′w/ξ has been kept in the fourth equation.
In the third and fifth equations, the terms depending on H belong to the right-hand sides since they
can be considered as a source once H has been calculated from the first equation. The same can be said
of the terms depending on w and F in the fourth and fifth equations: they are supposed to be known from
the simultaneous resolution of the second and third equations. The difficulty in solving the equations
thus lies with the two flexure equations for w and F which are linear with non-constant coefficients (all
other equations are linear - in their unknowns - with constant coefficients). Once these two core equations
have been solved, all other quantities are easily derived from them.
The bending rigidity D, defined by equation (21), characterizes the bending regime: the shell locally
bends in a similar way as a flat plate undergoing small deflections with negligible stretching. The pa-
rameter α, defined by equation (61), is the reciprocal of the extensional rigidity K and characterizes the
membrane regime of the shell in which bending moments are negligible and the load is mainly supported
by internal stresses tangent to the shell. Since D and K are respectively proportional to the third and
first power of the shell thickness, the membrane regime (in which the D-depending terms are neglected)
is obtained in the limit of an extremely thin shell. This observation and the fact that such a shell, lacking
rigidity, cannot support bending moments justify the use of the term ‘membrane’.
The weight of the various terms in the equations depends on two competing factors: the magnitude
of the coefficient multiplying the derivative and the number of derivatives. On the one hand, a coefficient
containing D will be smaller than a coefficient containing α−1 since αD/R2 ∼ ξ−1 is a small number
(see equation (22)). On the other hand, a large number of derivatives will increase the weight of the term
if the derived function has a small wavelength. The transition between the membrane and the bending
regimes thus depends on the wavelength of the load: if the load has a large wavelength (with respect to
the shell radius), the flexure of the shell will be well described by equations without the terms depending
on D (see section 5.1), whereas the flexure under loads of small wavelength is well described by equations
keeping only the D-depending terms with the largest number of derivatives (see section 5.2).
Stress functions are associated with membrane stretching and give negligible contributions in the
bending regime. Formulas (52) show that F (respectively H) plays the role of potential for the stress
resultants in the membrane regime when the load is transversal (respectively tangential toroidal). There
is no stress function associated with the tangential consoidal component of the load because it is identical
to Ω, the consoidal potential of the load.
The variation of the shell thickness has two effects. First, it couples the spherical harmonic modes
that are solutions to the equations for a shell of constant thickness. Second, the toroidal part remains
intertwined with the transversal and consoidal parts, whereas it decouples if the thickness is constant.
For example, the toroidal load is a source for the transverse deflection through the term B in the third
flexure equation. Furthermore, the stress function F is a source for the toroidal displacement in the fifth
flexure equation.
For numerical computations, the flexure equations (74)-(78) obtained in the thin shell approximation
(to which we can add the equations (73) for the stresses) are adequate. For this purpose, it is not useful
to keep small terms in 1/ξ since the theory rests on assumptions only true for a thin shell. In the rest of
the article, we will continue to work with the equations (49), (58), (66), (71) and (72) for more generality.
4.2 Covariance
Because of their tensorial form, the flexure equations (74)-(78) are covariant; this is also true of the more
general equations (49), (58), (66), (71) and (72). This property means that their form is valid in all
systems of coordinates on the sphere, though the tensor components (the covariant derivatives of scalar
functions, the metric and the antisymmetric tensor) will have a different expression in each system. The
scalar functions will also have a different dependence on the coordinates in each system. The covariance
of the final equations was expected. We indeed started with tensorial equations in section 2.1; their
restriction to the sphere in principle respected the tensoriality with respect to changes of coordinates
on the sphere. However the covariance of the two-dimensional theory was not made explicit until we
obtained the final equations. This property is thus a strong constraint on the form of the solution and
a check of its validity, though only necessary and not sufficient (other covariant terms may have been
ignored).
Another advantage of the covariant form is the facility to express the final equations in different
systems of coordinates (even non-orthogonal ones) with the aim of solving them. For example, the fi-
nite difference method in spherical coordinates (θ, ϕ) suffers from a very irregular grid and from pole
singularities. These problems can be avoided with the ‘cubed sphere’ coordinate system [Ronchi et al.,
1996]. Operators including covariant derivatives (here the Laplacian and the operators A and B) can ex-
pressed in any system of coordinates whose metric is known. Christoffel symbols and tensorial differential
operators can be computed with symbolic mathematical software.
4.3 Degree one
Displacements of degree one require special consideration. All differential operators acting on (w,F,H)
(that is ∆′, A and B) in the flexure equations (58) and (66) can be expressed in terms of the operators
Dij = ∇i∇j+gij (see section 3.2). The Dij have the interesting property that they give zero when acting
on spherical harmonics of degree one (see equation (30)). Therefore the degree-one terms in the spherical
harmonic expansion of w vanish from the flexure equations. The magnitude of the transverse deflection
of degree one neither depends on the load nor on the elastic properties of the spherical shell.
More generally, the homogeneous (q = Ω = V = 0) flexure equations (49), (58) and (66) are satisfied
by w and F being both of degree one. According to equations (71)-(72), the corresponding tangential
displacement is constrained by ∆S = −2w and ∆∆′T = 0, so that S and T are also of degree one,
with S = w. These conditions lead to vanishing strains (see equations (69)) which indicate a rigid
displacement. The total displacement is then given by
u = w r̂+ ∇̄w + ∇̄ × (T r̂) ,
where w and T are of degree one. In Appendix 7.6, we show that the first two terms represent a
rigid translation whereas the last term represents a rigid rotation. As expected, stresses vanish for
such displacements (see equation (73)). This freedom in translating or rotating the solution reflects the
freedom in the choice of the reference frame (in practice the reference frame is centered at the center of
the undeformed shell). The same freedom of translation is also found in the theory of deformations of a
spherical, radially stratified, gravitating solid [e.g. Farrell , 1972; Greff-Lefftz and Legros , 1997; Blewitt ,
2003].
What can we say about degree-one loading? Let us first examine what happens when the flexure
equations are projected on the spherical harmonics of degree one. Since the operator ∆′ annihilates the
degree one in any spherical harmonic expansion, terms of the form ∆′f vanish when they are projected
on the spherical harmonics of degree one. Moreover, the operators A and B also vanish in this projection
since they do not contain a degree-one term (see equations (116)-(117) of Appendix 7.8). Therefore the
degree-one component of the flexure equation (66) is identically zero, whereas the degree-one component
of the flexure equation (58) is
q + ∇̄ · qT
Yi = 0 (i = x, y, z) , (79)
where dω = sin θ dθ dϕ and Yi are the real spherical harmonics of degree one. The integral is taken over
the whole spherical surface. We have used the relation ∆Ω = −R ∇̄ ·qT derived from equations (45) and
(106).
The first term in the integrand of equation (79) is the projection on the Cartesian axes (x̂, ŷ, ẑ) of
the vector field q r̂:
(q Yx, q Yy, q Yz) = (q r̂ · x̂, q r̂ · ŷ, q r̂ · ẑ) ,
where we have used formulas (109) for the spherical harmonics.
The second term in equation (79) can be rewritten with the identity (102) and Gauss’ theorem (107):
∇̄ · qT
Yi = −
dω qT · ∇̄Yi (i = x, y, z) ,
where the Yi are considered as scalars. Since qT is orthogonal to r̂, the integrand is the projection of qT
on the Cartesian axes (x̂, ŷ, ẑ):
qT · ∇̄Yx,qT · ∇̄Yy,qT · ∇̄Yz
= (qT · x̂,qT · ŷ,qT · ẑ) ,
where we have used formulas (110) for the gradients of the spherical harmonics.
Recalling that the transverse load q was defined positive towards the center of the sphere (see equation
(26)), we define a total load vector q = −q r̂+qT . With the above results, we can rewrite the degree-one
projection (79) as
dω (q · x̂ ,q · ŷ ,q · ẑ) = (0, 0, 0) , (80)
which means that the integral (over the whole spherical surface) of the projection on the coordinate axes
of the total load vector vanishes. This result is the consequence of the static assumption in the equations
of motion (7), since a non-zero sum of the external forces would accelerate the sphere.
In practice, degree-one loads on planetary surfaces are essentially due to mass redistribution [Greff-Lefftz and Legros ,
1997] and have a tangential consoidal component (for example the gravitational force is not directed to-
ward the center of figure of the shell). If the shell thickness is variable, a non-zero Ω of degree one will
induce degrees higher than one in w and in S. If the shell thickness is constant, the degree-one load
drops from the flexure equations (58)-(66) and w is not affected. However, the degree-one Ω generates
(assuming a constant shell thickness) a degree-one tangential displacement through equation (71), so that
S 6= w. Whether the shell thickness is variable or not, a degree-one Ω thus generates a total displacement
which is not only a translation but also a tangential deformation [Blewitt , 2003], in which case stresses
do not vanish as shown by equations (73).
5 Limit cases
5.1 Membrane limit
A shell is in a membrane state of stress if bending moments (Mθ,Mϕ,Mθϕ) can be neglected, in analogy
with a membrane which cannot support bending moments. Equations (19) show that this is true if the
bending rigidity vanishes: D = 0. Consistency with equation (22) imposes the limit of infinite ξ (or
η = 1). With these approximations, the first flexure equation (58) for (w,F,H) becomes
∆′F = −Rq − 2Ω . (81)
The second flexure equation (66) for (w,F,H) becomes
∆′w = ∆′ (α∆′F )− (1 + ν)A(α ;F )− 2(1 + ν)B(α ;H) + (1− ν)∆′(αΩ) . (82)
If q is independent of w, F and w can be successively determined with spherical harmonic transforms
from equations (81) and (82) (though the right-hand side of the latter equation must be computed with
another method). If q has a linear dependence in w (such as when the sphere is filled with a fluid), w
can be eliminated between equations (81) and (82), so that F and w can also be computed in succession
(however the equation for F cannot be solved by a spherical harmonic transform). H is supposed to be
known since equation (49) is not modified and can be solved with spherical harmonics.
The stresses are obtained from equations (73) with the additional approximation of neglecting the
term in ζ:
(σ̂θθ, σ̂ϕϕ, σ̂θϕ) =
(O2F +Ω+ 2O3H,O1F +Ω− 2O3H,−O3F + (O2 −O1)H) .
Bending moments play a small role if the load has a large wavelength. In practice, the threshold at
which bending moments become significant can be evaluated from the constant thickness equation for w
(see section 5.3.1). One must compare the magnitudes of the terms in D and 1/α in the left-hand side
of equation (88). Bending moments are negligible (i.e. the term in D) if the spherical harmonic degree ℓ
of the transverse displacement w is such that
ℓ . k
, (83)
where k = (12(1− ν2))1/4 ≈ 1.8 if we take ν = 1/4. This threshold is about 10 for a planet with a radius
of 3400 km and a lithospheric thickness equal to 100 km. Flexure equations in the membrane limit of a
shell with constant thickness have been used by Sleep and Phillips [1985] to study the lithospheric stress
in the Tharsis region of the planet Mars.
5.2 Euclidean limit
The equation for the deflection of a rectangular plate with variable thickness has been derived by
Timoshenko and Woinowsky-Krieger [1964, p. 173]. We will check here that the Euclidean limit of
our equations gives the same answer.
Let us define the coordinates (x, y) = (Rϕ,R θ′), where θ′ = θ/2 − θ is the latitude. We work in a
small latitude band around the equator (so that θ′ ≪ 1) and in the limit of large spherical radius R and
large ξ (η = 1). Under this change of coordinate, each derivative introduces a factor R so that terms with
the largest number of derivatives dominate. In particular, covariant derivatives can be approximated by
usual derivatives. The surface Laplacian (104) can be approximated as follows:
∆ ≈ R2
≡ R2∆e .
We assume that there are no tangential loads (Ω = V = H = 0). The flexure equations (75)-(76) for the
transverse displacement become:
R4∆e (D∆ew)− (1 − ν)R4 Ae(D ;w) +R5∆eF = −R4 q , (84)
R4 ∆e (α∆eF )− (1 + ν)R4 Ae(α ;F )−R ∆ew = 0 , (85)
where the operator Ae is defined by
Ae(a ; b) = (∆e a) (∆e b)− a,i,j b,i,j
Equation (85) gives a relation between the magnitudes of F and w:
O(R3 ∆eF ) ∼ O(w/α) .
In the large R limit, equation (84) thus becomes
∆e (D∆ew) − (1− ν)Ae(D ;w) = −q ,
which is the equation derived by Timoshenko and Woinowsky-Krieger [1964]. This equation has been
used by Stark et al. [2003], Kirby and Swain [2004] and Pérez-Gussinyé et al. [2004] for the local analysis
of the lithosphere of the Earth. Its one-dimensional version, in which Ae vanishes, has been used by
Sandwell [1984] to describe the flexure of the oceanic lithosphere on Earth and by Stewart and Watts
[1997] to model the flexure at mountain ranges.
5.3 Shell with constant thickness
5.3.1 Displacements
If the thickness of the shell is constant, the toroidal part of the tangential displacement decouples. The
terms in B indeed drop from the flexure equations (58) and (66) so that the equations for (w,F ) depend
only on (q,Ω):
ηD∆′∆′w − (1− ν) ηD∆′w + ηR3 ∆′F = −R4 q +R3 ((1 − η)∆Ω− 2ηΩ) , (86)
∆′∆′F − (1 + ν)∆′F − 1
∆′w = −(1− ν)∆′Ω , (87)
where we used the property A(a ; b) = a∆′b valid for constant a. We eliminate F from these equations
and obtain a sixth order equation relating w to (q,Ω):
ηD∆∆′∆′w +
∆′w = −R4 (∆′ − 1− ν) q +R3
1 + ξ
∆′ − 1− ν
∆Ω . (88)
The elimination of ∆′F between equations (77) and (86) gives an equation relating the consoidal
displacement potential S to (w, q,Ω):
∆S = − 1
1 + ν
1 + ξ
∆∆′w − 2w − (1 − ν)R2α q + 1− ν
1 + ξ
Rα∆Ω .
Terms not including a Laplacian can be eliminated with equation (88), so that we obtain an explicit
solution for S in terms of (w, q,Ω):
1 + ξ
1− ν2
(∆ + 1 + ν)∆′w + w +R2α q − Rα
1 + ξ
(∆− ξ(1 + ν))Ω , (89)
where the integration constant has been set to zero.
Equations (49) and (72) give an equation for the toroidal displacement potential:
∆′ T = −2ηRα (1 + ν)V , (90)
where the integration constant has been set to zero. We have assumed that ∆V 6= 0, otherwise we get
∆′T = 0.
The differential equations given in this section can be solved with spherical harmonics so that the co-
efficients of the spherical harmonic expansions of (w, S, T ) can be expressed in terms of the corresponding
coefficients of the loads (q,Ω, V ) (see Kraus [1967], Turcotte et al. [1981], Banerdt [1986]).
5.3.2 Comparison with the literature
We now compare our equations for a shell of constant thickness with those found in the literature. The
formulas of Banerdt [1986] (taken from the work of Vlasov [1964]) are the most general:
∆3 + 4∆2
(∆ + 2)w = −R4 (∆ + 1− ν) q +R3
∆− 1− ν
∆Ω , (91)
(∆ + 2)χ =
Dξ(1− ν)
∆V . (92)
Banerdt’s notation is slightly different: his formulas are obtained with the substitutions ξ → ψ, Ω → RΩ
and V → RV . The normal rotation χ is proportional to the radial component (in a normalized basis) of
the curl of the tangential displacement:
(∇× v)r̂ .
The curl ∇× v is related to our surface curl (103) by
∇× v = ∇̄ × v + csc θ
(sin θv̂ϕ),θ − v̂θ,ϕ
With the formulas (68) and (104), we get ∇× v = −∆T r̂ so that equation (92) becomes
∆′ T = −2Rα(1 + ν)V . (93)
We see that Banerdt’s equations (91) and (93) coincide with our equations (88) and (90) in the limit
of large ξ (η = 1), with one exception: the bending term for w is written D(∆3 + 4∆2)w instead of
D∆∆′∆′w = (∆3 + 4∆2 + 4∆)w. This error has propagated in many articles and is of consequence for
the degree-one harmonic component, since it violates the static assumption and spoils the translation
invariance discussed in section 4.3. The impact on higher degrees is negligible. Because of this mistake,
many authors give a separate treatment to the first harmonic degree. Banerdt also gives formulas for the
tangential displacements in terms of consoidal and toroidal scalars (A,B) corresponding to our scalars
(S, T ): his formula (A10) is equivalent to our equation (89) in the limit of large ξ.
If we ignore temperature effects, Kraus’ first equation for (w,F ) is equivalent to our equation
(86) in the limit of large ξ, whereas his second equation for (w,F ) is equivalent to the combination
eq.(87)+ 1+ν
eq.(86) in the limit of large ξ [see Kraus , 1967, eq. 6.54h and 6.55d]. Note that the
definition of Kraus’ stress function F [Kraus, 1967, p. 243] differs from ours:
FKraus = F − k(1 − ν)
with k = 1. This freedom of redefining F for arbitrary k remains as long as D is constant. The flexure
equation for w, equation (88), is unaffected so that the solution for w is unchanged. In the final step,
Kraus makes a mistake when combining the two equations for (w,F ) and thus obtains a flexure equation
for w with the same error as in equation (91). Kraus does not include toroidal loading. The flexure
equation of Turcotte et al. [1981] is taken from Kraus [1967] without the tangential loading and is the
same as equation (91) with Ω = 0.
The flexure equation of Brotchie and Silvester [1969] is given in our notation by
D∆2w +
w = −R4 q , (94)
where q includes their term γw describing the response of the enclosed liquid. This equation can be
obtained from our equation (88) as follows: keep only the derivatives of the highest order in each term,
set Ω = 0, take the limit of large ξ (η = 1) and integrate. Brotchie and Silvester choose to work in
the approximation of a shallow shell and with axisymmetrical loading, solving their equation in polar
coordinates with Bessel-Kelvin functions. The reduction to fourth order in equation (94), the shallow
shell approximation and the axisymmetrical assumption are not justified nowadays since the full equation
(88) can be quickly solved with computer-generated spherical harmonics.
The contraction due to a transverse load of degree 0, w = −R2α(1 − ν)q/2, is equivalent to the
radial displacement computed by Love in the limit of a thin shell [Love, 1944, p.142]. However ad-
ditional assumptions about the initial state of stress and the internal density changes are necessary
[Willemann and Turcotte, 1982] so that the degree 0 is usually excluded from the analysis.
5.3.3 Breakdown of the third assumption of thin shell theory
With the spherical harmonic solutions of the equations for a shell of constant thickness, it is possible to
check the thin shell assumption stating that the transverse normal stress is negligible with respect to the
tangential normal stress. The magnitude of the former can be estimated by the load q (see definition
(26)) whereas the magnitude of the latter can be approximated with formulas (73) evaluated on the outer
surface:
(σ̂θθ + σ̂ϕϕ) |h
∆′F − Eh
4R2(1− ν)
∆′w .
where we have assumed the absence of tangential loads (Ω = 0) and the limit of large ξ. We can relate
σT to q by using the solution in spherical harmonics of equations (87) and (88). Since the thin shell
assumption is expected to fail for a load of sufficiently small wavelength, we assume that the spherical
harmonic degree ℓ is large. Assuming ℓ≫ 1, we obtain
(∆′F )ℓm ≈
wℓm ,
wℓm ≈ −
1 + ℓ
ξ(1−ν2)
qℓm ,
where the spherical harmonic coefficients are indexed by their degree ℓ and their order m. If the shell is
not in a membrane state of stress (see equation (83)), ℓ2 > 2R/h so that σT can be approximated by
(σT )lm ≈
ξ(1 + ν)
qlm .
The thin shell assumption holds if q < σT , that is if
3(1 + ν)
or λ >
3(1 + ν)
h , (95)
where λ is the load wavelength (λ ≈ 2πR/ℓ). We have
3(1 + ν) ≈ 1.9 and 2π/
3(1 + ν) ≈ 3.2 if we
take ν = 1/4. This condition on λ is consistent with the transition zone between the thin and thick shell
responses analyzed in Janes and Melosh [1990] and Zhong and Zuber [2000], but does not coincide with
the constraint given in Willemann and Turcotte [1982], which is ℓ < 2π
R/h (this last condition looks
more like the threshold (83) for the membrane regime).
Though the stress distribution is affected, the limit (95) on the degree ℓ is not important for the
displacements, since they tend to zero at small wavelengths. Therefore the theory does not break down
at short wavelength if one is interested in the computation of the gravity field associated to the transverse
deflection of the lithosphere.
6 Conclusion
The principal results of this article are the five flexure equations (74)-(78) governing the three displace-
ments of the thin spherical shell and the two auxiliary stress functions. Stresses are derived quantities
which can be obtained from equations (73). The shell thickness and Young’s modulus can vary, but
Poisson’s ratio must be constant. The loads acting on the shell can be of any type since we extend the
method of stress functions to include not only transverse and consoidal tangential loads, but also toroidal
tangential loads. The flexure equations can be solved one after the other, except the two equations
(75)-(76) for the transverse deflection w and the stress function F , which must be simultaneously solved.
Tangential loading is usually neglected when solving for the deflection because of its small effect. In that
case, it is sufficient to solve the two equations (75)-(76) with Ω = H = 0:
∆′ (D∆′w)− (1− ν)A(D ;w) +R3∆′F = −R4 q ,
∆′ (α∆′F )− (1 + ν)A(α ;F )− 1
∆′w = 0 .
However tangential loading must be taken into account when computing stress fields [Banerdt , 1986].
In the long-wavelength limit (i.e. membrane regime), all equations can be solved one after the other
because it is possible to solve for F before solving for the transverse deflection. If a small part of the shell
is considered, the flexure equations reduce to the equations governing the deflection of a flat plate with
variable thickness. If the shell thickness is constant, the flexure equations reduce to equations available
in the literature which can be completely solved with spherical harmonics. Our rigorous treatment
of the thin shell approximation has clarified the effect of the shell thickness on the flexure equations.
We emphasize the need to use the correct form for the equations (without the common mistake in the
differential operator acting on w) in order to have the correct properties for the degree-one deflection and
degree-one load.
We have also obtained two general properties of the flexure equations. First we have shown that
there is always a toroidal component in the tangential displacement if the shell thickness is variable.
Second we have proven that the degree-one harmonic components of the transverse deflection and of the
toroidal component of the tangential displacement do not depend on the elastic properties of the shell.
This property reflects the freedom under translations and rotations of the reference frame. Besides we
have shown that degree-one loads are constrained by the static assumption but can deform the shell and
generate stresses.
This article was dedicated to the theoretical treatment of the flexure of a thin elastic shell with
variable thickness. While the special case of constant thickness admits an analytical solution in terms of
spherical harmonics, the general flexure equations must be solved with numerical methods such as finite
differences, finite elements or pseudospectral methods. In a forthcoming paper, we will give a practical
method of solution and discuss applications to real cases.
Acknowledgments
M. Beuthe is supported by a PRODEX grant of the Belgian Science Federal Policy. The author thanks
Tim Van Hoolst for his help and Jeanne De Jaegher for useful comments. Special thanks are due to
Patrick Wu for his constructive criticisms which helped to improve the manuscript.
7 Appendix
7.1 Covariant, contravariant and normalized components
Tensors can be defined by their transformation law under changes of coordinates. The two types of
tensor components, namely covariant and contravariant components, transform in a reciprocal way under
changes of coordinates. Tensor components cannot be expressed in a normalized basis: the space must
have a coordinate vector basis (for contravariant components) and a dual basis (for covariant components)
which are not normalized.
The only exception is a flat space with Cartesian coordinates, where covariant, contravariant and
normalized components are identical. Since the metric is the scalar product of the elements of the
coordinate vector basis, the covariant components are related to components defined in a normalized
basis (written with a hat) by
gii ûi ,
whereas the relation for contravariant components is
ûi .
The normalized Cartesian basis (x̂, ŷ, ẑ) is related to the normalized basis for spherical coordinates
(r̂, θ̂, ϕ̂) by
x̂ = cos θ cosϕ θ̂ − sinϕ ϕ̂+ sin θ cosϕ r̂ ,
ŷ = cos θ sinϕ θ̂ + cosϕ ϕ̂+ sin θ sinϕ r̂ , (96)
ẑ = − sin θ θ̂ + cos θ r̂ .
7.2 Covariant derivatives
Usual derivatives are indicated by a ‘comma’:
vi,j =
Covariant derivatives (defined below) are indicated by a ‘bar’ or by the operator ∇i:
vi|j = ∇j vi .
The former notation emphasizes the tensorial character of the covariant derivative since the covariant
derivative adds a covariant index to the vector. The latter notation is more adapted when we are interested
by the properties of the operator.
The covariant derivative of a scalar function f is equal to the usual derivative, f|i = f,i, and is itself
a covariant vector: f|i = vi. Covariant derivatives on covariant and contravariant vector components are
defined by
vi|j = vi,j − Γkij vk , (97)
vi |j = v
,j + Γ
k , (98)
where the summation on repeated indices is implicit. The symbols Γkij are the Christoffel symbols of the
second kind [Synge and Schild , 1978]. Their expressions for the metrics used in this article are given in
sections 7.3 and 7.4.
Covariant differentiation of higher order tensors is explained in Synge and Schild [1978] but we only
need the rule for a covariant tensor of second order:
σij|k = σij,k − Γlik σlj − Γljk σil .
If some of the indices of the tensor are contravariant, the rule is changed according to equation (98). The
covariant derivatives of the metric and of the inverse metric are zero: gij|k = 0 and g
7.3 Three-dimensional spherical geometry
The geometry of a thin spherical shell of average radius R can be described with coordinates θ, ϕ and
ζ, respectively representing the colatitude, longitude and radial coordinates. The radial coordinate ζ is
zero on the reference surface (i.e. the sphere of radius R) of the shell. The non-zero components of the
metric are given by
gθθ = (R + ζ)
gϕϕ = (R + ζ)
2 sin2 θ ,
gζζ = 1 .
The non-zero Christoffel symbols are given by
θθ = −(R+ ζ) ,
Γζϕϕ = −(R+ ζ) sin2 θ ,
Γθζθ = Γ
θζ = Γ
ζϕ = Γ
Γθϕϕ = − sin θ cos θ ,
ϕθ = Γ
θϕ = cot θ .
7.4 Two-dimensional spherical geometry
If θ and ϕ respectively represent the colatitude and longitude coordinates, the non-zero components of
the metric on the surface of the sphere are given by
gθθ = 1 ,
gϕϕ = sin
2 θ . (99)
The non-zero Christoffel symbols are given by
Γθϕϕ = − sin θ cos θ ,
ϕθ = Γ
θϕ = cot θ .
The double covariant derivatives of a scalar function f are thus given by
f|θ|θ = f,θ,θ ,
f|θ|ϕ = f|ϕ|θ = f,θ,ϕ − cot θ f,ϕ ,
f|ϕ|ϕ = f,ϕ,ϕ + sin θ cos θ f,θ .
An antisymmetric tensor εij is defined by
εij ≡
det gij ε̄ij ,
where ε̄ij is the antisymmetric symbol invariant under coordinate transformations: ε̄θϕ = −ε̄ϕθ = 1,
ε̄θθ = ε̄ϕϕ = 0 (ε̄ij is usually called a tensor density; Synge and Schild [1978] call it a relative tensor of
weight -1). The non-zero covariant components of εij are given for the metric of the spherical surface by
εθϕ = −εϕθ = sin θ .
The non-zero contravariant components, εij = gikgjlεkl, are given by
εθϕ = −εϕθ = csc θ .
The covariant derivative of the tensor εij is zero: εij|k = 0.
7.5 Gradient, divergence, curl and Laplacian
Various differential operators on the surface of the sphere can be constructed with covariant derivatives.
In this section, f and t are scalar functions defined on the sphere and v is a vector tangent to the sphere.
Backus [1986] gives more details on surface operators and on Helmholtz’s theorem.
As mentioned in Appendix 7.2, the covariant derivative of a scalar function f defined on the sphere
is a covariant vector tangent to the sphere whose components are f,θ and f,ϕ. The surface gradient of f
is the same vector with its components expressed in the normalized basis (θ̂, ϕ̂):
∇̄f = f,θ θ̂ + csc θ f,ϕ ϕ̂ . (100)
The contraction of the covariant derivative with the components of a vector v yields a scalar:
vi |i = v
,θ + cot θ v
θ + vϕ,ϕ .
The surface divergence is the corresponding operation on the vector with its components expressed in the
normalized basis (θ̂, ϕ̂):
∇̄ · v = csc θ
(sin θ v̂θ),θ + v̂ϕ,ϕ
. (101)
Since the result is a scalar, vi
= ∇̄ · v. A useful identity is
∇̄ · (f v) = ∇̄f · v + f ∇̄ · v . (102)
The contraction of the antisymmetric tensor εij with the covariant derivative of a scalar t yields the
covariant components of a vector v:
vi = g
jk εik t,j .
The components are given for the metric (99) by vθ = csc θ t,ϕ and vϕ = − sin θ t,θ. If t is considered
as the radial component of the radial vector t = t r̂ (the covariant radial component is equal to the
normalized one), vi are the non-zero covariant components of the three-dimensional curl of t, which is
tangent to the sphere. This fact justifies the definition of the surface curl of t, which is equal to the
vector v but with components given in the normalized basis (θ̂, ϕ̂):
∇̄ × t = csc θ t,ϕ θ̂ − t,θ ϕ̂ . (103)
The contraction of the double covariant derivative acting on a scalar f defines the surface Laplacian:
∆f = gij f|i|j
= f,θ,θ + cot θ f,θ + csc
2 θ f,ϕ,ϕ . (104)
The surface Laplacian can also be seen as the composition of the surface divergence with the surface
gradient: ∆f = ∇̄ · ∇̄f .
According to Helmholtz’s theorem, a vector tangent to the sphere can be written as the sum of the
surface gradient of a scalar f and the surface curl of a radial vector t r̂:
v = ∇̄f + ∇̄ × (t r̂) . (105)
While t is always called the toroidal scalar (or potential) for v, there is no standard terminology for f .
Backus [1986] calls f the consoidal scalar for v. Some authors [e.g. Banerdt , 1986] call f the poloidal
potential for v. The origin of this use lies in the theory of mantle convection, in which plate tectonics are
assumed to be driven by mantle flow. Under the assumption of an incompressible mantle fluid, the velocity
field of the fluid is solenoidal, i.e. its 3-dimensional divergence vanishes. In such a case, the velocity field
can be decomposed into a poloidal part (∇×∇×(P r̂)) and a toroidal part (∇×(Q r̂)), where differential
operators are 3-dimensional [Backus , 1986]. If the velocity field is tangent to the spherical surface, the
poloidal component at the surface is also the consoidal component [Forte and Peltier , 1987]. However the
fields for which we use Helmholtz’s theorem, i.e. the tangential surface load and the tangential surface
displacement, do not belong to 3-dimensional solenoidal vector fields. We thus prefer to use the term
‘consoidal’.
The surface divergence of v depends only on the consoidal scalar f :
∇̄ · v = ∆f . (106)
The two-dimensional version of Gauss theorem is
dω ∇̄ · v = 0 . (107)
where dω = sin θ dθ dϕ and the integral is taken over the whole spherical surface. It can be proven with
formula (101).
7.6 Rigid displacements
At the surface of a sphere subjected to deformation, the displacement u of a point can be expressed with
the help of Helmholtz’s theorem (105) in terms of three scalar functions (w, S, T ) depending on θ and ϕ:
u = w r̂+ ∇̄S + ∇̄ × (T r̂) . (108)
Strains (and stresses) vanish for rigid displacements. Equations (69) show that strains vanish when
(w, S, T ) are of degree one, with S = w (recall that the operatorsOi annihilate the degree one). Assuming
these conditions, we now show that utransl = w r̂ + ∇̄w represents a rigid translation whereas urot =
∇̄ × (t r̂) represents a rigid rotation of the sphere. We choose as basis the real spherical harmonics of
degree one which form the components of the radial unit vector in Cartesian coordinates:
(Yx, Yy, Yz) = (sin θ cosϕ, sin θ sinϕ, cos θ)
= (x̂, ŷ, ẑ) · r̂ . (109)
We need the surface gradient of the real spherical harmonics which can be computed with formulas (96)
and (100):
∇̄Yx, ∇̄Yy, ∇̄Yz
= (x̂− sin θ cosϕ r̂, ŷ − sin θ sinϕ r̂, ẑ− cos θ r̂) . (110)
If the expansion of w in the degree-one basis is w = a Yx + b Yy + c Yz, then
w r̂+ ∇̄w = a x̂+ b ŷ + c ẑ ,
so that utransl is indeed a rigid translation of the sphere.
If the expansion of T in the degree-one basis is T = a′Yx + b
′Yy + c
′Yz, then
∇̄ × (T r̂) = a′
− sinϕ θ̂ − cos θ cosϕ ϕ̂
cosϕ θ̂ − cos θ sinϕ ϕ̂
+ c′ sin θ ϕ̂ ,
so that urot includes a rigid rotation of the sphere, with (a
′, b′, c′) being the angles of rotation around
the axes (x̂, ŷ, ẑ), respectively. Though urot seems to include a uniform radial expansion, one should
recall that linearized strain-displacement equations are not valid for large displacements. Since strains
vanish, the radial expansion is not physical and urot represents a pure rotation. Finite deformations are
for example discussed in Love [1944][pp. 66-73] and Sokolnikoff [1956][pp. 29-33].
7.7 Differential identities for the operators Oi
The differential operators Oi defined by equations (15) satisfy differential identities useful when obtaining
the flexure equations. They are special cases of differential identities valid in curved spaces. The presence
of curvature makes the parallel transport of vectors path-dependent; this property quantifies the curvature
of space and can be expressed as the lack of commutativity of the covariant derivatives of a vector v:
vi|j|k − vi|k|j = Riljk vl , (111)
where Riljk are the covariant components of the Riemann tensor. On the sphere, the Riemann tensor
has only one independent component that is non-zero, Rθϕθϕ = − sin2 θ. Other components are related
by the symmetries Rαβγδ = −Rβαγδ = −Rαβδγ = Rγδαβ.
The substitution of f,i to vi in the commutation relation (111) provides two differential identities
satisfied by double covariant derivatives acting on scalar functions:
csc θ f|ϕ|ϕ
− csc θ f|ϕ|θ,ϕ − cos θ f|θ|θ + sin θ f,θ = 0 ,
f|θ|θ,ϕ − f|ϕ|θ,θ − cot θ f|ϕ|θ + f,ϕ = 0 .
The replacement in the above equations of the double covariant derivatives by the normalized dif-
ferential operators (31) yields the following identities:
(sin θO2f),θ − (O3f),ϕ − cos θO1f = 0 (I1) , (112)
(sin θO3f),θ − (O1f),ϕ + cos θO3f = 0 (I2) . (113)
These identities can also be directly checked with the definitions (15) of the operators Oi.
The identities (I1)-(I2) can be differentiated to generate identities of higher order. A first useful
identity is obtained from sin θ(I1),θ − (I2),ϕ = 0:
csc2 θ
sin2 θ (O2f),θ
+ (O1f),ϕ,ϕ − 2
sin θ (O3f),ϕ
− cot θ (O1f),θ + 2O1f = ∆
′f . (114)
A second useful identity is obtained from (I1),ϕ + sin θ(I2),θ = 0:
csc2 θ
sin2 θ (O3f),θ
− (O3f),ϕ,ϕ +
sin θ ((O2 −O1) f),ϕ
+ cot θ (O3f),θ − 2O3f = 0 . (115)
7.8 No degree one in operators A and B
We want to prove that the operators A and B defined by equations (33) and (37) do not have any
degree-one term in their spherical harmonic expansion:
dω A(a ; b)Y ∗1p = 0 (p = −1, 0, 1) , (116)
dω B(a ; b)Y ∗1p = 0 (p = −1, 0, 1) , (117)
where (a, b) are arbitrary scalar functions on the sphere, dω = sin θ dθ dϕ and the integral is taken over
the whole spherical surface.
This property is not a straightforward consequence of constructing A and B with Dij as a building
block. Although A and B can be factored into terms without degree one (such as Dija or ∆′a), the
product of the factors may contain degree-one terms in its spherical harmonic expansion.
Without loss of generality, we can prove the above identities with the arguments (a, b) being spherical
harmonics of given degree and order. The general result is then obtained by superposition. Let a and b
be spherical harmonics of order m and n: a ∼ eimϕ and b ∼ einϕ (we will not use their harmonic degree
in the proof). All derivatives with respect to ϕ in the operators A and B can then be replaced with the
rules a,ϕ → ima and b,ϕ → inb. The integral over ϕ in equations (116)-(117) gives
dϕ ei(m+n−p)ϕ = 2π δm+n−p,0 ,
so that the integral is zero unless n = p−m.
First consider the case p = 0 (n = −m), that is the projection on the zonal spherical harmonic of
degree one. We thus have to calculate
dθ A0 and
dθ B0 with
A0 ≡ sin θ cos θA(a ; b) ,
B0 ≡ sin θ cos θ B(a ; b) .
The trick consists in rewriting the integrands as total derivatives:
cos2 θ a,θ b,θ + cot θ
sin2 θ −m2
(ab),θ +
sin2 θ +m2 csc2 θ cos 2θ
B0 = −im
cos θ a,θ b,θ + sin θ a b,θ − csc θ cos2 θ a,θ b+ cos θ a b,θ,θ
The sought integrals are thus given by
dθ A0 =
a,θ b,θ −m2 cot θ (ab),θ +m2 csc2 θ ab
dθ B0 = −im
cos θ a,θ b,θ − csc θ a,θ b+ cos θ a b,θ,θ
where we have dropped the terms containing at least one power of sin θ which vanish at the limits; we
have also replaced cos2 θ and cos 2θ by their value at the limits. The remaining terms can be evaluated
by recalling the dependence in sin θ of the spherical harmonics: a = (sin θ)|m| a0 and b = (sin θ)
|m| b0,
where a0 and b0 are polynomials in cos θ. The only non-zero terms at the limits of the integrals are those
for |m| = 1, in which case we have at the limits: a,θb,θ = a0b0, cot θ (ab),θ = 2a0b0, csc2 θ ab = a0b0,
csc θ a,θb = cos θ a0b0, ab,θ,θ = 0. However these terms cancel in the sums so that the integrals vanish for
all m. This completes the proof for the case p = 0.
Now consider the case p = ±1 (n = −m ± 1), that is the projections on the sectoral spherical
harmonics of degree one. We thus have to calculate
dθ A±1 and
dθ B±1 with
A±1 ≡ sin2 θA(a ; b) ,
B±1 ≡ sin2 θB(a ; b) .
We again write the integrands as total derivatives:
A±1 =
sin θ cos θ a,θ b,θ −
sin θ cos θ +m2
(ab),θ −
cos2 θ ∓ 2m
a,θ b
+ sin2 θ a b,θ + 2m(m∓ 1) cot θ ab
B±1 = −i
(m∓ 1) sin θ a,θ b,θ +m sin θ a b,θ,θ − (m∓ 1) cos θ a,θ b
−m cos θ a b,θ −m(1∓m) csc θ ab
The sought integrals are thus given by
dθ A±1 =
−m2 (ab),θ − (1∓ 2m) a,θ b+ 2m(m∓ 1) cot θ ab
dθ B±1 = i
(m∓ 1) cos θ a,θ b+m cos θ a b,θ +m(1∓m) csc θ ab
where we have dropped the terms containing at least one power of sin θ and replaced cos2 θ by its value
at the limits. The remaining terms can be evaluated as in the case p = 0, but with a = (sin θ)|m| a0 and
b = (sin θ)|m∓1| b0. All terms give zero at the limits of the integrals for all values of m. This completes
the proof for the case p = ±1. We have thus proven the identities (116)-(117).
References
Anderson, F. S., and S. E. Smrekar (2006), Global mapping of crustal and lithospheric thickness on Venus,
J. Geophys. Res., 111, E08006, doi:10.1029/2004JE002395.
Arkani-Hamed, J. (1998), The lunar mascons revisited, J. Geophys. Res., 103, 3709–3739.
Arkani-Hamed, J. (2000), Strength of Martian lithosphere beneath large volcanoes, J. Geophys. Res., 105, 26,713–
26,732, doi:10.1029/2000JE001267.
Arkani-Hamed, J., and L. Riendler (2002), Stress differences in the Martian lithosphere: Constraints on the
thermal state of Mars, J. Geophys. Res., 107, 2–1, doi:10.1029/2002JE001851.
Backus, G. (1986), Poloidal and toroidal fields in geomagnetic field modeling, Rev. Geophys., 24, 75–109.
Banerdt, W. B. (1986), Support of long-wavelength loads on Venus and implications for internal structure, J. Geo-
phys. Res., 91, 403–419.
Banerdt, W. B., and M. P. Golombek (2000), Tectonics of the Tharsis Region of Mars: Insights from MGS
Topography and Gravity, in Lunar and Planetary Institute Conference Abstracts, p. 2038.
Banerdt, W. B., R. S. Saunders, R. J. Phillips, and N. H. Sleep (1982), Thick shell tectonics on one-plate planets
- Applications to Mars, J. Geophys. Res., 87, 9723–9733.
Banerdt, W. B., M. P. Golombek, and K. L. Tanaka (1992), Stress and tectonics on Mars, in Mars, edited by
Kieffer, H. H., Jakosky, B. M.,Snyder, C. W. and Matthews, M. S., pp. 249–297, University of Arizona Press,
Tucson.
Belleguic, V., P. Lognonné, and M. Wieczorek (2005), Constraints on the Martian lithosphere from gravity and
topography data, J. Geophys. Res., 110, E11005, doi:10.1029/2005JE002437.
Blewitt, G. (2003), Self-consistency in reference frames, geocenter definition, and surface loading of the solid
Earth, J. Geophys. Res., 108, 2103, doi:10.1029/2002JB002082.
Brotchie, J. F. (1971), Flexure of a liquid-filled spherical shell in a radial gravity field, Mod. Geol., 3, 15–23.
Brotchie, J. F., and R. Silvester (1969), On crustal flexure, J. Geophys. Res., 74, 5249–5252.
Comer, R. P., S. C. Solomon, and J. W. Head (1985), Mars - Thickness of the lithosphere from the tectonic
response to volcanic loads, Rev. of Geophys., 23, 61–92.
Crosby, A., and D. McKenzie (2005), Measurements of the elastic thickness under ancient lunar terrain, Icarus,
173, 100–107, doi:10.1016/j.icarus.2004.07.017.
Farrell, W. E. (1972), Deformation of the Earth by Surface Loads, Rev. Geophys. Space Phys., 10, 761–797.
Forte, A. M., and W. R. Peltier (1987), Plate tectonics and aspherical earth structure: The importance of
poloidal-toroidal coupling, J. Geophys. Res., 92, 3645–3680.
Greff-Lefftz, M., and L. Legros (1997), Some remarks about the degree-one deformation of the Earth, Geo-
phys. J. Int., 131, 699–723, doi:10.1111/j.1365-246X.1997.tb06607.x.
Hall, J. L., S. C. Solomon, and J. W. Head (1986), Elysium region, Mars - Tests of lithospheric loading models
for the formation of tectonic features, J. Geophys. Res., 91, 11,377–11,392.
Janes, D. M., and H. J. Melosh (1990), Tectonics of planetary loading - A general model and results, J. Geo-
phys. Res., 95, 21,345–21,355.
Janle, P., and D. Jannsen (1986), Isostatic gravity and elastic bending models of Olympus Mons, Mars, Annales
Geophysicae, 4, 537–546.
Janle, P., and D. Jannsen (1988), Tepev Mons on Venus: Morphology and Elastic Bending Models, Earth Moon
Planets, 41, 127–139.
Johnson, C. L., and D. T. Sandwell (1994), Lithospheric flexure on Venus, Geophys. J. Int., 119, 627–647, doi:
10.1111/j.1365-246X.1994.tb00146.x.
Johnson, C. L., S. C. Solomon, J. W. Head, R. J. Phillips, D. E. Smith, and M. T. Zuber (2000), Lithospheric
Loading by the Northern Polar Cap on Mars, Icarus, 144, 313–328, doi:10.1006/icar.1999.6310.
Kirby, J. F., and C. J. Swain (2004), Global and local isostatic coherence from the wavelet transform, Geo-
phys. Res. Lett., 31, L24608, doi:10.1029/2004GL021569.
Kraus, H. (1967), Thin Elastic Shells, John Wiley, New York.
Latychev, K., J. X. Mitrovica, M. E. Tamisiea, J. Tromp, and R. Moucha (2005), Influence of lithospheric thickness
variations on 3-D crustal velocities due to glacial isostatic adjustment, Geophys. Res. Lett., 32, L01304, doi:
10.1029/2004GL021454.
Lawrence, K. P., and R. J. Phillips (2003), Gravity/topography admittance inversion on Venus using niching
genetic algorithms, Geophys. Res. Lett., 30, 1–1, doi:10.1029/2003GL017515.
Love, A. E. H. (1944), A treatise on the mathematical theory of elasticity, 4th edition, Dover, New York.
Lowry, A. R., and S. Zhong (2003), Surface versus internal loading of the Tharsis rise, Mars, J. Geophys. Res.,
108, 5099, doi:10.1029/2003JE002111.
McGovern, P. J., S. C. Solomon, D. E. Smith, M. T. Zuber, M. Simons, M. A. Wieczorek, R. J. Phillips, G. A.
Neumann, O. Aharonson, and J. W. Head (2002), Localized gravity/topography admittance and correlation
spectra on Mars: Implications for regional and global evolution, J. Geophys. Res., 107, 5136, doi:10.1029/
2002JE001854.
McGovern, P. J., S. C. Solomon, D. E. Smith, M. T. Zuber, M. Simons, M. A. Wieczorek, R. J. Phillips, G. A.
Neumann, O. Aharonson, and J. W. Head (2004), Correction to “Localized gravity/topography admittance and
correlation spectra on Mars: Implications for regional and global evolution”, J. Geophys. Res., 109, E07007,
doi:10.1029/2004JE002286.
McKenzie, D., D. N. Barnett, and D.-N. Yuan (2002), The relationship between Martian gravity and topography,
Earth Planet. Sci. Lett., 195, 1–16.
Métivier, L., M. Greff-Lefftz, and M. Diament (2006), Mantle lateral variations and elastogravitational deforma-
tions - I. Numerical modelling, Geophys. J. Int., 167, 1060–1076, doi:10.1111/j.1365-246X.2006.03159.x.
Novozhilov, V. V. (1964), Thin shell theory, 2nd edition, Noordhoff, Groningen.
Pérez-Gussinyé, M., A. R. Lowry, A. B. Watts, and I. Velicogna (2004), On the recovery of effective elastic
thickness using spectral methods: Examples from synthetic data and from the Fennoscandian Shield, J. Geo-
phys. Res., 109, B10409, doi:10.1029/2003JB002788.
Phillips, R. J., M. T. Zuber, S. C. Solomon, M. P. Golombek, B. M. Jakosky, W. B. Banerdt, D. E. Smith,
R. M. E. Williams, B. M. Hynek, O. Aharonson, and S. A. Hauck (2001), Ancient Geodynamics and Global-
Scale Hydrology on Mars, Science, 291, 2587–2591, doi:10.1126/science.1058701.
Ranalli, G. (1987), Rheology of the Earth, Allen and Unwin, Boston.
Ronchi, C., R. Iacono, and P. S. Paolucci (1996), The cubed sphere: a new method for the solution of partial
differential equations in spherical geometry, J. Comput. Phys., 124, 93–114.
Sandwell, D. T. (1984), Thermomechanical evolution of oceanic fracture zones, J. Geophys. Res., 89, 11,401–
11,413.
Sandwell, D. T., C. L. Johnson, F. Bilotti, and J. Suppe (1997), Driving Forces for Limited Tectonics on Venus,
Icarus, 129, 232–244, doi:10.1006/icar.1997.5721.
Searls, M. L., W. B. Banerdt, and R. J. Phillips (2006), Utopia and Hellas basins, Mars: Twins separated at
birth, J. Geophys. Res., 111, E08005, doi:10.1029/2005JE002666.
Simons, M., S. C. Solomon, and B. H. Hager (1997), Localization of gravity and topography: Constraints on the
tectonics and mantle dynamics of Venus, Geophys. J. Int., 131, 24–44, doi:10.1111/j.1365-246X.1997.tb00593.x.
Sleep, N. H., and R. J. Phillips (1985), Gravity and lithospheric stress on the terrestrial planets with reference to
the Tharsis region of Mars, J. Geophys. Res., 90, 4469–4489.
Sokolnikoff, I. S. (1956), Mathematical Theory of Elasticity, McGraw-Hill, New York.
Solomon, S. C. (1978), On volcanism and thermal tectonics on one-plate planets, Geophys. Res. Lett., 5, 461–464.
Solomon, S. C., and J. W. Head (1979), Vertical movement in mare basins - Relation to mare emplacement, basin
tectonics, and lunar thermal history, J. Geophys. Res., 84, 1667–1682.
Solomon, S. C., and J. W. Head (1982), Evolution of the Tharsis province of Mars: The importance of heteroge-
neous lithospheric thickness and volcanic construction, J. Geophys. Res., 87, 9755–9774.
Solomon, S. C., and J. W. Head (1990), Heterogeneities in the thickness of the elastic lithosphere of Mars:
Constraints on heat flow and internal dynamics, J. Geophys. Res., 95, 11,073–11,083.
Stark, C. P., J. Stewart, and C. J. Ebinger (2003), Wavelet transform mapping of effective elastic thickness and
plate loading: Validation using synthetic data and application to the study of southern African tectonics,
J. Geophys. Res., 108, 2558, doi:10.1029/2001JB000609.
Stewart, J., and A. B. Watts (1997), Gravity anomalies and spatial variations of flexural rigidity at mountain
ranges, J. Geophys. Res., 102, 5327–5352, doi:10.1029/96JB03664.
Sugano, T., and K. Heki (2004), Isostasy of the Moon from high-resolution gravity and topography data: Impli-
cation for its thermal history, Geophys. Res. Lett., 31, L24703, doi:10.1029/2004GL022059.
Synge, J. L., and J. N. Schild (1978), Tensor calculus, Dover, New York.
Tanimoto, T. (1998), State of stress within a bending spherical shell and its implications for subducting lithosphere,
Geophys. J. Int., 134, 199–206, doi:10.1046/j.1365-246x.1998.00554.x.
Thurber, C. H., and M. N. Toksöz (1978), Martian lithospheric thickness from elastic flexure theory, Geo-
phys. Res. Lett., 5, 977–980.
Timoshenko, S., and S. Woinowsky-Krieger (1964), Theory of Plates and Shells (2nd edition), McGraw-Hill, New
York.
Turcotte, D. L., R. J. Willemann, W. F. Haxby, and J. Norberry (1981), Role of membrane stresses in the support
of planetary topography, J. Geophys. Res., 86, 3951–3959.
Turcotte, D. L., R. Shcherbakov, B. D. Malamud, and A. B. Kucinskas (2002), Is the Martian crust also the
Martian elastic lithosphere?, J. Geophys. Res., 107, 5091, doi:10.1029/2001JE001594.
Vlasov, V. Z. (1964), General theory of shells and its applications in engineering, NASA Tech. Trans. TT F-99,
1-886.
Wang, H., and P. Wu (2006), Effects of lateral variations in lithospheric thickness and mantle viscosity on glacially
induced surface motion on a spherical, self-gravitating Maxwell Earth, Earth Planet. Sci. Lett., 244, 576–589,
doi:10.1016/j.epsl.2006.02.026.
Wieczorek, M. A. (2007), The gravity and topography of the terrestrial planets, in Treatise on Geophysics (in
press), edited by Schubert, G., Elsevier, Amterdam.
Wieczorek, M. A., and F. J. Simons (2005), Localized spectral analysis on the sphere, Geophys. J. Int., 162,
655–675, doi:10.1111/j.1365-246X.2005.02687.x.
Willemann, R. J., and D. L. Turcotte (1982), The role of lithospheric stress in the support of the Tharsis rise,
J. Geophys. Res., 87, 9793–9801.
Zhong, S. (2002), Effects of lithosphere on the long-wavelength gravity anomalies and their implications for the
formation of the Tharsis rise on Mars, J. Geophys. Res., 107, 5054–+, doi:10.1029/2001JE001589.
Zhong, S., and J. H. Roberts (2003), On the support of the Tharsis Rise on Mars, Earth Planet. Sci. Lett., 214,
1–9, doi:10.1016/S0012-821X(03)00384-4.
Zhong, S., and M. T. Zuber (2000), Long-wavelength topographic relaxation for self-gravitating planets and
implications for the time-dependent compensation of surface topography, J. Geophys. Res., 105, 4153–4164,
doi:10.1029/1999JE001075.
Zhong, S., A. Paulson, and J. Wahr (2003), Three-dimensional finite-element modelling of Earth’s viscoelastic
deformation: effects of lateral variations in lithospheric thickness, Geophys. J. Int., 155, 679–695, doi:10.1046/
j.1365-246X.2003.02084.x.
	Introduction
	Fundamental equations of elasticity
	Three-dimensional elasticity theory
	Spherical shell
	Assumptions of the thin shell theory
	Strain-displacement equations
	Hooke's law
	Equilibrium equations
	Resolution
	Available methods
	Differential operators
	Transverse displacement
	Resolution of the equations of equilibrium
	Compatibility relation
	Tangential displacements
	Stresses
	Flexure equations and their properties
	Thin shell approximation
	Covariance
	Degree one
	Limit cases
	Membrane limit
	Euclidean limit
	Shell with constant thickness
	Displacements
	Comparison with the literature
	Breakdown of the third assumption of thin shell theory
	Conclusion
	Appendix
	Covariant, contravariant and normalized components
	Covariant derivatives
	Three-dimensional spherical geometry
	Two-dimensional spherical geometry
	Gradient, divergence, curl and Laplacian
	Rigid displacements
	Differential identities for the operators Oi
	No degree one in operators A and B

  Planetary topography can either be modeled as a load supported by the
lithosphere, or as a dynamical effect due to lithospheric flexure caused by
mantle convection. In both cases the response of the lithosphere to external
forces can be calculated with the theory of thin elastic plates or shells. On
one-plate planets the spherical geometry of the lithospheric shell plays an
important role in the flexure mechanism. So far the equations governing the
deformations and stresses of a spherical shell have only been derived under the
assumption of a shell of constant thickness. However local studies of gravity
and topography data suggest large variations in the thickness of the
lithosphere. In this article we obtain the scalar flexure equations governing
the deformations of a thin spherical shell with variable thickness or variable
Young's modulus. The resulting equations can be solved in succession, except
for a system of two simultaneous equations, the solutions of which are the
transverse deflection and an associated stress function. In order to include
bottom loading generated by mantle convection, we extend the method of stress
functions to include loads with a toroidal tangential component. We further
show that toroidal tangential displacement always occurs if the shell thickness
varies, even in the absence of toroidal loads. We finally prove that the
degree-one harmonic components of the transverse deflection and of the toroidal
tangential displacement are independent of the elastic properties of the shell
and are associated with translational and rotational freedom. The flexure
equations for a shell of variable thickness are useful not only for the
prediction of the gravity signal in local admittance studies, but also for the
construction of stress maps in tectonic analysis.

Introduction
Terrestrial planets are flattened spheres only at first sight: close-up views reveal a rich topography with
unique characteristics for each planet. Unity in diversity is found by studying the support mechanism
for topographic deviations from the hydrostatic planetary shape. A simple mechanism, called isostasy,
postulates mountains floating with iceberg-like roots in the higher density mantle. Another simple mech-
anism assumes that mountains stand on a rigid membrane called the mechanical lithosphere. The two
simple models predict very different gravity signals, from very weak in the first to very strong in the
second. The truth lies in-between: an encompassing model views topographic structures as loads on an
elastic shell with finite rigidity, called the elastic lithosphere (the elastic lithosphere is a subset of the
mechanical lithosphere). The rigidity depends on the elastic properties of the rocks and on the apparent
elastic thickness of the lithosphere. The latter parameter is the main objective of many studies, since its
value can be related to the lithospheric composition and temperature. Another important application of
the model of lithospheric flexure is the determination of stress maps, which can then be compared with
the observed distribution of tectonic features.
http://arxiv.org/abs/0704.1627v2
Besides the important assumption of elasticity, the model of lithospheric flexure is often simplified
by two approximations. The first one states that the area to be analyzed is sufficiently small so that the
curved lithospheric shell can be modeled as a flat plate. The second approximation states that the shell
(or plate) is thin, which means that deformations are small and that the elastic thickness is small with
respect to the wavelength of the load.
On Earth, the model of lithospheric flexure has been very successful for the understanding of the
topography of the oceanic plates, whereas the analysis of continental plates is fraught with difficulties
due to the very old and complex structure of the continents. As far as we know, plate tectonics do not
occur at the present time on other terrestrial planets, which are deemed one-plate planets [Solomon ,
1978], although the term single-shell would fit better because of the curvature. A single shell can support
loads of much larger extent and greater weight, the best example of which is the huge Tharsis volcanic
formation covering a large portion of Mars. Indeed such a load cannot be supported by the bending
moments present in a thin flat plate, whereas it can be supported by stresses tangent to the shell: the
shell acts as a membrane.
The first application of thin shell theory to a planet was done by Brotchie and Silvester [1969] for
the lithosphere of the Earth (see also Brotchie [1971]). However the approximation of flat plate theory
was seen to be sufficient when it was understood that the lithosphere of the Earth is broken into several
plates [Tanimoto, 1998]. Brotchie and Silvester [1969] do not consider tangential loads and their flexure
equation only includes dominant terms in derivatives; their equation is thus a special case of the equations
of Kraus [1967] and Vlasov [1964] discussed below. The articles reviewed hereafter use the thin shell
theory of Kraus or Vlasov unless mentioned otherwise.
On the Moon, Solomon and Head [1979] used Brotchie’s equation to study displacement and stress
in mare basins. Turcotte et al. [1981] estimated gravity-topography ratios for the mascons and discussed
the type of stress supporting topographic loads. Arkani-Hamed [1998] modeled the support of mascons
with Brotchie’s equation. Sugano and Heki [2004] and Crosby and McKenzie [2005] estimated the elastic
thickness of the lithosphere from Lunar Prospector data.
On Mars, the dominance of the Tharsis rise in the topography led to numerous applications of
the theory of thin elastic shells to lithospheric flexure. Thurber and Toksöz [1978], Comer et al. [1985],
Hall et al. [1986] and Janle and Jannsen [1986] used Brotchie’s equation to estimate the lithospheric
thickness under Martian volcanoes. Turcotte et al. [1981] studied the transition between bending and
membrane regimes. Willemann and Turcotte [1982] analyzed the lithospheric support of Tharsis. Sleep and Phillips
[1985] analyzed the membrane stress distribution on the whole surface. Banerdt et al. [1992], Banerdt and Golombek
[2000] and Phillips et al. [2001] used a model of lithospheric flexure including membrane support, bending
stresses and tangent loads [Banerdt , 1986] in order to study the global stress distribution. Arkani-Hamed
[2000] determined the elastic thickness beneath large volcanoes with Brotchie’s equation whereas Johnson et al.
[2000] estimated the elastic thickness beneath the North Polar Cap. Using local admittance analysis with
spatiospectral methods, McGovern et al. [2002] determined the elastic thickness at various locations [see
also McGovern et al., 2004]. McKenzie et al. [2002] made local admittance analyses of line-of-sight grav-
ity data both with flat plates and with spherical shell models. Turcotte et al. [2002] used the spherical
shell formula for a one-dimensional wavelet analysis of the admittance in order to determine the average
elastic thickness of the lithosphere. Zhong and Roberts [2003] and Lowry and Zhong [2003] studied the
support of the Tharsis rise with an hybrid model including the flexure of a thin elastic shell as well as the
internal loading of a thermal plume in the mantle. Belleguic et al. [2005] determined the elastic thickness
and the density beneath large volcanoes. Searls et al. [2006] investigated the elastic thickness and the
density beneath the Utopia and Hellas basins.
Venus is considered as a one-plate planet but does not have giant volcanic or tectonic structures
comparable to Tharsis. The spherical shell model has thus not been used as often for Venus as for Mars.
Banerdt [1986] studied the global stress distribution. Janle and Jannsen [1988] and Johnson and Sandwell
[1994] used Brotchie’s equation to estimate the lithospheric thickness in various locations. Sandwell et al.
[1997] computed global strain trajectories for comparison with observed tectonics. Lawrence and Phillips
[2003] inverted the admittance in order to estimate the elastic thickness and mantle density anomalies
over two lowland regions and one volcanic rise. Anderson and Smrekar [2006] established a global map
of the elastic thickness based on local admittance analysis. Mercury’s topography and gravity fields are
not yet known well enough to warrant the application of a thin elastic shell model. We refer to Wieczorek
[2007] for a review of recent results regarding the lithosphere of terrestrial planets.
The elastic thickness of the lithosphere is not at all homogeneous over the surface of a planet. For
example, McGovern et al. [2004] (for Mars) and Anderson and Smrekar [2006] (for Venus) find litho-
spheric thickness variations of more than 100 km. The former study explains the variation in lithospheric
thickness in terms of different epochs of loading, as the lithosphere is thickening with time. Other studies
however inferred that spatial variations in lithospheric thickness on Mars are as important as temporal
variations [Solomon and Head , 1982, 1990; Comer et al., 1985]. Thin spherical shell models have always
been applied with a constant elastic thickness for the whole lithosphere. Local studies are done by win-
dowing both the data (gravity and topography) and the model predictions for gravity [Simons et al.,
1997; Wieczorek and Simons , 2005]. The assumptions behind these methods are that the elastic thick-
ness is constant within the window and that the area outside the window can be neglected. The first
assumption is of course true for a small enough window, but the size of the window is limited from below
by the resolution of the data [Wieczorek and Simons , 2005]. Even if the first assumption were true, the
second assumption is violated in two ways (unless the elastic thickness is spatially constant). First, the
deformation of the shell within the window as well as the associated stress field are both modified if the
elastic thickness is changed in the area outside the window. Second, the value of the predicted gravity
field within the window depends on the shell deflection outside the window.
These reasons make it interesting to develop a model of the lithospheric flexure for a spherical shell
of variable thickness. Although a full inversion of the gravity and topography data is impractical with
such a model because of the huge size of the parameter space, other applications are of high interest. For
example a two-stage inversion can be considered: in the first stage a constant elastic thickness is assumed,
and the resulting values are used in the second stage as a starting point for an inversion with variable
elastic thickness (the parameter space can also be constrained with an a priori). Moreover this model
can be used to produce synthetic data and thus allows us to check the validity of inversions assuming a
constant elastic thickness. Finally, stress and strain fields can be computed for given variations of the
elastic thickness, with the aim of comparing stress and strain maps with tectonic features.
General equations governing the deformations of a thin elastic shell have been given by various
authors [e.g. Love, 1944; Vlasov , 1964; Kraus, 1967]. However the possibility of variable shell thickness
is only considered at the early stage where the strain-displacement relationships, Hooke’s law and the
equilibrium equations are separately derived for the thin shell. The combination of these three sets of
equations into a unique equation for the transverse deflection is made under the restriction of constant
elastic thickness, ‘owing to the analytical complications which would otherwise arise’ [see Kraus, 1967,
p. 199].
This article is dedicated to the derivation of the minimum set of equations governing the deformations
of a thin elastic spherical shell with variable thickness. Using Kraus’ method of stress functions, we find
that the transverse deflection is the solution of a simultaneous system of two differential equations of the
fourth order. Contrary to the case of constant thickness, these equations cannot be combined due to
the presence of products of derivatives of the thickness and derivatives of the deflection or of the stress
function. In order to include bottom loading generated by mantle convection, we extend the method of
stress functions to include toroidal tangential loads, which were not considered by Kraus [1967]. Non-
toroidal loading is for example generated by tangential lithostatic forces whereas toroidal loading could
be due to mantle flow generating drag at the base of the lithosphere (mantle flow also produces non-
toroidal loading). With applications to tectonics in mind, we derive the equations relating the tangential
displacements and the stresses to the transverse deflection and the stress functions. We further show
that toroidal tangential displacement occurs even if there is no toroidal loading (unless the shell thickness
is constant). Finally we prove three properties specific to the degree-one harmonic components: (1)
the degree-one transverse deflection and the degree-one toroidal tangential displacement drop from the
elasticity equations because they represent rigid displacements of the whole shell, (2) the transverse and
tangential components of degree-one loads are related so that the shell does not accelerate, (3) degree-
one loads can deform the shell and generate stresses. Though our aim is to introduce a variable shell
thickness, the final equations are also valid for a variable Young’s modulus (Poisson’s ratio must be kept
constant).
Another way to take into account variations of the lithospheric thickness consists in treating the
lithosphere as a three-dimensional spherical solid that is elastic [Métivier et al., 2006] or viscoelastic [e.g.
Zhong et al., 2003; Latychev et al., 2005; Wang and Wu, 2006]. The resulting equations are exact (no
thin shell approximation) and can be solved with finite element methods. The thin shell assumption
is probably satisfied for known planetary lithospheres; in case of doubt, it is advisable to compare the
results of thin shell theory with three-dimensional models assuming constant elastic thickness: static
thick shell models [e.g. Banerdt et al., 1982; Janes and Melosh , 1990; Arkani-Hamed and Riendler , 2002]
or time-dependent viscoelastic models [e.g. Zhong and Zuber , 2000; Zhong , 2002]. The advantage of thin
shell equations is their two-dimensional character, making them much easier to program and quicker to
solve on a computer. Solving faster either gives access to finer two-dimensional grids or allows to examine
a larger parameter space.
We choose to work within the formalism of differential calculus on curved surfaces, without which
the final equations would be cumbersome. Actually the only tool used in this formalism is the covariant
derivative, which can be seen by geophysicists as just a way of combining several terms into one ‘deriva-
tive’. All necessary formulas are given in the Appendix. The possibility of writing a differential equation
in terms of covariant derivatives (in a tensorial form) also provides a consistency check. The presence of
derivatives that cannot be included into covariant derivatives is simply forbidden. This is not without
meaning for differential equations of the fourth order including products of derivatives.
In section 2, we show how to obtain the strain-displacement relationships, Hooke’s law and the
equilibrium equations for a thin spherical shell. These equations are available in the literature for the
general case of a thin shell [e.g. Kraus, 1967], but we derive them for the spherical case in a simpler way,
starting directly with the metric for the spherical shell. We examine in detail the various approximations
made to obtain the thin shell theory of flexure, refraining until the end from taking the ‘thin shell’
quantitative limit in order to ascertain its influence on the final equations. In section 3, we use the
method of stress functions to obtain the flexure equations governing the displacements and the stresses.
In section 4, we give the final form of the flexure equations in the thin shell approximation. We also study
the covariance and the degree-one projection of the flexure equations. In section 5, we examine various
limit cases in which the flexure equations take a simpler form: the membrane limit, the Euclidean limit
and the limit of constant thickness.
2 Fundamental equations of elasticity
2.1 Three-dimensional elasticity theory
Linear elasticity theory is based on three sets of equations. We directly state them in tensorial form for
an isotropic material, since they are derived in Cartesian coordinates in many books [e.g. Ranalli , 1987;
Synge and Schild , 1978]. Recall that, in tensorial notation, there is an implicit summation on indices
that are repeated on the same side of an equation. The first set of equations includes strain-displacement
relationships:
ǫij =
(ui,j + uj,i) , (1)
where ǫij is the infinitesimal strain tensor and ui are the finite displacements. The ‘comma’ notation
denotes the spatial derivative (see Appendix 7.2).
The second set includes the constitutive equations of elasticity, or Hooke’s law, relating the strain
tensor and the stress tensor σij :
σij = λ ǫ δij + 2Gǫij , (2)
where ǫ = ǫ11 + ǫ22 + ǫ33 and δij = 1 if i= j, otherwise it equals zero. The parameter λ is known as the
first Lamé constant. The parameter G is known as the second Lamé constant, or the shear modulus, or
the modulus of rigidity. Boundary conditions are given by
σij nj = Ti , (3)
with nj being the normal unit vector of the surface element and Ti being the surface force per unit area.
The third set includes equations of motion which reduce to equilibrium equations for stresses if the
problem is static:
σij,j = 0 . (4)
Body forces, such as gravity, are assumed to be absent. Both strain and stress tensors are symmetric:
ǫij = ǫji and σij = σji.
These three-dimensional equations do not yet have the right form for the description of the defor-
mations of a two-dimensional spherical shell. Various methods have been used to generate appropriate
equations. Love [1944] and Timoshenko and Woinowsky-Krieger [1964] derive strain-displacements and
equilibrium equations directly on the surface of the sphere (the latter only for the special cases of no
bending or axisymmetrical loading). Sokolnikoff [1956] derives strain-displacement equations in three-
dimensional curvilinear coordinates using an arbitrary diagonal metric but states without proof the equi-
librium equations in curvilinear coordinates. Kraus [1967] uses Sokolnikoff’s form of strain-displacement
equations and derives equilibrium equations for an arbitrary two-dimensional surface using Hamilton’s
principle (i.e. virtual displacements).
Instead of directly deriving equations on the two-dimensional surface of the sphere, we will first
obtain their form in three-dimensional curvilinear coordinates and then restrict them to the surface of
the sphere. The first step can elegantly be done through the use of tensors [Synge and Schild , 1978].
Equations (1)-(4) are tensorial with respect to orthogonal transformations, but not with respect to other
coordinate transforms (one reason being the presence of usual derivatives). In other words, they are
only valid in Cartesian coordinates. In a three-dimensional Euclidean space, tensorial equations have
a simplified form in Cartesian coordinates because supplementary terms that make them tensorial with
respect to arbitrary coordinate transformations are zero. The missing terms can be reconstructed by using
a set of rules, such as the replacement of usual derivatives by covariant derivatives and the substitution
of tensorial contraction to sum on components. Correspondence rules lead to the following three sets of
equations:
ǫij =
ui|j + uj|i
, (5)
σij = λ ǫ gij + 2Gǫij , (6)
gjk σij|k = 0 , (7)
where ǫ = gkl ǫkl. The notation ui|j denotes the covariant derivative of ui (see Appendix 7.2). The
metric and its inverse are noted gij and g
ij , respectively. Tensorial components cannot be expressed in
a normalized basis (except for Cartesian coordinates) which is more common for physical interpretation
(see Appendix 7.1). Covariant components in equations (5)-(7) are related to components defined in a
normalized basis (written with a hat) by:
gii ûi ,
ǫij =
giigjj ǫ̂ij ,
σij =
giigjj σ̂ij ,
where there is no implicit summation on repeated indices.
In the next section we will introduce additional assumptions in order to restrict the equations to the
two-dimensional surface of a spherical shell.
2.2 Spherical shell
2.2.1 Assumptions of the thin shell theory
Suppose that the two first coordinates are the colatitude θ and longitude ϕ on the surface of the sphere,
whereas the third coordinate ζ is radial. R is the shell radius.
Assumptions of the thin shell theory are [see Kraus, 1967, chap. 2.2]:
1. The shell is thin (say less than one tenth of the radius of the sphere).
2. The deflections of the shell are small.
3. The transverse normal stress is negligible: σζζ = 0.
4. Normals to the reference surface of the shell remain normal to it and undergo no change of length
during deformation: ǫθζ = ǫϕζ = ǫζζ = 0.
The second assumption allows us to use linear equations to describe the deflections. The third and fourth
assumptions are not fully consistent: we refer to Kraus [1967] for more details. We will relax them in
the derivation of the equations for the deflection of a spherical shell. The crucial assumption is σζζ = 0
which is essential for the restriction of Hooke’s law to the two-dimensional shell. As we will see later, σζζ
cannot be zero since it is related to the non-zero transverse load (besides the fact that it is incompatible
with a vanishing transverse strain). What is absolutely necessary is that σζζ ≪ σii for i = (θ, ϕ). In
section 5.3.3, we will show that this condition is satisfied if the wavelength of the load is much larger
than the thickness of the shell.
The reference surface is the middle surface of the shell. With the aim of integrating out the third
coordinate, a coordinate system is chosen so that the radial coordinate ζ is zero on the reference surface.
The metric is given by
ds2 = (R+ ζ)
dθ2 + sin2 θ dϕ2
+ dζ2 . (8)
Christoffel symbols necessary for the computation of the covariant derivatives are given in Appendix 7.3.
2.2.2 Strain-displacement equations
With the metric (8), the strain-displacement equations (5) become
ǫ̂θθ =
(ûθ,θ + ûζ) ,
ǫ̂ϕϕ =
(csc θ ûϕ,ϕ + cot θ ûθ + ûζ) ,
ǫ̂θϕ =
(csc θ ûθ,ϕ − cot θ ûϕ + ûϕ,θ) , (9)
ǫ̂ζζ = ûζ,ζ ,
ǫ̂θζ =
((R + ζ) ûθ,ζ − ûθ + ûζ,θ) ,
ǫ̂ϕζ =
((R + ζ) ûϕ,ζ − ûϕ + csc θ ûζ,ϕ) ,
where all quantities are given in a normalized basis. The fourth assumption of the thin shell theory
implies that the displacements are linearly distributed across the thickness of the shell, with the transverse
displacement being constant. Displacements can thus be expanded to first order in ζ:
(ûθ, ûϕ, ûζ) = (vθ + ζβθ, vϕ + ζβϕ, w) . (10)
The coefficients (vθ , vϕ, w, βθ, βϕ) are independent of ζ: (vθ, vϕ, w) represent the components of the
displacement vector of a point on the reference surface, whereas (βθ, βϕ) represent the rotations of tangents
to the reference surface oriented along the tangent axes. We determine βθ and βϕ by applying the fourth
assumption of the thin shell theory and the expansion (10) to the last two strain-displacement equations
(vθ − w,θ) ,
(vϕ − csc θ w,ϕ) . (11)
The substitution of the expansion (10) into the fourth strain-displacement equation (ǫ̂ζζ = ûζ,ζ) leads to
the condition ǫ̂ζζ = 0, as postulated by the thin shell theory.
After the insertion of the expansion formulas (10) and (11), the first three strain-displacement equa-
tions (9) become
ǫ̂θθ = ǫ
1 + ζ/R
κ0θ ,
ǫ̂ϕϕ = ǫ
1 + ζ/R
κ0ϕ , (12)
2 ǫ̂θϕ = γ
1 + ζ/R
The extensional strains ǫ0θ, ǫ
ϕ and γ
θϕ are defined by
ǫ0θ =
(vθ,θ + w) ,
ǫ0ϕ =
(csc θ vϕ,ϕ + cot θ vθ + w) , (13)
γ0θϕ =
(vϕ,θ − cot θ vϕ + csc θ vθ,ϕ) .
They represent the normal and shearing strains of the reference surface. The flexural strains κ0θ, κ
ϕ and
τ0 are given by
κ0θ = −
O1 w ,
κ0ϕ = −
O2 w , (14)
τ0 = −
O3 w .
They represent the changes in curvature and the torsion of the reference surface during deformation
[Kraus, 1967]. The differential operators O1,2,3 are defined by
+ 1 ,
O2 = csc2 θ
+ cot θ
+ 1 , (15)
O3 = csc θ
− cot θ
The zero upper index in (ǫ0θ, ǫ
θϕ, κ
0) refers to the reference surface and appears in order to
follow Kraus’ notation.
In Love’s theory, one makes the approximation 1
∼= 1R in equations (12). However this does not
simplify the calculations when the shell is a sphere because it has only one radius of curvature (for a
shell with two radii of curvature, this approximation is a great simplification). Our choice to keep the
factor 1
leads to the same results as the theory of Flügge-Lur’e-Byrne, explained in Kraus [1967,
chap. 3.3a] or Novozhilov [1964, p. 53], in which this factor is expanded up to second order. In any
case the choice between the approximation or the expansion of this factor does not affect the equations
(derived in the next section) relating the stress and moment resultants to the strains: they are the same
is approximated to zeroth order, expanded to second order or fully kept.
2.2.3 Hooke’s law
When the metric is diagonal and the basis is normalized, Hooke’s law (6) becomes
σ̂ii = λ
ǫ̂kk + 2G ǫ̂ii (i = j) ,
σ̂ij = 2G ǫ̂ij (i 6= j) .
There is no implicit summation on repeated indices.
The third assumption of the thin shell theory, σ̂ζζ = 0, can be used to eliminate ǫ̂ζζ from Hooke’s
σ̂θθ =
1− ν2
(ǫ̂θθ + ν ǫ̂ϕϕ) ,
σ̂ϕϕ =
1− ν2
(ǫ̂ϕϕ + ν ǫ̂θθ) .
Young’s modulus E and Poisson’s ratio ν are related to Lamé parameters by
(1 + ν)(1 − 2ν)
2(1 + ν)
In principle, the fourth assumption of the thin shell theory, ǫ̂θζ = ǫ̂ϕζ = 0, leads to σ̂θζ = σ̂ϕζ = 0 but
non-vanishing values must be retained for purposes of equilibrium.
Figure 1: Stress resultants and stress couples acting on a small element of the shell. The directions of the
stress resultants (simple arrows) and the rotation sense of the stress couples (double arrows) correspond
to positive components (tensile stress is positive). Loads (qθ, qϕ, q) act on the reference surface.
The substitution of the expansion (12) into the thin shell approximation of Hooke’s law gives
σ̂θθ =
1− ν2
ǫ0θ + νǫ
1 + ζ/R
κ0θ + ν κ
σ̂ϕϕ =
1− ν2
ǫ0ϕ + νǫ
1 + ζ/R
κ0ϕ + ν κ
, (16)
σ̂θϕ =
2(1 + ν)
γ0θϕ +
1 + ζ/R
with ǫ0θ, ǫ
θϕ, κ
ϕ and τ
0 defined by equations (13) and (14). The expressions for σ̂θζ and σ̂ϕζ will
not be needed.
We now integrate the stress distributions across the thickness h of the shell (see Figure 1). The stress
resultants and couples obtained in this way are defined per unit of arc length on the reference surface:
∫ h/2
σ̂ii (1 + ζ/R) dζ (i = θ, ϕ) ,
Nθϕ = Nϕθ =
∫ h/2
σ̂θϕ (1 + ζ/R) dζ ,
∫ h/2
σ̂iζ (1 + ζ/R) dζ (i = θ, ϕ) , (17)
∫ h/2
σ̂ii (1 + ζ/R) ζ dζ (i = θ, ϕ) ,
Mθϕ = Mϕθ =
∫ h/2
σ̂θϕ (1 + ζ/R) ζ dζ .
We evaluate these integrals with the expansion (16). The tangential stress resultants are
Nθ = K
ǫ0θ + ν ǫ
Nϕ = K
ǫ0ϕ + ν ǫ
, (18)
Nθϕ = K
γ0θϕ .
Explicit expressions for the transverse shearing stress resultants Qi are not needed since these quantities
will be determined from the equilibrium equations. The moment resultants are
Mθ = D
κ0θ + ν κ
ǫ0θ + ν ǫ
Mϕ = D
κ0ϕ + ν κ
ǫ0ϕ + ν ǫ
, (19)
Mθϕ = D
The extensional rigidity K and the bending rigidity D are defined by
1− ν2
, (20)
12(1− ν2)
. (21)
Their dimensionless ratio ξ is a large number,
ξ = R2
, (22)
the inverse of which will serve as an expansion parameter for thin shell theory.
2.2.4 Equilibrium equations
With the metric (8), the components θ, ϕ and ζ of the equilibrium equations (7) respectively become
(R+ ζ)
(sin θ σ̂θθ),θ + σ̂θϕ,ϕ − cos θ σ̂ϕϕ + sin θ σ̂ζθ
+ sin θ
(R+ ζ)
= 0 ,
(R + ζ)
(sin θ σ̂θϕ),θ + σ̂ϕϕ,ϕ + cos θ σ̂θϕ + sin θ σ̂ζϕ
+ sin θ
(R+ ζ)
= 0 ,
(R + ζ)
(sin θ σ̂θζ),θ + σ̂ϕζ,ϕ − sin θ (σ̂θθ + σ̂ϕϕ)
+ sin θ
(R+ ζ)
= 0 ,
where the equations have been multiplied by sin θ (R + ζ), (R + ζ) and sin θ (R + ζ)2, respectively. The
stress components are given in a normalized basis.
The integration on ζ of these three equations in the range [−h/2, h/2] yields the equilibrium equations
for the forces:
(sin θ Nθ),θ +Nθϕ,ϕ − cos θ Nϕ + sin θ Qθ +R sin θ qθ = 0 , (23)
(sin θNθϕ),θ +Nϕ,ϕ + cos θ Nθϕ + sin θ Qϕ +R sin θ qϕ = 0 , (24)
(sin θ Qθ),θ +Qϕ,ϕ − sin θ (Nθ +Nϕ)−R sin θ q = 0 , (25)
where qθ and qϕ are the components of the tangential load vector per unit area of the reference surface:
(R+ ζ)
= R2 qi (i = θ, ϕ) .
We choose the convention that tensile stresses are positive (see Figure 1). The transverse load per unit
area of the reference surface is noted q and is taken to be positive toward the center of the sphere:
(R+ ζ)
= −R2 q . (26)
The first two equilibrium equations for the stresses can also be multiplied by ζ before the integration
to yield the equilibrium equations for the moments:
(sin θMθ),θ +Mθϕ,ϕ − cos θMϕ −R sin θ Qθ = 0 , (27)
(sin θMθϕ),θ +Mϕ,ϕ + cos θMθϕ −R sin θ Qϕ = 0 . (28)
We have neglected small terms in
ζ (R+ ζ)
where i = (θ, ϕ). A third equilibrium equation for
the moments exists but has the form of an identity: Mθϕ =Mϕθ.
3 Resolution
3.1 Available methods
At this stage the elastic theory for a thin spherical shell involves 17 equations: six strain-displacement
relationships (13)-14), six stress-strain relations (18)-(19) making Hooke’s law, and five equilibrium equa-
tions (23)-(23) and (27)-(28). There are 17 unknowns: six strain components (ǫ0θ, ǫ
θϕ, κ
three displacements (w, vθ , vϕ), three tangential stress resultants (Nθ, Nϕ, Nθϕ), two transverse shearing
stress resultants (Qθ, Qϕ), and three moment resultants (Mθ,Mϕ,Mθϕ). The three equations (16) are
also needed if the tangential stresses (σ̂θθ, σ̂θϕ, σ̂ϕϕ) are required. The quantities of primary interest to
us are the transverse deflection and the tangential stresses (sometimes tangential strain is preferred, as
in Sandwell et al. [1997] or Banerdt and Golombek [2000]). We thus want to find the minimum set of
equations that must be solved to determine these quantities.
We are aware of two methods of resolution [Novozhilov , 1964, p. 66]. In the first one, we insert the
strain-displacement relationships into Hooke’s law, and substitute in turn Hooke’s law into the equilibrium
equations. This method yields three simultaneous differential equations for the displacements. Once the
displacements are known, it is possible to compute the strains and the stresses. The second method
supplements the equilibrium equations with the equations of compatibility [Novozhilov , 1964, p. 27] that
relate the partial derivatives of the strain components. It is then convenient to introduce the so-called
stress functions [Kraus , 1967, p. 243], without direct physical interpretation, which serve to define the
stress resultants without introducing the tangential displacements. Equations relating the transverse
displacement and the stress functions are then found by applying the third equation of equilibrium and
the third equation of compatibility (it is also possible to use all three equations of compatibility in order
to directly solve for the stress and moment resultants). Once the transverse displacement and the stress
functions are known, stresses can be computed.
If the shell thickness is constant, the deformations of a thin spherical shell can be completely calcu-
lated with both methods. If the shell thickness is variable, the three equations governing displacements,
obtained with the first method, cannot be decoupled and are not easy to solve. Kraus’ method with
stress functions leads to a system of three equations (relating the transverse displacement and the two
stress functions), in which the first equation is decoupled and solved before the other two. This method
thus provides a system of equations much easier to solve and will be chosen in this article.
When solving the equations, one usually assumes from the beginning the large ξ limit, i.e. 1+ ξ ∼= ξ
where ξ is defined by equation (22). We will only take this limit at the end of the resolution. This
procedure will not complicate the computations, since we have to compute anyway many new terms
because of the variable shell thickness.
3.2 Differential operators
We will repeatedly encounter the operators Oi which intervene in the expressions (14) for the flexural
strains (κ0θ, κ
0). Since we are looking for scalar equations, we need to find out how the operators Oi
can be combined in order to yield scalar expressions, i.e. expressions that are invariant with respect to
changes of coordinates on the sphere.
The first thing is to relate the Oi to tensorial operators. Starting from the covariant derivatives on the
sphere ∇i, we construct the following tensorial differential operators of the second degree in derivatives:
Dij = ∇i∇j + gij , (29)
where ∇i denotes the covariant derivative (see Appendix 7.2). These operators give zero when applied
on spherical harmonics of degree one (considered as scalars):
Dij Y1m = 0 (m = −1, 0, 1) . (30)
This property can be explicitly checked on the spherical harmonics (109) with the metric and the formulas
for the double covariant derivatives given in Appendix 7.4.
In two-dimensional spherical coordinates (θ, ϕ), the three operators Oi defined by equations (15) are
related to the operators Dij acting on a scalar function f through
O1 f = Dθθ f ,
O2 f = csc2 θDϕϕ f , (31)
O3 f = csc θDθϕ f .
The operators O1,2,3f actually correspond to normalized Dijf , i.e. Dijf/(
giigjj).
The usual derivatives of the operators Oi satisfy the useful identities (112)-(113) which are proven
in Appendix 7.7. These identities are the consequence of the path dependence of the parallel transport
of vectors on the curved surface of the sphere.
Invariant expressions are built by contracting all indices of the differential operators in their tensorial
form. The indices can be contracted with the inverse metric gij or with the antisymmetric tensor εij (see
Appendix 7.4), which should not be confused with the strain tensor ǫij . In the following, a and b are
scalar functions on the sphere.
At degree 2, the only non-zero contraction of the Dij is related to the Laplacian (104):
∆′a ≡ gij Dij a
= (∆+ 2) a (32)
= (O1 +O2) a
= a,θ,θ + cot θ a,θ + csc
2 θ a,ϕ,ϕ + 2 a .
At degree 4, a scalar expression symmetric in (a, b) is given by
A(a ; b) ≡ [∆′a][∆′b]− [Dij a][Dij b] ,
= [∆ a][∆ b]− [∇i∇j a][∇i∇j b] + [∆ a] b+ a [∆ b] + 2 a b (33)
= [O1 a][O2 b] + [O2 a][O1 b]− 2 [O3 a][O3 b]
= (a,θ,θ + a)
csc2 θ b,ϕ,ϕ + cot θ b,θ + b
csc2 θ a,ϕ,ϕ + cot θ a,θ + a
(b,θ,θ + b)
−2 csc2 θ (a,θ,ϕ − cot θ a,ϕ) (b,θ,ϕ − cot θ b,ϕ) .
where upper indices are raised with the inverse metric: Dij = gikgjlDkl. The action of an operator does
not extend beyond the brackets enclosing it.
If a is constant, A(a ; b) = a∆′b. It is useful to define an associated operator A0 that gives zero if
its first argument is constant:
A0(a ; b) = A(a ; b)− a [∆′b] . (34)
A scalar expression of degree 4 antisymmetric in (a, b) is given by
B1(a ; b) ≡ gij εkl [Dik a] [Djl b]
= gij εkl [∇i∇k a] [∇j∇l b] (35)
= [(O1 −O2) a] [O3b]− [O3a] [(O1 −O2) b]
= csc θ
a,θ,θ − csc2 θ a,ϕ,ϕ − cot θ a,θ
(b,θ,ϕ − cot θ b,ϕ)
− csc θ (a,θ,ϕ − cot θ a,ϕ)
b,θ,θ − csc2 θ b,ϕ,ϕ − cot θ b,θ
We will also need another operator of degree 4:
B2(a ; b) ≡ εij [∇ia ] [∇j ∆′b]
= csc θ
a,θ [∆
′b],ϕ − a,ϕ [∆
′b],θ
. (36)
The sum of the operators B1 and B2 is noted B:
B(a ; b) = B1(a ; b) + B2(a ; b) . (37)
If either a or b is constant, B(a ; b) = 0.
The operators A and B have an interesting property proven in Appendix 7.8: for arbitrary scalar
functions a and b, A(a ; b) and B(a ; b) do not have a degree-one term in their spherical harmonic expansion.
This is not true of A0, B1 and B2.
3.3 Transverse displacement
3.3.1 Resolution of the equations of equilibrium
The first step consists in finding expressions for the moment resultants (Mθ,Mϕ,Mθϕ) in terms of the
transverse displacement and the stress resultants. The extensional strains (ǫ0θ, ǫ
θϕ) can be eliminated
from the equations for stress and moment resultants (18)-(19). The flexural strains (κ0θ, κ
0) depend
on the transverse displacement w through equations (14). We thus obtain
Mθ = −
(O1 + νO2)w +
Mϕ = −
(O2 + νO1)w +
Nϕ , (38)
Mθϕ = −
(1− ν)O3 w +
Nθϕ ,
where the parameter ξ is defined by equation (22).
The second step consists in solving the equilibrium equations for moments in order to find the
transverse shearing stress resultants (Qθ, Qϕ). We substitute expressions (38) into equations (27)-(28).
Knowing that the stress resultants satisfy the equilibrium equations (23)-(24), we obtain new expressions
for Qθ and Qϕ (identities (112)-(113) are helpful):
Qθ = −
(D∆′w),θ − (1− η)Rqθ
(1− ν) (D,θ O2 − csc θD,ϕO3) w −
(η,θNθ + csc θ η,ϕNθϕ) , (39)
Qϕ = −
csc θ (D∆′w),ϕ − (1− η)Rqϕ
(1− ν) (−D,θ O3 + csc θD,ϕ O1) w −
(η,θ Nθϕ + csc θ η,ϕNϕ) . (40)
The operator ∆′ is defined by equation (32) and η is a parameter close to 1:
1 + ξ
. (41)
The third step consists in finding expressions for the stress resultants (Nθ, Nϕ, Nθϕ) in terms of
stress functions by solving the first two equilibrium equations (23)-(24). Let us define the following linear
combinations of the stress and moment resultants:
(Pθ, Pϕ, Pθϕ) =
Mθ, Nϕ +
Mϕ, Nθϕ +
. (42)
We observe that these linear combinations satisfy simplified equations of equilibrium:
(sin θ Pθ),θ + Pθϕ,ϕ − cos θ Pϕ +R sin θ qθ = 0 ,
(sin θ Pθϕ),θ + Pϕ,ϕ + cos θ Pθϕ +R sin θ qϕ = 0 . (43)
Comparing these equilibrium equations with the identities (112)-(113), we see that the homogeneous
equations (i.e. equations (43) with a zero tangential load qθ = qϕ = 0) are always satisfied if
(Pθ, Pϕ, Pθϕ) = (O2,O1,−O3)F , (44)
where F is an auxiliary function called stress function. For the moment, this function is completely
arbitrary apart from being scalar and differentiable.
Particular solutions of the full equations (43) can be found if we express the tangential load qT =
qθθ̂ + qϕϕ̂ in terms of the surface gradient of a scalar potential Ω (consoidal or poloidal component) and
the surface curl of a vector potential V r̂ (toroidal component):
qT = −
∇̄Ω + 1
∇̄ × (V r̂) . (45)
Surface operators are defined in Appendix 7.5, where the terms consoidal/poloidal are also discussed.
The covariant components of qT are (qθ, sin θqϕ) and can be expressed as − 1RΩ,i +
gjkεikV,j , which
gives
qθ = −
Ω,θ +
R sin θ
V,ϕ ,
sin θ qϕ = −
Ω,ϕ −
sin θ
V,θ . (46)
If the tangential load is consoidal (V = 0), a particular solution of equations (43) is given by
(Pθ, Pϕ, Pθϕ) = (1, 1, 0)Ω . (47)
If the tangential load is toroidal (Ω = 0), a particular solution of equations (43) is given by
(Pθ, Pϕ, Pθϕ) = (2O3,−2O3,O2 −O1)H , (48)
where we have introduced a second stress function H which satisfies the constraint
∆′H = −V + V0 , (49)
where V0 is a constant (identities (112)-(113) are useful). This equation allows us to determine the stress
function H if the toroidal source V is known.
The general solution of the equations (43) is given by the sum of the general solution (44) of the
homogeneous equations and the two particular solutions (47)-(48) of the full equations:
(Pθ, Pϕ, Pθϕ) = (O2F +Ω + 2O3H,O1F +Ω− 2O3H,−O3F + (O2 −O1)H) . (50)
The stress resultants (Nθ, Nϕ, Nθϕ) can now be obtained from (Pθ, Pϕ, Pθϕ) by using equations (42) and
(38):
(Nθ, Nϕ, Nθϕ) = η (Pθ, Pϕ, Pθϕ) + η
(O1 + νO2,O2 + νO1, (1− ν)O3)w , (51)
which finally give
Nθ = η
O2 F +
(∆′ − (1− ν)O2)w +Ω+ 2O3H
Nϕ = η
O1 F +
(∆′ − (1− ν)O1)w +Ω− 2O3H
, (52)
Nθϕ = η
−O3 F + (1 − ν)
O3 w + (O2 −O1)H
The fourth step consists in expressing the third equation of equilibrium (25) in terms of the transverse
displacement w and the stress functions F and H . For this purpose, it is handy to express the transverse
shearing stress resultants (Qθ, Qϕ) in terms of (w,F,H). We thus substitute Nϕ, Nθ, Nθϕ, given by
equations (52), into the expressions for Qθ and Qϕ, given by equations (39)-(40):
Qθ = −
(ηD∆′w),θ +
(ηD),θ O2 − csc θ (ηD),ϕ O3
− (η,θ O2 − csc θ η,ϕO3)F + Ω,θ − (ηΩ),θ − csc θ V,ϕ
− 2 (η,θ O3 + csc θ η,ϕ O2)H + csc θ (η,ϕ ∆′H − η(∆′H),ϕ) , (53)
Qϕ = −
csc θ (ηD∆′w),ϕ +
− (ηD),θ O3 + csc θ (ηD),ϕ O1
− (−η,θ O3 + csc θ η,ϕ O1)F + csc θ (Ω,ϕ − (ηΩ),ϕ) + V,θ
+2 (η,θ O1 + csc θ η,ϕ O3)H − (η,θ ∆′H − η(∆′H),θ) . (54)
The various terms present in Qθ and Qϕ can be classified into generic types according to their
differential structure. Each generic type will contribute in a characteristic way to the third equation
of equilibrium. It is worthwhile to compute the generic contribution of each type before using the full
expressions of the stress resultants with their multiple terms. Identities (112)-(113) are helpful in this
calculation. We thus write the third equation of equilibrium (25) as follows:
I (Qθ ;Qϕ ; 0)− (Nθ +Nϕ)−Rq = 0 , (55)
where the differential operator I is defined for arbitrary expressions (X,Y, Z) by
I(X ;Y ;Z) = csc θ (sin θX),θ + csc θ Y,ϕ − cot θ Z . (56)
This operator must be evaluated for the following generic types present in (Qθ, Qϕ):
I (a,θ ; csc θ a,ϕ ; 0) = ∆a ,
I (− csc θ a,ϕ ; a,θ ; 0) = 0 ,
I (a,θ O2b− csc θ a,ϕ O3b ;−a,θ O3b+ csc θ a,ϕO1b ; 0) = A0(a ; b) , (57)
I (a,θ O3b+ csc θ a,ϕ O2b ;−a,θ O1b− csc θ a,ϕO3b ; 0) = B1(a ; b) ,
I (− csc θ (a,ϕ∆′b − a[∆′b],ϕ) ; a,θ ∆′b− a[∆′b],θ ; 0) = 2B2(a ; b) ,
where a and b are scalar functions. The operators A0, B1 and B2 are defined in section 3.2.
We now substitute (Qθ, Qϕ) and (Nϕ, Nθ) into the third equation of equilibrium (55) and use formulas
(57). We thus obtain the first of the differential equations that relate w and the stress functions F and
∆′ (ηD∆′w)− (1 − ν)A(ηD ;w) +R3 A(η ;F ) + 2R3 B(η ;H) = −R4 q +R3 (∆Ω−∆′(ηΩ)) . (58)
The operators ∆′, A and B are defined in section 3.2.
3.3.2 Compatibility relation
A second equation relating the transverse displacement and the stress functions comes from the compat-
ibility relation which is derived by eliminating (vθ, vϕ) from the strain-displacement equations (13):
sin θ γ0θϕ,ϕ
sin2 θ ǫ0ϕ,θ
+ ǫ0θ,ϕ,ϕ − sin θ cos θ ǫ0θ,θ + 2 sin2 θ ǫ0θ −
sin2 θ
∆′w . (59)
The strain components are related to the stress resultants through Hooke’s law (18). The substitution of
equations (18) into the compatibility equation (59) gives
∆′ (α (Nθ +Nϕ))−
∆′w − (1 + ν)J (αNθ ;αNϕ ;αNθϕ) = 0 . (60)
where α is the reciprocal of the extensional rigidity:
α ≡ 1
K(1− ν2)
. (61)
For arbitrary expressions (X,Y, Z), J (X ;Y ;Z) is defined by
J (X ;Y ;Z) = csc2 θ
sin2 θX,θ
+ Y,ϕ,ϕ + 2 (sin θ Z,ϕ),θ
− cot θ Y,θ + 2 Y . (62)
As in the case of the third equation of equilibrium, it is practical to classify the terms present in
(Nθ, Nϕ, Nθϕ) into generic types and evaluate separately their contribution to the equation of compati-
bility. There are three types of terms for which we must evaluate the operator J (identities (112)-(115)
are helpful in this calculation):
J (a ; a ; 0) = ∆′a ,
J (aO2b ; aO1b ;−aO3b) = A(a ; b) , (63)
J (2aO3b ;−2aO3b ; a (O2 −O1)b) = 2B(a ; b) ,
where a and b are scalar functions. The operators ∆′, A and B are defined in section 3.2.
We now evaluate J in equation (60) with expressions (52) and formulas (63):
J (αNθ ;αNϕ ;αNθϕ) = A(ηα ;F ) +
∆′ (ηαD∆′w)− 1
(1− ν)A(ηαD ;w)
+ ∆′ (ηαΩ) + 2B(ηα ;H) . (64)
The term A(ηαD ;w) can be rewritten with the help of the following equality:
1− ν2
A(ηαD ;w) = ∆′w −A(η ;w) , (65)
since (1 − ν2)ηαD/R2 = η/ξ and (η/ξ),i = −η,i.
We finally substitute Nθ + Nϕ, given by equations (52), into the compatibility equation (60) and
use expressions (64)-(65). We thus obtain the second of the differential equations that relate w and the
stress functions F and H :
∆′ (ηα∆′F )− (1 + ν)A(ηα ;F )−
A(η ;w) − 2(1 + ν)B(ηα ;H) = −(1− ν)∆′(ηαΩ) . (66)
3.4 Tangential displacements
Assuming that the flexure equations (58) and (66) for the transverse displacement w and the stress
functions (F,H) have been solved, we now show how to calculate the tangential displacements.
In analogy with the decomposition of the tangential load in equations (45)-(46), the tangential
displacement can be separated into consoidal and toroidal components:
v = ∇̄S + ∇̄ × (T r̂) , (67)
where S and T are the consoidal and toroidal scalars, respectively. The covariant components of v are
(vθ, sin θ vϕ) and can be expressed as S,i + g
jkεikT,j (see Appendix 7.5), which gives
vθ = S,θ + csc θ T,ϕ ,
sin θ vϕ = S,ϕ − sin θ T,θ . (68)
The strain-displacement equations (13) become
ǫ0θ =
((O1 − 1)S +O3 T + w) ,
ǫ0ϕ =
((O2 − 1)S −O3 T + w) , (69)
γ0θϕ =
(2O3S + (O2 −O1)T ) .
The stress resultants (Nθ, Nϕ, Nθϕ) given by equations (18) become
(∆S + (1 + ν)w − (1− ν) ((O2 − 1)S −O3 T )) ,
(∆S + (1 + ν)w − (1− ν) ((O1 − 1)S +O3 T )) , (70)
Nθϕ =
(2O3 S + (O2 −O1)T ) .
The toroidal potential T cancels in the sum Nθ +Nϕ:
Nθ +Nϕ =
(1 + ν) (∆S + 2w) .
The consoidal displacement potential S can thus be related to (w,F,Ω) by using expressions (52) for the
stress resultants:
∆S = Rηα (1− ν) (∆′F + 2Ω) + η
∆′w − 2w , (71)
where α is defined by equation (61).
It is more difficult to extract the toroidal displacement potential T . When the shell thickness is
constant, decoupled equations for the displacements can be found by suitable differentiation and combi-
nation of the three equilibrium equations (23)-(25) for the stress resultants. This method does not work
if the shell thickness is variable because the resulting equations are coupled. The trick consists in relating
the tangential displacements to (w,F,H) by the way of equations similar to the homogeneous part of the
first two equilibrium equations, but with (Nθ, Nϕ, Nθϕ) replaced by
(Nθ, Nϕ, Nθϕ), so that derivatives
do not mix with derivatives of T . We will thus calculate the following expression in two different
ways (from equations (52) and (70)):
= csc2 θ (−X,ϕ + sin θ Y,θ) ,
where X and Y are defined by
X = sin θ I
Nθϕ ;
Y = sin θ I
Nθϕ ;
Nϕ ;−
The operator I is defined by equation (56).
As before, it is easier to begin with the evaluation of the operator Z for generic contributions:
Z(a ; a ; 0) = 0 ,
Z(aO2b ; aO1b ;−aO3b) = −B(a ; b) ,
Z(2aO3b ;−2aO3b ; a(O2 −O1)b) = 2A(a ; b)−∆′ (a∆′b) ,
On the one hand, the evaluation Z for (Nθ, Nϕ, Nθϕ) given by equations (70) gives
−1− ν
∆∆′ T .
On the other hand, the evaluation of Z for (Nθ, Nϕ, Nθϕ) given by equations (52) gives
+ (1− ν)B
+ 2RA
The equality of the two previous formulas yields the sought equation for T :
∆∆′ T = 2R (1 + ν) (B (ηα ;F )− 2A (ηα ;H) + ∆′ (ηα∆′H)) + 2B (η ;w) . (72)
We have used the relation (η/ξ),i = −η,i.
The equation for T shows that toroidal displacement always occurs when the shell thickness varies.
The right-hand side of equation (72) only vanishes when two conditions are met: (1) there is no toroidal
source (so that the stress function H vanishes) and (2) the shell thickness is constant (so that the terms
in B vanish).
3.5 Stresses
Stresses can be computed from (w,F,H) and Ω by substituting equations (14), (18) and (52) into equa-
tions (16):
σ̂θθ =
(O2 F +Ω+ 2O3H) +
R(1− ν2)
(∆′ − (1− ν)O2)w ,
σ̂ϕϕ =
(O1 F +Ω− 2O3H) +
R(1− ν2)
(∆′ − (1− ν)O1)w , (73)
σ̂θϕ =
(−O3 F + (O2 −O1)H) +
R(1 + ν)
O3w .
Stresses at the surface are obtained by setting ζ = h/2.
4 Flexure equations and their properties
4.1 Thin shell approximation
The flexure equations derived in section 3 already include several assumptions of the thin shell theory,
but not yet the first one stating that the shell is thin. Of course, the three other assumptions can be
seen to be consequences of the first one [see Kraus, 1967, chap. 2.2], but we have not imposed in a
quantitative way the thinness condition on the equations. We thus impose the limit of small h/R or,
equivalently, the limit of large ξ (defined by equation (22)) on the flexure equations for (H,F,w, S, T ).
This procedure amounts to expand η ≈ 1 − 1/ξ (neglecting terms in 1/ξ wherever appropriate) and to
neglect the derivatives of η in equations (49), (58), (66), (71) and (72). We thus obtain the final flexure
equations for the displacements (w, S, T ) and for the stress functions (F,H):
∆′H = −V + V0 , (74)
∆′ (D∆′w)− (1 − ν)A(D ;w) +R3∆′F = −R4 q − 2R3Ω+ R
∆Ω , (75)
∆′ (α∆′F )− (1 + ν)A(α ;F )− 1
∆′w = −(1− ν)∆′(αΩ) + 2(1 + ν)B(α ;H) , (76)
∆S = Rα (1 − ν) (∆′F + 2Ω) + 1
∆′w − 2w , (77)
∆∆′ T = 2R (1 + ν) (B (α ;F )− 2A (α ;H) + ∆′ (α∆′H)) . (78)
Recall that the differential operators ∆′, A and B are defined by equations (32), (33) and (37). The
potential pairs (Ω, V ) and (S, T ) are related to the tangential load and displacement by equations (45)
and (67), respectively. In the second equation, the term R3∆Ω/ξ has been kept since it could be large if
Ω has a short wavelength. For the same reason, the term ∆′w/ξ has been kept in the fourth equation.
In the third and fifth equations, the terms depending on H belong to the right-hand sides since they
can be considered as a source once H has been calculated from the first equation. The same can be said
of the terms depending on w and F in the fourth and fifth equations: they are supposed to be known from
the simultaneous resolution of the second and third equations. The difficulty in solving the equations
thus lies with the two flexure equations for w and F which are linear with non-constant coefficients (all
other equations are linear - in their unknowns - with constant coefficients). Once these two core equations
have been solved, all other quantities are easily derived from them.
The bending rigidity D, defined by equation (21), characterizes the bending regime: the shell locally
bends in a similar way as a flat plate undergoing small deflections with negligible stretching. The pa-
rameter α, defined by equation (61), is the reciprocal of the extensional rigidity K and characterizes the
membrane regime of the shell in which bending moments are negligible and the load is mainly supported
by internal stresses tangent to the shell. Since D and K are respectively proportional to the third and
first power of the shell thickness, the membrane regime (in which the D-depending terms are neglected)
is obtained in the limit of an extremely thin shell. This observation and the fact that such a shell, lacking
rigidity, cannot support bending moments justify the use of the term ‘membrane’.
The weight of the various terms in the equations depends on two competing factors: the magnitude
of the coefficient multiplying the derivative and the number of derivatives. On the one hand, a coefficient
containing D will be smaller than a coefficient containing α−1 since αD/R2 ∼ ξ−1 is a small number
(see equation (22)). On the other hand, a large number of derivatives will increase the weight of the term
if the derived function has a small wavelength. The transition between the membrane and the bending
regimes thus depends on the wavelength of the load: if the load has a large wavelength (with respect to
the shell radius), the flexure of the shell will be well described by equations without the terms depending
on D (see section 5.1), whereas the flexure under loads of small wavelength is well described by equations
keeping only the D-depending terms with the largest number of derivatives (see section 5.2).
Stress functions are associated with membrane stretching and give negligible contributions in the
bending regime. Formulas (52) show that F (respectively H) plays the role of potential for the stress
resultants in the membrane regime when the load is transversal (respectively tangential toroidal). There
is no stress function associated with the tangential consoidal component of the load because it is identical
to Ω, the consoidal potential of the load.
The variation of the shell thickness has two effects. First, it couples the spherical harmonic modes
that are solutions to the equations for a shell of constant thickness. Second, the toroidal part remains
intertwined with the transversal and consoidal parts, whereas it decouples if the thickness is constant.
For example, the toroidal load is a source for the transverse deflection through the term B in the third
flexure equation. Furthermore, the stress function F is a source for the toroidal displacement in the fifth
flexure equation.
For numerical computations, the flexure equations (74)-(78) obtained in the thin shell approximation
(to which we can add the equations (73) for the stresses) are adequate. For this purpose, it is not useful
to keep small terms in 1/ξ since the theory rests on assumptions only true for a thin shell. In the rest of
the article, we will continue to work with the equations (49), (58), (66), (71) and (72) for more generality.
4.2 Covariance
Because of their tensorial form, the flexure equations (74)-(78) are covariant; this is also true of the more
general equations (49), (58), (66), (71) and (72). This property means that their form is valid in all
systems of coordinates on the sphere, though the tensor components (the covariant derivatives of scalar
functions, the metric and the antisymmetric tensor) will have a different expression in each system. The
scalar functions will also have a different dependence on the coordinates in each system. The covariance
of the final equations was expected. We indeed started with tensorial equations in section 2.1; their
restriction to the sphere in principle respected the tensoriality with respect to changes of coordinates
on the sphere. However the covariance of the two-dimensional theory was not made explicit until we
obtained the final equations. This property is thus a strong constraint on the form of the solution and
a check of its validity, though only necessary and not sufficient (other covariant terms may have been
ignored).
Another advantage of the covariant form is the facility to express the final equations in different
systems of coordinates (even non-orthogonal ones) with the aim of solving them. For example, the fi-
nite difference method in spherical coordinates (θ, ϕ) suffers from a very irregular grid and from pole
singularities. These problems can be avoided with the ‘cubed sphere’ coordinate system [Ronchi et al.,
1996]. Operators including covariant derivatives (here the Laplacian and the operators A and B) can ex-
pressed in any system of coordinates whose metric is known. Christoffel symbols and tensorial differential
operators can be computed with symbolic mathematical software.
4.3 Degree one
Displacements of degree one require special consideration. All differential operators acting on (w,F,H)
(that is ∆′, A and B) in the flexure equations (58) and (66) can be expressed in terms of the operators
Dij = ∇i∇j+gij (see section 3.2). The Dij have the interesting property that they give zero when acting
on spherical harmonics of degree one (see equation (30)). Therefore the degree-one terms in the spherical
harmonic expansion of w vanish from the flexure equations. The magnitude of the transverse deflection
of degree one neither depends on the load nor on the elastic properties of the spherical shell.
More generally, the homogeneous (q = Ω = V = 0) flexure equations (49), (58) and (66) are satisfied
by w and F being both of degree one. According to equations (71)-(72), the corresponding tangential
displacement is constrained by ∆S = −2w and ∆∆′T = 0, so that S and T are also of degree one,
with S = w. These conditions lead to vanishing strains (see equations (69)) which indicate a rigid
displacement. The total displacement is then given by
u = w r̂+ ∇̄w + ∇̄ × (T r̂) ,
where w and T are of degree one. In Appendix 7.6, we show that the first two terms represent a
rigid translation whereas the last term represents a rigid rotation. As expected, stresses vanish for
such displacements (see equation (73)). This freedom in translating or rotating the solution reflects the
freedom in the choice of the reference frame (in practice the reference frame is centered at the center of
the undeformed shell). The same freedom of translation is also found in the theory of deformations of a
spherical, radially stratified, gravitating solid [e.g. Farrell , 1972; Greff-Lefftz and Legros , 1997; Blewitt ,
2003].
What can we say about degree-one loading? Let us first examine what happens when the flexure
equations are projected on the spherical harmonics of degree one. Since the operator ∆′ annihilates the
degree one in any spherical harmonic expansion, terms of the form ∆′f vanish when they are projected
on the spherical harmonics of degree one. Moreover, the operators A and B also vanish in this projection
since they do not contain a degree-one term (see equations (116)-(117) of Appendix 7.8). Therefore the
degree-one component of the flexure equation (66) is identically zero, whereas the degree-one component
of the flexure equation (58) is
q + ∇̄ · qT
Yi = 0 (i = x, y, z) , (79)
where dω = sin θ dθ dϕ and Yi are the real spherical harmonics of degree one. The integral is taken over
the whole spherical surface. We have used the relation ∆Ω = −R ∇̄ ·qT derived from equations (45) and
(106).
The first term in the integrand of equation (79) is the projection on the Cartesian axes (x̂, ŷ, ẑ) of
the vector field q r̂:
(q Yx, q Yy, q Yz) = (q r̂ · x̂, q r̂ · ŷ, q r̂ · ẑ) ,
where we have used formulas (109) for the spherical harmonics.
The second term in equation (79) can be rewritten with the identity (102) and Gauss’ theorem (107):
∇̄ · qT
Yi = −
dω qT · ∇̄Yi (i = x, y, z) ,
where the Yi are considered as scalars. Since qT is orthogonal to r̂, the integrand is the projection of qT
on the Cartesian axes (x̂, ŷ, ẑ):
qT · ∇̄Yx,qT · ∇̄Yy,qT · ∇̄Yz
= (qT · x̂,qT · ŷ,qT · ẑ) ,
where we have used formulas (110) for the gradients of the spherical harmonics.
Recalling that the transverse load q was defined positive towards the center of the sphere (see equation
(26)), we define a total load vector q = −q r̂+qT . With the above results, we can rewrite the degree-one
projection (79) as
dω (q · x̂ ,q · ŷ ,q · ẑ) = (0, 0, 0) , (80)
which means that the integral (over the whole spherical surface) of the projection on the coordinate axes
of the total load vector vanishes. This result is the consequence of the static assumption in the equations
of motion (7), since a non-zero sum of the external forces would accelerate the sphere.
In practice, degree-one loads on planetary surfaces are essentially due to mass redistribution [Greff-Lefftz and Legros ,
1997] and have a tangential consoidal component (for example the gravitational force is not directed to-
ward the center of figure of the shell). If the shell thickness is variable, a non-zero Ω of degree one will
induce degrees higher than one in w and in S. If the shell thickness is constant, the degree-one load
drops from the flexure equations (58)-(66) and w is not affected. However, the degree-one Ω generates
(assuming a constant shell thickness) a degree-one tangential displacement through equation (71), so that
S 6= w. Whether the shell thickness is variable or not, a degree-one Ω thus generates a total displacement
which is not only a translation but also a tangential deformation [Blewitt , 2003], in which case stresses
do not vanish as shown by equations (73).
5 Limit cases
5.1 Membrane limit
A shell is in a membrane state of stress if bending moments (Mθ,Mϕ,Mθϕ) can be neglected, in analogy
with a membrane which cannot support bending moments. Equations (19) show that this is true if the
bending rigidity vanishes: D = 0. Consistency with equation (22) imposes the limit of infinite ξ (or
η = 1). With these approximations, the first flexure equation (58) for (w,F,H) becomes
∆′F = −Rq − 2Ω . (81)
The second flexure equation (66) for (w,F,H) becomes
∆′w = ∆′ (α∆′F )− (1 + ν)A(α ;F )− 2(1 + ν)B(α ;H) + (1− ν)∆′(αΩ) . (82)
If q is independent of w, F and w can be successively determined with spherical harmonic transforms
from equations (81) and (82) (though the right-hand side of the latter equation must be computed with
another method). If q has a linear dependence in w (such as when the sphere is filled with a fluid), w
can be eliminated between equations (81) and (82), so that F and w can also be computed in succession
(however the equation for F cannot be solved by a spherical harmonic transform). H is supposed to be
known since equation (49) is not modified and can be solved with spherical harmonics.
The stresses are obtained from equations (73) with the additional approximation of neglecting the
term in ζ:
(σ̂θθ, σ̂ϕϕ, σ̂θϕ) =
(O2F +Ω+ 2O3H,O1F +Ω− 2O3H,−O3F + (O2 −O1)H) .
Bending moments play a small role if the load has a large wavelength. In practice, the threshold at
which bending moments become significant can be evaluated from the constant thickness equation for w
(see section 5.3.1). One must compare the magnitudes of the terms in D and 1/α in the left-hand side
of equation (88). Bending moments are negligible (i.e. the term in D) if the spherical harmonic degree ℓ
of the transverse displacement w is such that
ℓ . k
, (83)
where k = (12(1− ν2))1/4 ≈ 1.8 if we take ν = 1/4. This threshold is about 10 for a planet with a radius
of 3400 km and a lithospheric thickness equal to 100 km. Flexure equations in the membrane limit of a
shell with constant thickness have been used by Sleep and Phillips [1985] to study the lithospheric stress
in the Tharsis region of the planet Mars.
5.2 Euclidean limit
The equation for the deflection of a rectangular plate with variable thickness has been derived by
Timoshenko and Woinowsky-Krieger [1964, p. 173]. We will check here that the Euclidean limit of
our equations gives the same answer.
Let us define the coordinates (x, y) = (Rϕ,R θ′), where θ′ = θ/2 − θ is the latitude. We work in a
small latitude band around the equator (so that θ′ ≪ 1) and in the limit of large spherical radius R and
large ξ (η = 1). Under this change of coordinate, each derivative introduces a factor R so that terms with
the largest number of derivatives dominate. In particular, covariant derivatives can be approximated by
usual derivatives. The surface Laplacian (104) can be approximated as follows:
∆ ≈ R2
≡ R2∆e .
We assume that there are no tangential loads (Ω = V = H = 0). The flexure equations (75)-(76) for the
transverse displacement become:
R4∆e (D∆ew)− (1 − ν)R4 Ae(D ;w) +R5∆eF = −R4 q , (84)
R4 ∆e (α∆eF )− (1 + ν)R4 Ae(α ;F )−R ∆ew = 0 , (85)
where the operator Ae is defined by
Ae(a ; b) = (∆e a) (∆e b)− a,i,j b,i,j
Equation (85) gives a relation between the magnitudes of F and w:
O(R3 ∆eF ) ∼ O(w/α) .
In the large R limit, equation (84) thus becomes
∆e (D∆ew) − (1− ν)Ae(D ;w) = −q ,
which is the equation derived by Timoshenko and Woinowsky-Krieger [1964]. This equation has been
used by Stark et al. [2003], Kirby and Swain [2004] and Pérez-Gussinyé et al. [2004] for the local analysis
of the lithosphere of the Earth. Its one-dimensional version, in which Ae vanishes, has been used by
Sandwell [1984] to describe the flexure of the oceanic lithosphere on Earth and by Stewart and Watts
[1997] to model the flexure at mountain ranges.
5.3 Shell with constant thickness
5.3.1 Displacements
If the thickness of the shell is constant, the toroidal part of the tangential displacement decouples. The
terms in B indeed drop from the flexure equations (58) and (66) so that the equations for (w,F ) depend
only on (q,Ω):
ηD∆′∆′w − (1− ν) ηD∆′w + ηR3 ∆′F = −R4 q +R3 ((1 − η)∆Ω− 2ηΩ) , (86)
∆′∆′F − (1 + ν)∆′F − 1
∆′w = −(1− ν)∆′Ω , (87)
where we used the property A(a ; b) = a∆′b valid for constant a. We eliminate F from these equations
and obtain a sixth order equation relating w to (q,Ω):
ηD∆∆′∆′w +
∆′w = −R4 (∆′ − 1− ν) q +R3
1 + ξ
∆′ − 1− ν
∆Ω . (88)
The elimination of ∆′F between equations (77) and (86) gives an equation relating the consoidal
displacement potential S to (w, q,Ω):
∆S = − 1
1 + ν
1 + ξ
∆∆′w − 2w − (1 − ν)R2α q + 1− ν
1 + ξ
Rα∆Ω .
Terms not including a Laplacian can be eliminated with equation (88), so that we obtain an explicit
solution for S in terms of (w, q,Ω):
1 + ξ
1− ν2
(∆ + 1 + ν)∆′w + w +R2α q − Rα
1 + ξ
(∆− ξ(1 + ν))Ω , (89)
where the integration constant has been set to zero.
Equations (49) and (72) give an equation for the toroidal displacement potential:
∆′ T = −2ηRα (1 + ν)V , (90)
where the integration constant has been set to zero. We have assumed that ∆V 6= 0, otherwise we get
∆′T = 0.
The differential equations given in this section can be solved with spherical harmonics so that the co-
efficients of the spherical harmonic expansions of (w, S, T ) can be expressed in terms of the corresponding
coefficients of the loads (q,Ω, V ) (see Kraus [1967], Turcotte et al. [1981], Banerdt [1986]).
5.3.2 Comparison with the literature
We now compare our equations for a shell of constant thickness with those found in the literature. The
formulas of Banerdt [1986] (taken from the work of Vlasov [1964]) are the most general:
∆3 + 4∆2
(∆ + 2)w = −R4 (∆ + 1− ν) q +R3
∆− 1− ν
∆Ω , (91)
(∆ + 2)χ =
Dξ(1− ν)
∆V . (92)
Banerdt’s notation is slightly different: his formulas are obtained with the substitutions ξ → ψ, Ω → RΩ
and V → RV . The normal rotation χ is proportional to the radial component (in a normalized basis) of
the curl of the tangential displacement:
(∇× v)r̂ .
The curl ∇× v is related to our surface curl (103) by
∇× v = ∇̄ × v + csc θ
(sin θv̂ϕ),θ − v̂θ,ϕ
With the formulas (68) and (104), we get ∇× v = −∆T r̂ so that equation (92) becomes
∆′ T = −2Rα(1 + ν)V . (93)
We see that Banerdt’s equations (91) and (93) coincide with our equations (88) and (90) in the limit
of large ξ (η = 1), with one exception: the bending term for w is written D(∆3 + 4∆2)w instead of
D∆∆′∆′w = (∆3 + 4∆2 + 4∆)w. This error has propagated in many articles and is of consequence for
the degree-one harmonic component, since it violates the static assumption and spoils the translation
invariance discussed in section 4.3. The impact on higher degrees is negligible. Because of this mistake,
many authors give a separate treatment to the first harmonic degree. Banerdt also gives formulas for the
tangential displacements in terms of consoidal and toroidal scalars (A,B) corresponding to our scalars
(S, T ): his formula (A10) is equivalent to our equation (89) in the limit of large ξ.
If we ignore temperature effects, Kraus’ first equation for (w,F ) is equivalent to our equation
(86) in the limit of large ξ, whereas his second equation for (w,F ) is equivalent to the combination
eq.(87)+ 1+ν
eq.(86) in the limit of large ξ [see Kraus , 1967, eq. 6.54h and 6.55d]. Note that the
definition of Kraus’ stress function F [Kraus, 1967, p. 243] differs from ours:
FKraus = F − k(1 − ν)
with k = 1. This freedom of redefining F for arbitrary k remains as long as D is constant. The flexure
equation for w, equation (88), is unaffected so that the solution for w is unchanged. In the final step,
Kraus makes a mistake when combining the two equations for (w,F ) and thus obtains a flexure equation
for w with the same error as in equation (91). Kraus does not include toroidal loading. The flexure
equation of Turcotte et al. [1981] is taken from Kraus [1967] without the tangential loading and is the
same as equation (91) with Ω = 0.
The flexure equation of Brotchie and Silvester [1969] is given in our notation by
D∆2w +
w = −R4 q , (94)
where q includes their term γw describing the response of the enclosed liquid. This equation can be
obtained from our equation (88) as follows: keep only the derivatives of the highest order in each term,
set Ω = 0, take the limit of large ξ (η = 1) and integrate. Brotchie and Silvester choose to work in
the approximation of a shallow shell and with axisymmetrical loading, solving their equation in polar
coordinates with Bessel-Kelvin functions. The reduction to fourth order in equation (94), the shallow
shell approximation and the axisymmetrical assumption are not justified nowadays since the full equation
(88) can be quickly solved with computer-generated spherical harmonics.
The contraction due to a transverse load of degree 0, w = −R2α(1 − ν)q/2, is equivalent to the
radial displacement computed by Love in the limit of a thin shell [Love, 1944, p.142]. However ad-
ditional assumptions about the initial state of stress and the internal density changes are necessary
[Willemann and Turcotte, 1982] so that the degree 0 is usually excluded from the analysis.
5.3.3 Breakdown of the third assumption of thin shell theory
With the spherical harmonic solutions of the equations for a shell of constant thickness, it is possible to
check the thin shell assumption stating that the transverse normal stress is negligible with respect to the
tangential normal stress. The magnitude of the former can be estimated by the load q (see definition
(26)) whereas the magnitude of the latter can be approximated with formulas (73) evaluated on the outer
surface:
(σ̂θθ + σ̂ϕϕ) |h
∆′F − Eh
4R2(1− ν)
∆′w .
where we have assumed the absence of tangential loads (Ω = 0) and the limit of large ξ. We can relate
σT to q by using the solution in spherical harmonics of equations (87) and (88). Since the thin shell
assumption is expected to fail for a load of sufficiently small wavelength, we assume that the spherical
harmonic degree ℓ is large. Assuming ℓ≫ 1, we obtain
(∆′F )ℓm ≈
wℓm ,
wℓm ≈ −
1 + ℓ
ξ(1−ν2)
qℓm ,
where the spherical harmonic coefficients are indexed by their degree ℓ and their order m. If the shell is
not in a membrane state of stress (see equation (83)), ℓ2 > 2R/h so that σT can be approximated by
(σT )lm ≈
ξ(1 + ν)
qlm .
The thin shell assumption holds if q < σT , that is if
3(1 + ν)
or λ >
3(1 + ν)
h , (95)
where λ is the load wavelength (λ ≈ 2πR/ℓ). We have
3(1 + ν) ≈ 1.9 and 2π/
3(1 + ν) ≈ 3.2 if we
take ν = 1/4. This condition on λ is consistent with the transition zone between the thin and thick shell
responses analyzed in Janes and Melosh [1990] and Zhong and Zuber [2000], but does not coincide with
the constraint given in Willemann and Turcotte [1982], which is ℓ < 2π
R/h (this last condition looks
more like the threshold (83) for the membrane regime).
Though the stress distribution is affected, the limit (95) on the degree ℓ is not important for the
displacements, since they tend to zero at small wavelengths. Therefore the theory does not break down
at short wavelength if one is interested in the computation of the gravity field associated to the transverse
deflection of the lithosphere.
6 Conclusion
The principal results of this article are the five flexure equations (74)-(78) governing the three displace-
ments of the thin spherical shell and the two auxiliary stress functions. Stresses are derived quantities
which can be obtained from equations (73). The shell thickness and Young’s modulus can vary, but
Poisson’s ratio must be constant. The loads acting on the shell can be of any type since we extend the
method of stress functions to include not only transverse and consoidal tangential loads, but also toroidal
tangential loads. The flexure equations can be solved one after the other, except the two equations
(75)-(76) for the transverse deflection w and the stress function F , which must be simultaneously solved.
Tangential loading is usually neglected when solving for the deflection because of its small effect. In that
case, it is sufficient to solve the two equations (75)-(76) with Ω = H = 0:
∆′ (D∆′w)− (1− ν)A(D ;w) +R3∆′F = −R4 q ,
∆′ (α∆′F )− (1 + ν)A(α ;F )− 1
∆′w = 0 .
However tangential loading must be taken into account when computing stress fields [Banerdt , 1986].
In the long-wavelength limit (i.e. membrane regime), all equations can be solved one after the other
because it is possible to solve for F before solving for the transverse deflection. If a small part of the shell
is considered, the flexure equations reduce to the equations governing the deflection of a flat plate with
variable thickness. If the shell thickness is constant, the flexure equations reduce to equations available
in the literature which can be completely solved with spherical harmonics. Our rigorous treatment
of the thin shell approximation has clarified the effect of the shell thickness on the flexure equations.
We emphasize the need to use the correct form for the equations (without the common mistake in the
differential operator acting on w) in order to have the correct properties for the degree-one deflection and
degree-one load.
We have also obtained two general properties of the flexure equations. First we have shown that
there is always a toroidal component in the tangential displacement if the shell thickness is variable.
Second we have proven that the degree-one harmonic components of the transverse deflection and of the
toroidal component of the tangential displacement do not depend on the elastic properties of the shell.
This property reflects the freedom under translations and rotations of the reference frame. Besides we
have shown that degree-one loads are constrained by the static assumption but can deform the shell and
generate stresses.
This article was dedicated to the theoretical treatment of the flexure of a thin elastic shell with
variable thickness. While the special case of constant thickness admits an analytical solution in terms of
spherical harmonics, the general flexure equations must be solved with numerical methods such as finite
differences, finite elements or pseudospectral methods. In a forthcoming paper, we will give a practical
method of solution and discuss applications to real cases.
Acknowledgments
M. Beuthe is supported by a PRODEX grant of the Belgian Science Federal Policy. The author thanks
Tim Van Hoolst for his help and Jeanne De Jaegher for useful comments. Special thanks are due to
Patrick Wu for his constructive criticisms which helped to improve the manuscript.
7 Appendix
7.1 Covariant, contravariant and normalized components
Tensors can be defined by their transformation law under changes of coordinates. The two types of
tensor components, namely covariant and contravariant components, transform in a reciprocal way under
changes of coordinates. Tensor components cannot be expressed in a normalized basis: the space must
have a coordinate vector basis (for contravariant components) and a dual basis (for covariant components)
which are not normalized.
The only exception is a flat space with Cartesian coordinates, where covariant, contravariant and
normalized components are identical. Since the metric is the scalar product of the elements of the
coordinate vector basis, the covariant components are related to components defined in a normalized
basis (written with a hat) by
gii ûi ,
whereas the relation for contravariant components is
ûi .
The normalized Cartesian basis (x̂, ŷ, ẑ) is related to the normalized basis for spherical coordinates
(r̂, θ̂, ϕ̂) by
x̂ = cos θ cosϕ θ̂ − sinϕ ϕ̂+ sin θ cosϕ r̂ ,
ŷ = cos θ sinϕ θ̂ + cosϕ ϕ̂+ sin θ sinϕ r̂ , (96)
ẑ = − sin θ θ̂ + cos θ r̂ .
7.2 Covariant derivatives
Usual derivatives are indicated by a ‘comma’:
vi,j =
Covariant derivatives (defined below) are indicated by a ‘bar’ or by the operator ∇i:
vi|j = ∇j vi .
The former notation emphasizes the tensorial character of the covariant derivative since the covariant
derivative adds a covariant index to the vector. The latter notation is more adapted when we are interested
by the properties of the operator.
The covariant derivative of a scalar function f is equal to the usual derivative, f|i = f,i, and is itself
a covariant vector: f|i = vi. Covariant derivatives on covariant and contravariant vector components are
defined by
vi|j = vi,j − Γkij vk , (97)
vi |j = v
,j + Γ
k , (98)
where the summation on repeated indices is implicit. The symbols Γkij are the Christoffel symbols of the
second kind [Synge and Schild , 1978]. Their expressions for the metrics used in this article are given in
sections 7.3 and 7.4.
Covariant differentiation of higher order tensors is explained in Synge and Schild [1978] but we only
need the rule for a covariant tensor of second order:
σij|k = σij,k − Γlik σlj − Γljk σil .
If some of the indices of the tensor are contravariant, the rule is changed according to equation (98). The
covariant derivatives of the metric and of the inverse metric are zero: gij|k = 0 and g
7.3 Three-dimensional spherical geometry
The geometry of a thin spherical shell of average radius R can be described with coordinates θ, ϕ and
ζ, respectively representing the colatitude, longitude and radial coordinates. The radial coordinate ζ is
zero on the reference surface (i.e. the sphere of radius R) of the shell. The non-zero components of the
metric are given by
gθθ = (R + ζ)
gϕϕ = (R + ζ)
2 sin2 θ ,
gζζ = 1 .
The non-zero Christoffel symbols are given by
θθ = −(R+ ζ) ,
Γζϕϕ = −(R+ ζ) sin2 θ ,
Γθζθ = Γ
θζ = Γ
ζϕ = Γ
Γθϕϕ = − sin θ cos θ ,
ϕθ = Γ
θϕ = cot θ .
7.4 Two-dimensional spherical geometry
If θ and ϕ respectively represent the colatitude and longitude coordinates, the non-zero components of
the metric on the surface of the sphere are given by
gθθ = 1 ,
gϕϕ = sin
2 θ . (99)
The non-zero Christoffel symbols are given by
Γθϕϕ = − sin θ cos θ ,
ϕθ = Γ
θϕ = cot θ .
The double covariant derivatives of a scalar function f are thus given by
f|θ|θ = f,θ,θ ,
f|θ|ϕ = f|ϕ|θ = f,θ,ϕ − cot θ f,ϕ ,
f|ϕ|ϕ = f,ϕ,ϕ + sin θ cos θ f,θ .
An antisymmetric tensor εij is defined by
εij ≡
det gij ε̄ij ,
where ε̄ij is the antisymmetric symbol invariant under coordinate transformations: ε̄θϕ = −ε̄ϕθ = 1,
ε̄θθ = ε̄ϕϕ = 0 (ε̄ij is usually called a tensor density; Synge and Schild [1978] call it a relative tensor of
weight -1). The non-zero covariant components of εij are given for the metric of the spherical surface by
εθϕ = −εϕθ = sin θ .
The non-zero contravariant components, εij = gikgjlεkl, are given by
εθϕ = −εϕθ = csc θ .
The covariant derivative of the tensor εij is zero: εij|k = 0.
7.5 Gradient, divergence, curl and Laplacian
Various differential operators on the surface of the sphere can be constructed with covariant derivatives.
In this section, f and t are scalar functions defined on the sphere and v is a vector tangent to the sphere.
Backus [1986] gives more details on surface operators and on Helmholtz’s theorem.
As mentioned in Appendix 7.2, the covariant derivative of a scalar function f defined on the sphere
is a covariant vector tangent to the sphere whose components are f,θ and f,ϕ. The surface gradient of f
is the same vector with its components expressed in the normalized basis (θ̂, ϕ̂):
∇̄f = f,θ θ̂ + csc θ f,ϕ ϕ̂ . (100)
The contraction of the covariant derivative with the components of a vector v yields a scalar:
vi |i = v
,θ + cot θ v
θ + vϕ,ϕ .
The surface divergence is the corresponding operation on the vector with its components expressed in the
normalized basis (θ̂, ϕ̂):
∇̄ · v = csc θ
(sin θ v̂θ),θ + v̂ϕ,ϕ
. (101)
Since the result is a scalar, vi
= ∇̄ · v. A useful identity is
∇̄ · (f v) = ∇̄f · v + f ∇̄ · v . (102)
The contraction of the antisymmetric tensor εij with the covariant derivative of a scalar t yields the
covariant components of a vector v:
vi = g
jk εik t,j .
The components are given for the metric (99) by vθ = csc θ t,ϕ and vϕ = − sin θ t,θ. If t is considered
as the radial component of the radial vector t = t r̂ (the covariant radial component is equal to the
normalized one), vi are the non-zero covariant components of the three-dimensional curl of t, which is
tangent to the sphere. This fact justifies the definition of the surface curl of t, which is equal to the
vector v but with components given in the normalized basis (θ̂, ϕ̂):
∇̄ × t = csc θ t,ϕ θ̂ − t,θ ϕ̂ . (103)
The contraction of the double covariant derivative acting on a scalar f defines the surface Laplacian:
∆f = gij f|i|j
= f,θ,θ + cot θ f,θ + csc
2 θ f,ϕ,ϕ . (104)
The surface Laplacian can also be seen as the composition of the surface divergence with the surface
gradient: ∆f = ∇̄ · ∇̄f .
According to Helmholtz’s theorem, a vector tangent to the sphere can be written as the sum of the
surface gradient of a scalar f and the surface curl of a radial vector t r̂:
v = ∇̄f + ∇̄ × (t r̂) . (105)
While t is always called the toroidal scalar (or potential) for v, there is no standard terminology for f .
Backus [1986] calls f the consoidal scalar for v. Some authors [e.g. Banerdt , 1986] call f the poloidal
potential for v. The origin of this use lies in the theory of mantle convection, in which plate tectonics are
assumed to be driven by mantle flow. Under the assumption of an incompressible mantle fluid, the velocity
field of the fluid is solenoidal, i.e. its 3-dimensional divergence vanishes. In such a case, the velocity field
can be decomposed into a poloidal part (∇×∇×(P r̂)) and a toroidal part (∇×(Q r̂)), where differential
operators are 3-dimensional [Backus , 1986]. If the velocity field is tangent to the spherical surface, the
poloidal component at the surface is also the consoidal component [Forte and Peltier , 1987]. However the
fields for which we use Helmholtz’s theorem, i.e. the tangential surface load and the tangential surface
displacement, do not belong to 3-dimensional solenoidal vector fields. We thus prefer to use the term
‘consoidal’.
The surface divergence of v depends only on the consoidal scalar f :
∇̄ · v = ∆f . (106)
The two-dimensional version of Gauss theorem is
dω ∇̄ · v = 0 . (107)
where dω = sin θ dθ dϕ and the integral is taken over the whole spherical surface. It can be proven with
formula (101).
7.6 Rigid displacements
At the surface of a sphere subjected to deformation, the displacement u of a point can be expressed with
the help of Helmholtz’s theorem (105) in terms of three scalar functions (w, S, T ) depending on θ and ϕ:
u = w r̂+ ∇̄S + ∇̄ × (T r̂) . (108)
Strains (and stresses) vanish for rigid displacements. Equations (69) show that strains vanish when
(w, S, T ) are of degree one, with S = w (recall that the operatorsOi annihilate the degree one). Assuming
these conditions, we now show that utransl = w r̂ + ∇̄w represents a rigid translation whereas urot =
∇̄ × (t r̂) represents a rigid rotation of the sphere. We choose as basis the real spherical harmonics of
degree one which form the components of the radial unit vector in Cartesian coordinates:
(Yx, Yy, Yz) = (sin θ cosϕ, sin θ sinϕ, cos θ)
= (x̂, ŷ, ẑ) · r̂ . (109)
We need the surface gradient of the real spherical harmonics which can be computed with formulas (96)
and (100):
∇̄Yx, ∇̄Yy, ∇̄Yz
= (x̂− sin θ cosϕ r̂, ŷ − sin θ sinϕ r̂, ẑ− cos θ r̂) . (110)
If the expansion of w in the degree-one basis is w = a Yx + b Yy + c Yz, then
w r̂+ ∇̄w = a x̂+ b ŷ + c ẑ ,
so that utransl is indeed a rigid translation of the sphere.
If the expansion of T in the degree-one basis is T = a′Yx + b
′Yy + c
′Yz, then
∇̄ × (T r̂) = a′
− sinϕ θ̂ − cos θ cosϕ ϕ̂
cosϕ θ̂ − cos θ sinϕ ϕ̂
+ c′ sin θ ϕ̂ ,
so that urot includes a rigid rotation of the sphere, with (a
′, b′, c′) being the angles of rotation around
the axes (x̂, ŷ, ẑ), respectively. Though urot seems to include a uniform radial expansion, one should
recall that linearized strain-displacement equations are not valid for large displacements. Since strains
vanish, the radial expansion is not physical and urot represents a pure rotation. Finite deformations are
for example discussed in Love [1944][pp. 66-73] and Sokolnikoff [1956][pp. 29-33].
7.7 Differential identities for the operators Oi
The differential operators Oi defined by equations (15) satisfy differential identities useful when obtaining
the flexure equations. They are special cases of differential identities valid in curved spaces. The presence
of curvature makes the parallel transport of vectors path-dependent; this property quantifies the curvature
of space and can be expressed as the lack of commutativity of the covariant derivatives of a vector v:
vi|j|k − vi|k|j = Riljk vl , (111)
where Riljk are the covariant components of the Riemann tensor. On the sphere, the Riemann tensor
has only one independent component that is non-zero, Rθϕθϕ = − sin2 θ. Other components are related
by the symmetries Rαβγδ = −Rβαγδ = −Rαβδγ = Rγδαβ.
The substitution of f,i to vi in the commutation relation (111) provides two differential identities
satisfied by double covariant derivatives acting on scalar functions:
csc θ f|ϕ|ϕ
− csc θ f|ϕ|θ,ϕ − cos θ f|θ|θ + sin θ f,θ = 0 ,
f|θ|θ,ϕ − f|ϕ|θ,θ − cot θ f|ϕ|θ + f,ϕ = 0 .
The replacement in the above equations of the double covariant derivatives by the normalized dif-
ferential operators (31) yields the following identities:
(sin θO2f),θ − (O3f),ϕ − cos θO1f = 0 (I1) , (112)
(sin θO3f),θ − (O1f),ϕ + cos θO3f = 0 (I2) . (113)
These identities can also be directly checked with the definitions (15) of the operators Oi.
The identities (I1)-(I2) can be differentiated to generate identities of higher order. A first useful
identity is obtained from sin θ(I1),θ − (I2),ϕ = 0:
csc2 θ
sin2 θ (O2f),θ
+ (O1f),ϕ,ϕ − 2
sin θ (O3f),ϕ
− cot θ (O1f),θ + 2O1f = ∆
′f . (114)
A second useful identity is obtained from (I1),ϕ + sin θ(I2),θ = 0:
csc2 θ
sin2 θ (O3f),θ
− (O3f),ϕ,ϕ +
sin θ ((O2 −O1) f),ϕ
+ cot θ (O3f),θ − 2O3f = 0 . (115)
7.8 No degree one in operators A and B
We want to prove that the operators A and B defined by equations (33) and (37) do not have any
degree-one term in their spherical harmonic expansion:
dω A(a ; b)Y ∗1p = 0 (p = −1, 0, 1) , (116)
dω B(a ; b)Y ∗1p = 0 (p = −1, 0, 1) , (117)
where (a, b) are arbitrary scalar functions on the sphere, dω = sin θ dθ dϕ and the integral is taken over
the whole spherical surface.
This property is not a straightforward consequence of constructing A and B with Dij as a building
block. Although A and B can be factored into terms without degree one (such as Dija or ∆′a), the
product of the factors may contain degree-one terms in its spherical harmonic expansion.
Without loss of generality, we can prove the above identities with the arguments (a, b) being spherical
harmonics of given degree and order. The general result is then obtained by superposition. Let a and b
be spherical harmonics of order m and n: a ∼ eimϕ and b ∼ einϕ (we will not use their harmonic degree
in the proof). All derivatives with respect to ϕ in the operators A and B can then be replaced with the
rules a,ϕ → ima and b,ϕ → inb. The integral over ϕ in equations (116)-(117) gives
dϕ ei(m+n−p)ϕ = 2π δm+n−p,0 ,
so that the integral is zero unless n = p−m.
First consider the case p = 0 (n = −m), that is the projection on the zonal spherical harmonic of
degree one. We thus have to calculate
dθ A0 and
dθ B0 with
A0 ≡ sin θ cos θA(a ; b) ,
B0 ≡ sin θ cos θ B(a ; b) .
The trick consists in rewriting the integrands as total derivatives:
cos2 θ a,θ b,θ + cot θ
sin2 θ −m2
(ab),θ +
sin2 θ +m2 csc2 θ cos 2θ
B0 = −im
cos θ a,θ b,θ + sin θ a b,θ − csc θ cos2 θ a,θ b+ cos θ a b,θ,θ
The sought integrals are thus given by
dθ A0 =
a,θ b,θ −m2 cot θ (ab),θ +m2 csc2 θ ab
dθ B0 = −im
cos θ a,θ b,θ − csc θ a,θ b+ cos θ a b,θ,θ
where we have dropped the terms containing at least one power of sin θ which vanish at the limits; we
have also replaced cos2 θ and cos 2θ by their value at the limits. The remaining terms can be evaluated
by recalling the dependence in sin θ of the spherical harmonics: a = (sin θ)|m| a0 and b = (sin θ)
|m| b0,
where a0 and b0 are polynomials in cos θ. The only non-zero terms at the limits of the integrals are those
for |m| = 1, in which case we have at the limits: a,θb,θ = a0b0, cot θ (ab),θ = 2a0b0, csc2 θ ab = a0b0,
csc θ a,θb = cos θ a0b0, ab,θ,θ = 0. However these terms cancel in the sums so that the integrals vanish for
all m. This completes the proof for the case p = 0.
Now consider the case p = ±1 (n = −m ± 1), that is the projections on the sectoral spherical
harmonics of degree one. We thus have to calculate
dθ A±1 and
dθ B±1 with
A±1 ≡ sin2 θA(a ; b) ,
B±1 ≡ sin2 θB(a ; b) .
We again write the integrands as total derivatives:
A±1 =
sin θ cos θ a,θ b,θ −
sin θ cos θ +m2
(ab),θ −
cos2 θ ∓ 2m
a,θ b
+ sin2 θ a b,θ + 2m(m∓ 1) cot θ ab
B±1 = −i
(m∓ 1) sin θ a,θ b,θ +m sin θ a b,θ,θ − (m∓ 1) cos θ a,θ b
−m cos θ a b,θ −m(1∓m) csc θ ab
The sought integrals are thus given by
dθ A±1 =
−m2 (ab),θ − (1∓ 2m) a,θ b+ 2m(m∓ 1) cot θ ab
dθ B±1 = i
(m∓ 1) cos θ a,θ b+m cos θ a b,θ +m(1∓m) csc θ ab
where we have dropped the terms containing at least one power of sin θ and replaced cos2 θ by its value
at the limits. The remaining terms can be evaluated as in the case p = 0, but with a = (sin θ)|m| a0 and
b = (sin θ)|m∓1| b0. All terms give zero at the limits of the integrals for all values of m. This completes
the proof for the case p = ±1. We have thus proven the identities (116)-(117).
References
Anderson, F. S., and S. E. Smrekar (2006), Global mapping of crustal and lithospheric thickness on Venus,
J. Geophys. Res., 111, E08006, doi:10.1029/2004JE002395.
Arkani-Hamed, J. (1998), The lunar mascons revisited, J. Geophys. Res., 103, 3709–3739.
Arkani-Hamed, J. (2000), Strength of Martian lithosphere beneath large volcanoes, J. Geophys. Res., 105, 26,713–
26,732, doi:10.1029/2000JE001267.
Arkani-Hamed, J., and L. Riendler (2002), Stress differences in the Martian lithosphere: Constraints on the
thermal state of Mars, J. Geophys. Res., 107, 2–1, doi:10.1029/2002JE001851.
Backus, G. (1986), Poloidal and toroidal fields in geomagnetic field modeling, Rev. Geophys., 24, 75–109.
Banerdt, W. B. (1986), Support of long-wavelength loads on Venus and implications for internal structure, J. Geo-
phys. Res., 91, 403–419.
Banerdt, W. B., and M. P. Golombek (2000), Tectonics of the Tharsis Region of Mars: Insights from MGS
Topography and Gravity, in Lunar and Planetary Institute Conference Abstracts, p. 2038.
Banerdt, W. B., R. S. Saunders, R. J. Phillips, and N. H. Sleep (1982), Thick shell tectonics on one-plate planets
- Applications to Mars, J. Geophys. Res., 87, 9723–9733.
Banerdt, W. B., M. P. Golombek, and K. L. Tanaka (1992), Stress and tectonics on Mars, in Mars, edited by
Kieffer, H. H., Jakosky, B. M.,Snyder, C. W. and Matthews, M. S., pp. 249–297, University of Arizona Press,
Tucson.
Belleguic, V., P. Lognonné, and M. Wieczorek (2005), Constraints on the Martian lithosphere from gravity and
topography data, J. Geophys. Res., 110, E11005, doi:10.1029/2005JE002437.
Blewitt, G. (2003), Self-consistency in reference frames, geocenter definition, and surface loading of the solid
Earth, J. Geophys. Res., 108, 2103, doi:10.1029/2002JB002082.
Brotchie, J. F. (1971), Flexure of a liquid-filled spherical shell in a radial gravity field, Mod. Geol., 3, 15–23.
Brotchie, J. F., and R. Silvester (1969), On crustal flexure, J. Geophys. Res., 74, 5249–5252.
Comer, R. P., S. C. Solomon, and J. W. Head (1985), Mars - Thickness of the lithosphere from the tectonic
response to volcanic loads, Rev. of Geophys., 23, 61–92.
Crosby, A., and D. McKenzie (2005), Measurements of the elastic thickness under ancient lunar terrain, Icarus,
173, 100–107, doi:10.1016/j.icarus.2004.07.017.
Farrell, W. E. (1972), Deformation of the Earth by Surface Loads, Rev. Geophys. Space Phys., 10, 761–797.
Forte, A. M., and W. R. Peltier (1987), Plate tectonics and aspherical earth structure: The importance of
poloidal-toroidal coupling, J. Geophys. Res., 92, 3645–3680.
Greff-Lefftz, M., and L. Legros (1997), Some remarks about the degree-one deformation of the Earth, Geo-
phys. J. Int., 131, 699–723, doi:10.1111/j.1365-246X.1997.tb06607.x.
Hall, J. L., S. C. Solomon, and J. W. Head (1986), Elysium region, Mars - Tests of lithospheric loading models
for the formation of tectonic features, J. Geophys. Res., 91, 11,377–11,392.
Janes, D. M., and H. J. Melosh (1990), Tectonics of planetary loading - A general model and results, J. Geo-
phys. Res., 95, 21,345–21,355.
Janle, P., and D. Jannsen (1986), Isostatic gravity and elastic bending models of Olympus Mons, Mars, Annales
Geophysicae, 4, 537–546.
Janle, P., and D. Jannsen (1988), Tepev Mons on Venus: Morphology and Elastic Bending Models, Earth Moon
Planets, 41, 127–139.
Johnson, C. L., and D. T. Sandwell (1994), Lithospheric flexure on Venus, Geophys. J. Int., 119, 627–647, doi:
10.1111/j.1365-246X.1994.tb00146.x.
Johnson, C. L., S. C. Solomon, J. W. Head, R. J. Phillips, D. E. Smith, and M. T. Zuber (2000), Lithospheric
Loading by the Northern Polar Cap on Mars, Icarus, 144, 313–328, doi:10.1006/icar.1999.6310.
Kirby, J. F., and C. J. Swain (2004), Global and local isostatic coherence from the wavelet transform, Geo-
phys. Res. Lett., 31, L24608, doi:10.1029/2004GL021569.
Kraus, H. (1967), Thin Elastic Shells, John Wiley, New York.
Latychev, K., J. X. Mitrovica, M. E. Tamisiea, J. Tromp, and R. Moucha (2005), Influence of lithospheric thickness
variations on 3-D crustal velocities due to glacial isostatic adjustment, Geophys. Res. Lett., 32, L01304, doi:
10.1029/2004GL021454.
Lawrence, K. P., and R. J. Phillips (2003), Gravity/topography admittance inversion on Venus using niching
genetic algorithms, Geophys. Res. Lett., 30, 1–1, doi:10.1029/2003GL017515.
Love, A. E. H. (1944), A treatise on the mathematical theory of elasticity, 4th edition, Dover, New York.
Lowry, A. R., and S. Zhong (2003), Surface versus internal loading of the Tharsis rise, Mars, J. Geophys. Res.,
108, 5099, doi:10.1029/2003JE002111.
McGovern, P. J., S. C. Solomon, D. E. Smith, M. T. Zuber, M. Simons, M. A. Wieczorek, R. J. Phillips, G. A.
Neumann, O. Aharonson, and J. W. Head (2002), Localized gravity/topography admittance and correlation
spectra on Mars: Implications for regional and global evolution, J. Geophys. Res., 107, 5136, doi:10.1029/
2002JE001854.
McGovern, P. J., S. C. Solomon, D. E. Smith, M. T. Zuber, M. Simons, M. A. Wieczorek, R. J. Phillips, G. A.
Neumann, O. Aharonson, and J. W. Head (2004), Correction to “Localized gravity/topography admittance and
correlation spectra on Mars: Implications for regional and global evolution”, J. Geophys. Res., 109, E07007,
doi:10.1029/2004JE002286.
McKenzie, D., D. N. Barnett, and D.-N. Yuan (2002), The relationship between Martian gravity and topography,
Earth Planet. Sci. Lett., 195, 1–16.
Métivier, L., M. Greff-Lefftz, and M. Diament (2006), Mantle lateral variations and elastogravitational deforma-
tions - I. Numerical modelling, Geophys. J. Int., 167, 1060–1076, doi:10.1111/j.1365-246X.2006.03159.x.
Novozhilov, V. V. (1964), Thin shell theory, 2nd edition, Noordhoff, Groningen.
Pérez-Gussinyé, M., A. R. Lowry, A. B. Watts, and I. Velicogna (2004), On the recovery of effective elastic
thickness using spectral methods: Examples from synthetic data and from the Fennoscandian Shield, J. Geo-
phys. Res., 109, B10409, doi:10.1029/2003JB002788.
Phillips, R. J., M. T. Zuber, S. C. Solomon, M. P. Golombek, B. M. Jakosky, W. B. Banerdt, D. E. Smith,
R. M. E. Williams, B. M. Hynek, O. Aharonson, and S. A. Hauck (2001), Ancient Geodynamics and Global-
Scale Hydrology on Mars, Science, 291, 2587–2591, doi:10.1126/science.1058701.
Ranalli, G. (1987), Rheology of the Earth, Allen and Unwin, Boston.
Ronchi, C., R. Iacono, and P. S. Paolucci (1996), The cubed sphere: a new method for the solution of partial
differential equations in spherical geometry, J. Comput. Phys., 124, 93–114.
Sandwell, D. T. (1984), Thermomechanical evolution of oceanic fracture zones, J. Geophys. Res., 89, 11,401–
11,413.
Sandwell, D. T., C. L. Johnson, F. Bilotti, and J. Suppe (1997), Driving Forces for Limited Tectonics on Venus,
Icarus, 129, 232–244, doi:10.1006/icar.1997.5721.
Searls, M. L., W. B. Banerdt, and R. J. Phillips (2006), Utopia and Hellas basins, Mars: Twins separated at
birth, J. Geophys. Res., 111, E08005, doi:10.1029/2005JE002666.
Simons, M., S. C. Solomon, and B. H. Hager (1997), Localization of gravity and topography: Constraints on the
tectonics and mantle dynamics of Venus, Geophys. J. Int., 131, 24–44, doi:10.1111/j.1365-246X.1997.tb00593.x.
Sleep, N. H., and R. J. Phillips (1985), Gravity and lithospheric stress on the terrestrial planets with reference to
the Tharsis region of Mars, J. Geophys. Res., 90, 4469–4489.
Sokolnikoff, I. S. (1956), Mathematical Theory of Elasticity, McGraw-Hill, New York.
Solomon, S. C. (1978), On volcanism and thermal tectonics on one-plate planets, Geophys. Res. Lett., 5, 461–464.
Solomon, S. C., and J. W. Head (1979), Vertical movement in mare basins - Relation to mare emplacement, basin
tectonics, and lunar thermal history, J. Geophys. Res., 84, 1667–1682.
Solomon, S. C., and J. W. Head (1982), Evolution of the Tharsis province of Mars: The importance of heteroge-
neous lithospheric thickness and volcanic construction, J. Geophys. Res., 87, 9755–9774.
Solomon, S. C., and J. W. Head (1990), Heterogeneities in the thickness of the elastic lithosphere of Mars:
Constraints on heat flow and internal dynamics, J. Geophys. Res., 95, 11,073–11,083.
Stark, C. P., J. Stewart, and C. J. Ebinger (2003), Wavelet transform mapping of effective elastic thickness and
plate loading: Validation using synthetic data and application to the study of southern African tectonics,
J. Geophys. Res., 108, 2558, doi:10.1029/2001JB000609.
Stewart, J., and A. B. Watts (1997), Gravity anomalies and spatial variations of flexural rigidity at mountain
ranges, J. Geophys. Res., 102, 5327–5352, doi:10.1029/96JB03664.
Sugano, T., and K. Heki (2004), Isostasy of the Moon from high-resolution gravity and topography data: Impli-
cation for its thermal history, Geophys. Res. Lett., 31, L24703, doi:10.1029/2004GL022059.
Synge, J. L., and J. N. Schild (1978), Tensor calculus, Dover, New York.
Tanimoto, T. (1998), State of stress within a bending spherical shell and its implications for subducting lithosphere,
Geophys. J. Int., 134, 199–206, doi:10.1046/j.1365-246x.1998.00554.x.
Thurber, C. H., and M. N. Toksöz (1978), Martian lithospheric thickness from elastic flexure theory, Geo-
phys. Res. Lett., 5, 977–980.
Timoshenko, S., and S. Woinowsky-Krieger (1964), Theory of Plates and Shells (2nd edition), McGraw-Hill, New
York.
Turcotte, D. L., R. J. Willemann, W. F. Haxby, and J. Norberry (1981), Role of membrane stresses in the support
of planetary topography, J. Geophys. Res., 86, 3951–3959.
Turcotte, D. L., R. Shcherbakov, B. D. Malamud, and A. B. Kucinskas (2002), Is the Martian crust also the
Martian elastic lithosphere?, J. Geophys. Res., 107, 5091, doi:10.1029/2001JE001594.
Vlasov, V. Z. (1964), General theory of shells and its applications in engineering, NASA Tech. Trans. TT F-99,
1-886.
Wang, H., and P. Wu (2006), Effects of lateral variations in lithospheric thickness and mantle viscosity on glacially
induced surface motion on a spherical, self-gravitating Maxwell Earth, Earth Planet. Sci. Lett., 244, 576–589,
doi:10.1016/j.epsl.2006.02.026.
Wieczorek, M. A. (2007), The gravity and topography of the terrestrial planets, in Treatise on Geophysics (in
press), edited by Schubert, G., Elsevier, Amterdam.
Wieczorek, M. A., and F. J. Simons (2005), Localized spectral analysis on the sphere, Geophys. J. Int., 162,
655–675, doi:10.1111/j.1365-246X.2005.02687.x.
Willemann, R. J., and D. L. Turcotte (1982), The role of lithospheric stress in the support of the Tharsis rise,
J. Geophys. Res., 87, 9793–9801.
Zhong, S. (2002), Effects of lithosphere on the long-wavelength gravity anomalies and their implications for the
formation of the Tharsis rise on Mars, J. Geophys. Res., 107, 5054–+, doi:10.1029/2001JE001589.
Zhong, S., and J. H. Roberts (2003), On the support of the Tharsis Rise on Mars, Earth Planet. Sci. Lett., 214,
1–9, doi:10.1016/S0012-821X(03)00384-4.
Zhong, S., and M. T. Zuber (2000), Long-wavelength topographic relaxation for self-gravitating planets and
implications for the time-dependent compensation of surface topography, J. Geophys. Res., 105, 4153–4164,
doi:10.1029/1999JE001075.
Zhong, S., A. Paulson, and J. Wahr (2003), Three-dimensional finite-element modelling of Earth’s viscoelastic
deformation: effects of lateral variations in lithospheric thickness, Geophys. J. Int., 155, 679–695, doi:10.1046/
j.1365-246X.2003.02084.x.
	Introduction
	Fundamental equations of elasticity
	Three-dimensional elasticity theory
	Spherical shell
	Assumptions of the thin shell theory
	Strain-displacement equations
	Hooke's law
	Equilibrium equations
	Resolution
	Available methods
	Differential operators
	Transverse displacement
	Resolution of the equations of equilibrium
	Compatibility relation
	Tangential displacements
	Stresses
	Flexure equations and their properties
	Thin shell approximation
	Covariance
	Degree one
	Limit cases
	Membrane limit
	Euclidean limit
	Shell with constant thickness
	Displacements
	Comparison with the literature
	Breakdown of the third assumption of thin shell theory
	Conclusion
	Appendix
	Covariant, contravariant and normalized components
	Covariant derivatives
	Three-dimensional spherical geometry
	Two-dimensional spherical geometry
	Gradient, divergence, curl and Laplacian
	Rigid displacements
	Differential identities for the operators Oi
	No degree one in operators A and B
