Geophysical Research Letters
RESEARCH LETTER
10.1002/2017GL073334
Key Points:• Inconsistencies in the deﬁnitionsof isostatic equilibrium lead todiscrepancies when the sphericalgeometry of the problem isconsidered• Previous work has likely
overestimated crustal thicknesses
for small and middle-sized bodies
by as much as a factor of 2

Supporting Information:• Supporting Information S1

Correspondence to:
D. J. Hemingway,
djheming@berkeley.edu

Citation:
Hemingway, D. J., and I. Matsuyama(2017), Isostatic equilibriumin spherical coordinates andimplications for crustal thickness
on the Moon, Mars, Enceladus,and elsewhere, Geophys. Res. Lett.,
44, doi:10.1002/2017GL073334.

Received 3 MAR 2017
Accepted 4 MAY 2017
Accepted article online 16 MAY 2017

Isostatic equilibrium in spherical coordinates and implications
for crustal thickness on the Moon, Mars,
Enceladus, and elsewhere
Douglas J. Hemingway1

and Isamu Matsuyama2

1 Department of Earth and Planetary Science, University of California Berkeley, Berkeley, California, USA, 2 Lunar andPlanetary Laboratory, University of Arizona, Tucson, Arizona, USA

Abstract Isostatic equilibrium is commonly deﬁned as the state achieved when there are no lateral
gradients in hydrostatic pressure, and thus no lateral ﬂow, at depth within the lower viscosity mantle
that underlies a planetary body’s outer crust. In a constant-gravity Cartesian framework, this deﬁnition isequivalent to the requirement that columns of equal width contain equal masses. Here we show, however,that this equivalence breaks down when the spherical geometry of the problem is taken into account.
Imposing the “equal masses” requirement in a spherical geometry, as is commonly done in the literature,leads to signiﬁcant lateral pressure gradients along internal equipotential surfaces and thus corresponds toa state of disequilibrium. Compared with the “equal pressures” model we present here, the equal masses
model always overestimates the compensation depth—by ∼27% in the case of the lunar highlands and by
nearly a factor of 2 in the case of Enceladus.
Plain Language Summary “Isostasy” is the principle that, just as an iceberg ﬂoats on the water,
crustal rocks eﬀectively ﬂoat on the underlying higher density mantle, which behaves essentially like a ﬂuid
on geologic timescales. Although there are subtle inconsistencies among the various ways isostasy canbe deﬁned, they have not been historically problematic for bodies like the Earth, where the crust is thincompared with the overall radius. When the thickness of the crust is a nonnegligible fraction of a planetarybody’s radius, however, it becomes important to take the spherical geometry into account. In this case, theinconsistencies in the deﬁnitions can lead to signiﬁcant discrepancies. Here we argue that one of the mostcommonly used approaches, which requires equal width columns to contain equal masses, always results inoverestimating the crustal thickness. In particular, we suggest that the lunar and Martian highlands crustal
thickness may have been overestimated by ∼27% and ∼10%, respectively, and that the ice shell thickness
for Saturn’s small icy moon Enceladus may have been overestimated by nearly a factor of 2.

1. Introduction
Rocky and icy bodies with radii larger than roughly 200 km typically have ﬁgures that are close to the expectation for hydrostatic equilibrium (i.e., the surface conforms roughly to a gravitational equipotential) because
their interiors are weak enough that they behave like ﬂuids on geologic timescales. Because of high eﬀective
viscosities in their cold exteriors, however, these bodies can maintain some nonhydrostatic topography, evenon long timescales. This nonhydrostatic topography may be supported in part by bending and membranestresses in the lithosphere [e.g., Turcotte et al., 1981], but over long timescales, and especially when considering broad topographic loads, or loads that formed at a time when the lithosphere was weak, the rocks may
yield until much of the support comes from buoyancy—that is, the crustal material essentially ﬂoats on thehigher density, lower viscosity mantle material beneath it. This is the classic picture of isostatic equilibrium,ﬁrst discussed by Pratt and Airy in the 1850s, and is often invoked as a natural mechanism by which gravity
anomalies associated with topography can be compensated [e.g., Heiskanen and Vening-Meinesz, 1958; Watts,
2001]. The two standard end-member models for isostatic compensation are Airy, involving lateral variations
in crustal thickness, and Pratt, involving lateral variations in crustal density.
©2017. American Geophysical Union.
All Rights Reserved.

HEMINGWAY AND MATSUYAMA

The problem of modeling Airy-type isostatic compensation can be framed as the need to compute the deﬂection of the interface between the crust and the underlying higher density, lower viscosity material (we address
ISOSTASY ON A SPHERE

1Geophysical Research Letters10.1002/2017GL073334

Pratt-type compensation in the supporting information section S2). Given the known surface topography (ht ),the Airy-compensated basal topography (hb ) can be computed as
hb = −ht

𝜌c
Δ𝜌

(1)

where 𝜌c is the density of the crustal material and Δ𝜌 is the density contrast at the crust/mantle interface.
The negative sign reﬂects the fact that the basal topography is inverted with respect to the surface topography if both ht and hb are taken as positive upward relief with respect to their respective reference levels (i.e.,the hypothetical equipotential surfaces to which the density interfaces would conform if the layers were allinviscid). This equation follows from requiring equal hydrostatic pressures at equal depths (or equivalently,
requiring equal masses in columns of equal width) and ensures that, regardless of the topography, there areno horizontal pressure gradients and thus there is no lateral ﬂow at depth within the ﬂuid mantle (there isalso no vertical ﬂow because vertical pressure gradients are balanced by gravity). Hence, neglecting mantle
dynamics and the slow relaxation of the crust itself, we have a state of equilibrium.
Equation (1) implicitly assumes a Cartesian geometry and a uniform gravity ﬁeld. However, for long wavelength loads or when the compensation depth is a substantial fraction of the body’s radius, it becomes
necessary to take into account the spherical geometry of the problem. In this case, the requirement of equalmasses in equal width columns leads to (section 2.2)
( )2
Rt
𝜌
hb = −ht c
(2)
Δ𝜌 Rb
where Rt and Rb are the mean radii corresponding to the top and bottom of the crust, respectively. This
expression (or its equivalent) is widely used in the literature [e.g., Jeﬀreys, 1976; Phillips and Lambeck, 1980;
Hager, 1983; Lambeck, 1988; Wieczorek and Phillips, 1997; Wieczorek and Zuber, 2004; Hemingway et al., 2013;
McKinnon, 2015; Wieczorek, 2015; Čadek et al., 2016, 2017]. However, as we show in section 2.3, this is not
equivalent to the requirement of equal pressures at equal depths, which instead leads to( )𝜌 c gthb = −ht(3)Δ𝜌 gb
where gt and gb are the mean gravitational accelerations at the top and bottom of the crust, respectively.
Although the distinction between “equal masses” and “equal pressures” isostasy has long been recognized[e.g., Lambert, 1930; Heiskanen and Vening-Meinesz, 1958], it has widely been ignored because the eﬀect isdeemed negligible in the case of the Earth, where the crustal thickness is small compared to the radius.
However, the diﬀerence between equations (2) and (3) becomes increasingly signiﬁcant as the compensation
depth becomes an increasingly large fraction of the total radius and can therefore be important for bodies
like the Moon, Mars, Ceres, Pluto, and the outer solar system’s many middle-sized moons.
Arguably, this basic picture of isostatic equilibrium suﬀers from some internal inconsistencies in that, on the
one hand, it assumes that the crust is stiﬀ or viscous enough that the topography does not relax away completely while on the other hand assuming that the crust is weak enough that it cannot support vertical shear
stresses, meaning that radial pressure gradients are the only available means of supporting the topographic
loads against gravity. Besides handling the spherical geometry properly, a fully self-consistent conception of
the problem would have to account for the internal stresses, the elastic and rheological behaviors of the crustand mantle, the nature of the topographic loads (i.e., where and when they were emplaced), and the system’stime-varying response to those loads.
Elastic stresses may prevent or at least slow the progression toward equilibrium, especially in the case of
relatively short-wavelength loads that deﬂect, but do not readily break, the lithosphere. Accordingly, many
authors construct analytical models based on thin elastic shell theory [e.g., Kraus, 1967; Turcotte et al., 1981;
Willemann and Turcotte, 1982; McGovern et al., 2002; Hemingway et al., 2013], wherein the loads are supported
by a wavelength-dependent combination of bending stresses, membrane stresses, and buoyancy (in which
the equal masses versus equal pressures distinction remains important). Still more sophisticated approaches
exist as well. Beuthe [2008], for example, develops a more generalized analytical elastic shell model that allows
for tangential loading and laterally variable elastic properties. Taking another approach, Belleguic et al. [2005]
solve the elastic gravitational problem numerically, accommodating the spherical geometry and the force
balances in a self-consistent manner.
HEMINGWAY AND MATSUYAMAISOSTASY ON A SPHERE2

Geophysical Research Letters10.1002/2017GL073334

In the limit of a weak lithosphere (the isostatic limit), however, elastic stresses do not play such a signiﬁcantrole in supporting the topography. Some authors thus deﬁne isostatic equilibrium as the state of minimum
deviatoric stresses within the lithosphere [e.g., Dahlen, 1982; Beuthe et al., 2016]. This state is achieved insuch models by splitting the crustal thickening (or thinning) into a suitable combination of surface and basal
loads—in reality, the applied loads may have been entirely at the surface, entirely at the base, or some combination of the two; the combination that yields the state of minimum deviatoric stresses is merely intended
to represent the ﬁnal stress state after the lithosphere has ﬁnished failing or deforming in response to theapplied loads. This approach aligns well with the basic concept of complete isostatic equilibrium in that it
involves supporting the topography mainly by buoyancy but with the additional advantage of maintaining
internal consistency—deviatoric stresses do not go precisely to zero, and can thus keep the topography fromrelaxing away completely. Whereas implementation of this solution is far from straightforward [e.g., Dahlen,
1982; Beuthe et al., 2016, and references therein], our simpliﬁed approach, in spite of its limitations, leads to a
result that closely matches the minimum deviatoric stress result of Dahlen [1982].
One further consideration is the fact that relaxation may continue even after the initial gross isostatic adjustments have taken place. Provided that a topographic load is broad, and that the underlying layer is much
weaker, the system will respond relatively rapidly at ﬁrst, on a timescale governed mainly by the viscosity of
the underlying weaker mantle, until reaching a quasi-static equilibrium in which the lateral ﬂow of that weakmaterial is reduced to nearly zero. Relaxation does, however, continue after this point and may not necessarily be negligible, especially when the base of the crust is relatively warm and ductile [e.g., Zhong, 1997; Zhong
and Zuber, 2000; McKenzie et al., 2000; Zhong, 2002]. Nevertheless, this latter stage of relaxation will usually
be slow compared with the timescale for reaching isostatic equilibrium, and so we will often use the word“equilibrium” without qualiﬁcation, even as we recognize that the system may be continuing to evolve at someslow rate following the initial isostatic adjustment. We stress, however, that this is merely an assumption andthat caution should be used in cases where the materials are likely to relax more rapidly.
Notwithstanding the above complicating factors, the basic concept of isostatic equilibrium, in which topographic loads are supported entirely by buoyancy (i.e., without appeal to elastic stresses), has been widelyand productively adopted as a useful approximation in Earth and planetary sciences. To the extent that such asimpliﬁed model remains desirable for analyses of planetary topography, it should at least be consistent with
its core principle of avoiding lateral gradients in hydrostatic pressure at depth. This paper’s modest goal is toshow that, when accounting for the spherical geometry, the equal pressures model, equation (3), provides avery good approximation that is consistent with this principle, while the commonly used equal masses model,
equation (2), does not.
In section 2, we show how we obtained equations (2) and (3), and we compare the two in terms of the resulting
internal pressure anomalies. In section 3, we show how the two diﬀerent conceptions of isostasy aﬀect spectral
admittance and geoid-to-topography ratio (GTR) models, addressing implications including crustal thickness
estimates for the speciﬁc examples of the lunar and Martian highlands, as well as the ice shell thickness on
Enceladus. Finally, we make concluding remarks in section 4.

2. Analysis
2.1. Framework
Consider a body consisting of concentric layers, each having uniform density, and with the layer densities
increasing monotonically inward. The shape of the ith layer can be expanded in spherical harmonics as
Hi (𝜃, 𝜙) = Ri +

∞
l
∑
∑

Hilm Ylm (𝜃, 𝜙)(4)

l=1 m=−lwhere 𝜃 and 𝜙 are the colatitude and longitude, respectively, Ylm (𝜃, 𝜙) are the spherical harmonic functions
for degree l and order m [e.g., Wieczorek, 2015], Ri is the mean outer radius of the ith layer, and where thecoeﬃcients Hilm describe the departure from spherical symmetry for the ith layer.
Each layer’s shape is primarily a ﬁgure determined by hydrostatic equilibrium but may include smaller additional nonhydrostatic topographic anomalies. Hence, we take the shape coeﬃcients to be the sum of theirhydnhhydrostatic and nonhydrostatic parts, Hilm = Hilm + Hilm
. Since isostatic equilibrium concerns providing supnh,port for the departures from hydrostatic equilibrium, it is only the nonhydrostatic topographic anomalies, Hilm
that are involved in the isostatic equations. To a good approximation, the hydrostatic equilibrium ﬁgure can
HEMINGWAY AND MATSUYAMAISOSTASY ON A SPHERE

3Geophysical Research Letters10.1002/2017GL073334

be described by a degree-2 spherical harmonic function. Hence, this complication generally does not apply
hydto the topographic relief at degrees 3 and higher, where Hilm = 0. A possible exception is fast-rotating bodies,
for which higher-order hydrostatic terms may be nonnegligible [Rambaux et al., 2015].
We assume that the outermost shell (the “crust”) does not relax on the timescale relevant for achieving isostatic equilibrium, whereas we take the layer below the crust (the “mantle”) to be inviscid. Given the observednhnh, we are concerned with ﬁnding the basal relief, Hblm
, required to deliver
topographic relief at the surface, Htlm
isostatic equilibrium. We consider the condition of isostatic equilibrium to be satisﬁed when there are no lateral variations in hydrostatic pressure along equipotential surfaces within the inviscid layer below the crust.
The hydrostatic pressure at radial position r is given by∞p (r, 𝜃, 𝜙) =

∫r

(
) ( )
𝜌 r′ , 𝜃, 𝜙 g r′ dr′(5)where 𝜌 is the spatially varying density, where g(r) = GM(r)∕r2 is the gravitational acceleration at radius r
and where M(r) is the enclosed mass at radius r. Here the small lateral variations in gravitational acceleration
are neglected. Although lateral variations in gravity can approach a few percent due to rotation and tidal
forces, this simpliﬁcation is justiﬁed on the grounds that the quantity of interest is often the ratio gt ∕gb , as inequation (3), for example, and this ratio may be regarded as laterally constant.
A datum equipotential surface with mean radius Rd can be approximated to ﬁrst order as()
ΔU Rd , 𝜃, 𝜙Ed (𝜃, 𝜙) = Rd −
( )
g Rd

(6)

( )
where g Rd is the mean gravitational acceleration at r = Rd and where ΔU (r, 𝜃, 𝜙) represents the lateralvariations in the potential (section S1.3), given byΔU (r, 𝜃, 𝜙) = Urot (r, 𝜃, 𝜙) + Utid (r, 𝜃, 𝜙) +

l
∞
∑
∑

Ulm (r)Ylm (𝜃, 𝜙)(7)

l=1 m=−lwhere Urot and Utid are the laterally varying rotational and (if applicable) tidal potentials, respectively, andwhere the coeﬃcients Ulm account for the gravitation associated with the topography and thus depend on
the layer shapes and densities and are given by⎧ ( Ri )l+2
Nr ≥ Ri⎪4𝜋Gr ∑Ulm (r) = −Δ𝜌i Hilm ⎨ ( r )l−1
2l + 1 i=1⎪ r
r < Ri⎩ Ri(8)where Δ𝜌i is the density contrast between layer i and the layer above it.
Below, we examine two distinct conceptions of the condition of Airy-type isostasy in spherical coordinates:(1) the requirement of equal masses in columns (or cones) of equal solid angle and (2) the requirement of theabsence of lateral pressure gradients at depth, where pressure is assumed to be hydrostatic. We use simplifying assumptions to obtain compact expressions for each case. We then evaluate these simple models bycomputing lateral pressure variations along the equipotential surface deﬁned by (6). A good model should
yield little or no lateral pressure gradients along this equipotential surface. For both models, we consider atwo-layer body having a crust with density 𝜌c and an underlying mantle with density 𝜌m , where 𝜌m > 𝜌c . For
clarity and simplicity in the following derivations, we assume that the body is not subjected to rotational orhydtidal deformation so that Hilm = 0. The top and bottom of the crust have mean radii Rt and Rb , respectively.
A portion of the body has some positive topographic anomaly at the top of the crust (ht > 0) and a corresponding compensating isostatic root (inverted topography) at the base of the crust (hb < 0) (Figure S2a). A
reference datum is deﬁned at an arbitrary internal radius Rd < Rb + hb .
2.2. Equal Masses in Equal Columns
The mass above radius r, in any given column, taken as a narrow wedge, or cone, is given by∞

M=

HEMINGWAY AND MATSUYAMA∫r()
𝜌 r′ , 𝜃, 𝜙 r′2 sin 𝜃d𝜃d𝜙dr′

ISOSTASY ON A SPHERE(9)4Geophysical Research Letters10.1002/2017GL073334

where 𝜃 and 𝜙 are colatitude and longitude, respectively. Equating the wedge mass in the absence of the topographic anomaly (left side of Figure S2a) with the wedge mass in the presence of the topographic anomaly(right side of Figure S2a), yields
RbΔ𝜌

∫Rb +hb

r2 dr = 𝜌c

Rt +ht

∫ Rt

r2 dr

where Δ𝜌 = 𝜌m − 𝜌c . After integrating, and some manipulation, we obtain
)(
)−1
( ) (
h2b
h2t
𝜌 c Rt 2
ht
hb
hb = −ht
+
+
1+
1+
Δ𝜌 Rb
Rt 3R2t
Rb 3R2b
If ||ht || ≪ Rt and ||hb || ≪ Rb , this expression reduces to equation (2)
( )2
Rt
𝜌
hb ≈ −ht cΔ𝜌 Rb
2.3. Equal Pressures at Depth
Equating the hydrostatic pressure in the absence of the topographic anomaly (left side of Figure S2a) withthe hydrostatic pressure in the presence of the topographic anomaly (right side of Figure S2a), in both cases
evaluated at r = Rd , we obtain
RbΔ𝜌∫Rb +hbg(r)dr = 𝜌c

Rt +htg(r)dr∫ Rt

where again Δ𝜌 = 𝜌m − 𝜌c .
If ||ht || ≪ Rt , then over the small radial distance between Rt and Rt + ht , the integrand on the right-hand side
has a nearly constant value of gt , the mean gravitational acceleration at r = Rt . Similarly, if ||hb || ≪ Rb , then on
the left-hand side, the integrand is always close to gb , the mean gravitational acceleration at r = Rb . Hence, if
the relief at the density interfaces is small, then it is a good approximation to write
Rb

Δ𝜌gb

∫Rb +hb

dr ≈ 𝜌c gt

leading to equation (3)hb ≈ −ht𝜌cΔ𝜌

(Rt +ht

dr

∫ Rt

gt
gb

)

Because it is often more convenient to specify 𝜌c ∕𝜌̄ (the ratio of the crustal density to the body’s bulk density),
it is useful to note that gt ∕gb is given by (section S1.5)()2
Rb ∕Rt
gt
=
(10)
((
)
)3
𝜌
gb
1 + R ∕R − 1 c
b

t

𝜌̄

Note that the mass anomalies associated with the topographic anomaly and its compensating isostatic root
will displace the datum equipotential surface slightly—an eﬀect that is captured in (6) but which we haveneglected in the derivation of equation (3). If the radial displacement of this equipotential surface is hd , thehydrostatic pressure at this depth (within the mantle) will be diﬀerent by approximately 𝜌m gd hd , where gd isthe mean gravitational acceleration on this datum surface.
For comparison, Turcotte et al. [1981] include the equivalent of this additional term (which they call 𝜌m ghg ) in
their pressure balance (their equation (3)), though they neglect the radial variation in gravity and the fact thatthe shape of this equipotential surface will vary with depth (i.e., they evaluate hg only at the exterior surface,
using their equation (25)). In the limit of complete isostatic compensation, their hg goes to zero (substitutetheir equation (28) into their equation (25)). Hence, in the isostatic limit, their equation (3) is identical to ours,
except that we also account for the radial variation in gravity.
In reality, due to the ﬁnite thickness of the crust, the displacement hd will not be precisely zero (it goes to zerofor Turcotte et al. [1981] owing to some approximations they make to simplify their equation (25)), but because
we are concerned only with relatively small departures from hydrostatic equilibrium, hd is minuscule, and, as
we show in the next section, in spite of neglecting the 𝜌m gd hd term in the above derivation, equation (3) isnevertheless an excellent approximation when the goal is to make internal equipotential surfaces isobaric.
HEMINGWAY AND MATSUYAMAISOSTASY ON A SPHERE

5Geophysical Research Letters10.1002/2017GL073334

Figure 1. Comparison of internal pressure anomalies (𝛿p∕p̄ ) for various basal topography solutions. Zero pressureanomaly means zero lateral pressure gradients along the equipotential surface Ed (𝜃, 𝜙). The black dotted line illustrates
the pressure anomaly resulting from imposing the surface topography Ht30 without imposing any compensating
isostatic root (i.e., with Hb30 = 0). The colored lines illustrate the pressure anomalies obtained when the isostatic root
topography (Hb30 ) is computed via equations (1) (Cartesian isostatic balance; blue dash-dotted line), (2) (equal mass inequal columns; red dashed line), and (3) (equal pressures at depth; solid gold line).

2.4. Comparison
In spite of the simpliﬁcations used to obtain equations (2) and (3), it is clear that the two results are not equivalent. To illustrate the diﬀerence, consider the case of a two-layer body (high viscosity crust, low viscosity
mantle) that is initially spherically symmetric (for simplicity, we again assume no tidal or rotational deforming
potentials). We impose some topography at the top of the crust, Ht (𝜃, 𝜙), and compute the amplitude of thecorresponding basal topography, Hb (𝜃, 𝜙), using (1), (2), or (3). In each case, we then use (5) to compute thehydrostatic pressure at depth. Again, we are ultimately concerned with eliminating pressure gradients alongequipotential surfaces at depth, not just at a speciﬁc radial position, so we compute internal pressure along
the equipotential surface deﬁned by (6).
Figure 1 illustrates an example in which the surface topography is described by a single nonzero coeﬃcient,
Ht30 , which is longitudinally symmetric, allowing us to plot the internal pressure anomalies on an internal reference equipotential surface as a function of colatitude only. For reference, when the basal topography, Hblm ,is zero, there are of course signiﬁcant lateral variations in pressure along the equipotential surface, meaning we have a state of disequilibrium (black dotted line in Figure 1). When the topography is compensatedaccording to equation (1), the pressure anomalies are reduced but not eliminated (blue dash-dotted line).
When the topography is compensated according to equation (2), the internal pressures change substantially,
but large lateral pressure gradients remain, and so we still have a state of disequilibrium (red dashed line).
When the topography is compensated according to equation (3), on the other hand, the lateral pressure gradients nearly vanish (solid gold line), as expected if the assumptions made in section 2.3 are reasonable.
Hence, only equation (3) describes a condition that is close to equilibrium. In this example, we arbitrarily set𝜌c = 1000 kg/m3 , 𝜌m = 3000 kg/m3 , Rt = 100 km, and Rb = 80 km, such that 𝜌c ∕𝜌̄ ≈ 0.49 and Rd = 50 km, andwe impose a topographic anomaly with amplitude Ht30 = 200 m, 1% of the mean crustal thickness. The fundamental conclusions are not, however, sensitive to these choices: compared with equation (3), equation (2)always gives rise to larger pressure anomalies.
When compensation depths are shallow, gt ≈ gb and Rt ≈ Rb , so that equations (2) and (3) both reduceto the usual Cartesian form of the isostatic balance. However, when compensation depths become nonnegligible fractions of the body’s total radius, equations (1), (2), and (3) begin to diverge. When the crustal
density is less than ∼70% of the body’s bulk density, then gt < gb (section S1.5 and Figure S1), meaning that
equation (1) generally overestimates the amplitude of the basal topography. When the crustal density is morethan ∼70% of the body’s bulk density (as is likely the case for Mars, for example), gt may be larger than gb , andso equation (1) could underestimate the amplitude of the basal topography. However, of the three equations,
HEMINGWAY AND MATSUYAMAISOSTASY ON A SPHERE

6

Geophysical Research Letters10.1002/2017GL073334(2) always yields the largest (most overestimated) isostatic roots because Rt > Rb and because, assuming
density does not increase with radius, 𝜌̄ ≤ 𝜌̄b (section S1.5).

3. Implications
3.1. Spectral AdmittanceIn combined studies of gravity and topography, it is common to use the spectral admittance as a means ofcharacterizing the degree or depth of compensation [e.g., Wieczorek, 2015]. The mass associated with any
surface topography (represented using spherical harmonic expansion coeﬃcients, Htlm ) produces a corresponding gravity anomaly. However, if the topography is compensated isostatically—that is, if there is somecompensating basal topography (Hblm )—the gravity anomaly can be reduced.
Using equation (S13), we can compute the surface gravity anomaly caused by the topography at the top andbottom of the crust, yielding(( )l+2 )Rbl+1
4𝜋G 𝜌c Htlm + Δ𝜌Hblm
glm =
(11)
2l + 1
Rt
where again, 𝜌c is the density of the crust, Δ𝜌 is the density contrast at the crust/mantle interface, and where
we have neglected any contributions that may arise from asymmetries on deeper density interfaces.
Taking the degree-l admittance, Zl , to be the ratio of gravitational acceleration (glm ) to topography (Htlm ), andassuming complete Airy compensation, with the basal topography (Hblm ) computed via the equal masses
model, equation (2), we have(( )l )
Rb
l+1
4𝜋G𝜌c 1 −Zl =(12)
2l + 1Rt
Equation (12) is commonly used to generate a model admittance spectrum under the assumption of complete Airy compensation. Comparison of the model admittance with the observed admittance, along with anassumption about the crustal density then allows for an estimate of the compensation depth, d = Rt − Rb .
However, when we instead compute the basal topography using the equal pressures equation (3), we obtain(( ) ( )l+2 )
gt
Rb
l+1
4𝜋G𝜌c 1 −Zl =(13)
2l + 1gb
Rt
where again gt ∕gb is given by equation (10).
Compared with equation (13), equation (12) will always lead to an overestimate of the compensation depth.
That is, at any given spherical harmonic degree, using equation (13) yields the same admittance with a smaller
compensation depth (Figure 2a). Equivalently, for any given compensation depth, the model admittancespectrum computed via equation (13) is larger than that obtained via equation (12) (Figure 2b). The discrepancy is always greatest at low spherical harmonic degrees (e.g., focusing on degree 3, and assuming that𝜌c ∕𝜌̄ = 0.6, would yield a compensation depth estimate that is roughly ∼ 50% too large) and vanishes in the
short-wavelength limit (e.g., the compensation depth overestimate reduces to < 5% for l > 50).
For clarity and simplicity, we have not included the ﬁnite amplitude (or terrain) correction [e.g., Wieczorek
and Phillips, 1998] in the above admittance equations. When the topographic relief is a nonnegligible fraction of the body’s radius, it may be important to include this eﬀect, which will in general lead to larger
admittances. However, the point of this paper is not so much to advocate the use of equation (13) in the admittance calculation but rather, more fundamentally, to advocate the use of equation (3) in computing the basal
topography.
It is worth emphasizing that the degree-2 admittance is complicated by the eﬀects of rotational andpossibly tidal deformation. A meaningful admittance calculation for degree-2 requires ﬁrst removing thetidal/rotational eﬀects from both the gravity and topography signals. Only the remaining, nonhydrostatic,
signals should then be used in the admittance calculation. Unfortunately, determination of the hydrostatic
components of the degree-2 gravity and topography signals requires knowledge of the body’s interior structure, which may not be readily available. In such cases, the easiest option would be to simply exclude thedegree-2 terms in the admittance analysis. Alternatively, one might appeal to self-consistency arguments toconstrain the internal structure and admittance simultaneously [e.g., Iess et al., 2014].
HEMINGWAY AND MATSUYAMAISOSTASY ON A SPHERE

7

Geophysical Research Letters10.1002/2017GL073334

Figure 2. Admittance assuming Airy compensation. (a) Admittance as a function of relative compensation depth (d∕R)
for various example spherical harmonic degrees. (b) Spectral admittance for various examples of relative compensation
depths. Dashed lines show admittance as computed via (12), which assumes equal masses in equal columns. Solidlines show admittance as computed via (13), which eliminates lateral pressure gradients at depth. The equal masses
conception of isostasy always leads to underestimating the admittance, especially at low spherical harmonic degrees(long wavelengths). In both panels, admittance is normalized to an assumed crustal density of 1000 kg/m3 (i.e., if thecrustal density is 2000 kg/m3 , all admittance values double). Equation (13) also depends weakly on the internal density
structure, which is here arbitrarily deﬁned by 𝜌c ∕𝜌̄ = 0.6.

3.2. Geoid-To-Topography Ratio
A closely related concept is the geoid-to-topography ratio (GTR), which has been used to estimate regional
crustal thicknesses in situations where local isostasy can be reasonably expected [e.g., Wieczorek and Phillips,
1997; 2004]. Wieczorek and Phillips [1997] showed that the GTR is primarily a function of crustal thickness andcan be computed from a compensation model according to
GTR = Rt

lmax
∑

(14)

Wl Ql

l=lminwhere Wl is a weighting coeﬃcient for degree l and Ql is a transfer function relating the degree-l gravitational
potential and topography coeﬃcientsC
Ql = lm(15)
Hlm
The weighting coeﬃcients reﬂect the fact that the geoid is most strongly aﬀected by the longest wavelengths(lowest spherical harmonic degrees) and are constructed based on the topographic power spectrum, Shh (l) =∑l
H2 , according tom=−l lmlmax∑
Wl = Shh (l)∕Shh (i)(16)
i=lmin[Wieczorek, 2015]. Ql may be regarded as another expression for the spectral admittance (Zl ), except that it
employs dimensionless gravitational potential coeﬃcients rather than acceleration, and so we denote it herewith a distinct symbol (also in accord with Wieczorek [2015]).
Neglecting the eﬀects of topography on boundaries other than the surface and the crust/mantle interface,
we can use equation (S12) to rewrite (15) as
(
)(
( ) )
𝜌c
Δ𝜌Hblm Rb l+2
3
1+
Ql =
(17)
2l + 1 Rt 𝜌̄𝜌c Htlm
Rt
Assuming complete Airy compensation, with the basal topography (Hblm ) computed via the equal masses
equation (2), we then have

GTR =

lmax
∑
l=lminHEMINGWAY AND MATSUYAMA(
Wl

32l + 1ISOSTASY ON A SPHERE

)(𝜌 )
c

𝜌̄

(

(
1−

Rb
Rt

)l )

(18)

8

Geophysical Research Letters10.1002/2017GL073334

Figure 3. Geoid-to-Topography Ratio (GTR) as a function of relative compensation depth (d∕R). Blue dotted line shows
GTR computed via (20), using the linear dipole moment approximation. Red dashed line shows GTR computed via (18),which assumes equal masses in equal columns. Solid gold line shows GTR computed via (19), which avoids lateralpressure gradients at depth. The internal density structure is again arbitrarily deﬁned by 𝜌c ∕𝜌̄ = 0.6. The sum in (14) istaken from lmin = 3 to lmax = 70. The weighting coeﬃcients are obtained from (16) by assuming a synthetic power
spectrum deﬁned by Shh = Al−1.5 , where A is an arbitrary constant. [cf. Wieczorek and Phillips, 1997, Figure 3a; Wieczorek
and Zuber, 2004, Figure 1].

If we instead compute the basal topography using the equal pressures equation (3), we obtain(( ) ( )l+2 )
lmax
)(𝜌 )
(
∑
gt
Rb
3
c
1−
GTR =
Wl
2l
+
1
𝜌
̄
g
Rt
b
l=l

(19)

min

For reference, the linear dipole moment approximation [Ockendon and Turcotte, 1977; Haxby and Turcotte,
1978] can be written as)( )(𝜌 )(
R
3
c
1− b
GTR =
(20)
2
𝜌̄
Rt
Each model thus suggests a diﬀerent relationship between the GTR and the compensation depth (Figure 3).
For shallow compensation depths (i.e., less than ∼ 4% of the body’s radius assuming 𝜌c ∕𝜌̄ = 0.6), the equalpressures conception of isostasy and the linear dipole moment approximation give similar results. For deeper
compensation depths, the dipole moment approach begins to overestimate the GTR. In all cases, the equalmasses approach underestimates the GTR and therefore leads to an overestimate of the compensation depth(Figure 3).
3.3. Application to the Moon, Mars, and Icy Satellites
Here we consider a few realistic examples to illustrate how crustal thickness estimates diﬀer when one
adopts the equal pressures rather than the equal masses model. Note that the equal pressures-based crustal
thickness values discussed in this section should not be taken as deﬁnitive new estimates. There are many
subtleties to the interpretation of gravity and topography data that we have ignored here. The tools discussed
in sections 3.1 and 3.2 will comprise only one component of any meaningful analysis of planetary crusts.
Wieczorek and Zuber [2004], for instance, provide a comprehensive analysis that incorporates geochemical
and mechanical equilibrium considerations to complement their GTR analysis. An updated estimate of theMartian highlands crustal thickness would require careful consideration of a wide range of relevant factors
and an exploration of the permissible parameter space. Here we wish only to illustrate, using a few speciﬁc
examples, the importance of adjusting the admittance and GTR components of the analysis to incorporate
the equal pressures isostatic equilibrium model rather than the equal masses model.
For the case of the nearside lunar highlands, Wieczorek and Phillips [1997] obtained geoid-to-topography
ratios (GTRs) of roughly 14–34 m/km. Taking the case of a single layer crust (Wieczorek and Phillips [1997]
also considered dual-layer crusts), with a density of 2900 kg/m3 (𝜌c ∕𝜌̄ ≈ 0.87), this yields a crustal thickness estimate of roughly 22–61 km when the topography is assumed to be in isostatic equilibrium in the
HEMINGWAY AND MATSUYAMAISOSTASY ON A SPHERE

9Geophysical Research Letters10.1002/2017GL073334

equal masses sense. Adopting the equal pressures model instead leads to crustal thickness estimates of
18–48 km, suggesting that the equal masses model overestimates the crustal thickness by up to ∼27% in thiscase (section S3.1 and Figure S6a). For the Martian highlands, Wieczorek and Zuber [2004] obtained GTRs ofroughly 13–19 m/km, corresponding to crustal thicknesses of roughly 48–73 km, assuming a crustal densityof 2900 kg/m3 (𝜌c ∕𝜌̄ ≈ 0.74) and adopting the equal masses approach. The equal pressures model instead
leads to crustal thicknesses of roughly 44–66 km, not as dramatically diﬀerent as in the case of the lunar highlands, but still indicating that the equal masses model overestimates the crustal thickness by ∼10% in the caseof the Martian highlands (section S3.1 and Figure S6b).
For icy bodies, the ice shell’s density can be a considerably smaller fraction of the bulk density, leading tosmaller gt ∕gb ratios and therefore even more pronounced diﬀerences between the equal masses and equalpressures isostasy models (Figure S3). In the case of Europa, for example, a crustal density of 930 kg/m3 corresponds to 𝜌c ∕𝜌̄ ≈ 0.31, leading the crustal thickness estimates to diﬀer by a factor of roughly 2 at the lowest
spherical harmonic degrees. For Enceladus (𝜌c ∕𝜌̄ ≈ 0.58, assuming 𝜌c = 930 kg/m3 ), where the degree-2 anddegree-3 gravity terms have been measured based on a series of Cassini ﬂybys, Iess et al. [2014] were able toobtain a degree-3 admittance of 14.0±2.8 mGal/km, which allows for a crustal thickness estimate of 30±6 km,
adopting the equal masses model. Adopting the equal pressures model instead leads to a remarkably diﬀerent
estimate of just 17 ± 4 km (section S3.2 and Figure S7).

4. Conclusions
To the extent that isostatic equilibrium is a useful model for the state of mature planetary crusts, where broad
topographic loads are supported mainly by buoyancy, it should be taken to mean a state in which hydrostatic (or lithostatic) pressures are equal along equipotential surfaces within the relatively low viscosity mantle.
However, it is common in the literature to deﬁne isostatic equilibrium as the requirement that columnsof equal width contain equal masses. Whereas these two deﬁnitions would be equivalent in a Cartesianframework, we have shown here that they are not equivalent in a spherical geometry (section 2). We havedemonstrated that adopting the equal masses model leads to lateral pressure gradients that can be nearly aslarge (though opposite in sign) as if there were no isostatic compensation at all (Figure 1). We further showedthat the equal masses model leads to an overestimate of either the compensating basal topography in thecase of Airy compensation (section 2) or the compensating lateral crustal density variations in the case of Prattcompensation (section S2).
In combined studies of gravity and topography, using an equal masses model leads to an overestimate of thecompensation depth (Figures 2 and S4). The discrepancy is always most signiﬁcant at the lowest sphericalharmonic degrees (longest wavelengths) and increases as the crustal density becomes a smaller fraction of
the body’s bulk density. As examples, we showed that, in the case of the lunar and Martian highlands, theequal masses model could overestimate the crustal thicknesses by ∼27% and ∼10%, respectively. For the caseof Enceladus, where the compensation depth may be on the order of 10% of the radius and where the ice
shell density is roughly 58% of the bulk density, the equal masses model may overestimate the shell thickness
by nearly a factor of 2. In the case of asymmetric loads (odd harmonics), we additionally note that the equalmasses and equal pressures models will lead to distinct center of mass center of ﬁgure oﬀsets, a factor that
could be important for smaller bodies.
Whereas, for the sake of clarity, we have focused here on the end-member case of complete isostatic equilibrium (purely buoyant support), the distinction between equal masses and equal pressures remains importantfor models in which the topography is supported by a combination of both buoyancy and elastic ﬂexure—a
topic that is beyond the scope of this paper. While we acknowledge the limitations of the very concept of isostatic equilibrium (see section 1), our goal here is merely to ensure that isostasy models at least correspondto what they are intended to mean—no lateral ﬂow at depth when topographic loads are supported entirely
by buoyancy. That is, in order to be consistent with the basic principle of isostasy, we must be sure to use theequal pressures model presented here and not the equal masses model. Beyond this simple picture, a fully
self-consistent model of a planetary crust and its topography requires consideration of its loading history (i.e.,
where and when the loads were emplaced), the state of internal stresses (and failures) through time, and thepotentially time-varying rheology of the relevant materials, within both the crust and the underlying mantle.
Such models could be highly valuable but only where suﬃcient clues are available to meaningfully constrain
HEMINGWAY AND MATSUYAMAISOSTASY ON A SPHERE

10Geophysical Research Letters10.1002/2017GL073334

these many factors. In the absence of such information, the condition of isostatic equilibrium, as we havepresented it here, is likely to remain a useful model, at least as a reference end-member case.
Acknowledgments
This work was initially motivated bya discussion with Bill McKinnon andalso beneﬁted from exchanges with
Bruce Buﬀet, Anton Ermakov, Roger Fu,
Michael Manga, Tushar Mittal, Francis
Nimmo, Gabriel Tobie, and especially
Mikael Beuthe. We thank Mark
Wieczorek and Dave Stevenson for
constructive reviews that improved
the manuscript. All data are publicly
available, as described in the text.
Financial support was provided by theMiller Institute for Basic Research in
Science at the University of California
Berkeley and the NASA Gravity
Recovery and Interior Laboratory
Guest Scientist Program.

HEMINGWAY AND MATSUYAMA

References
Belleguic, V., P. Lognonné, and M. Wieczorek (2005), Constraints on the Martian lithosphere from gravity and topography data, J. Geophys.
Res., 110, E11005, doi:10.1029/2005JE002437.
Beuthe, M. (2008), Thin elastic shells with variable thickness for lithospheric ﬂexure of one-plate planets, Geophys. J. Int., 172(2), 817–841,
doi:10.1111/j.1365-246X.2007.03671.x.
Beuthe, M., A. Rivoldini, and A. Trinh (2016), Enceladus’s and Dione’s ﬂoating ice shells supported by minimum stress isostasy, Geophys. Res.
Lett., 43, 10,088–10,096, doi:10.1002/2016GL070650.
Čadek, O., et al. (2016), Enceladus’s internal ocean and ice shell constrained from Cassini gravity, shape and libration data, Geophys. Res. Lett.,
43, 5653–5660, doi:10.1002/2016GL068634.
Čadek, O., M. Běhounková, G. Tobie, and G. Choblet (2017), Viscoelastic relaxation of Enceladus’s ice shell, Icarus, 291, 31–35,
doi:10.1016/j.icarus.2017.03.011.
Dahlen, F. A. (1982), Isostatic geoid anomalies on a sphere, J. Geophys. Res., 87(B5), 3943–3947, doi:10.1029/JB087iB05p03943.
Hager, B. H. (1983), Global isostatic geoid anomalies for plate and boundary layer models of the lithosphere, Earth Planet. Sci. Lett., 63(1),97–109, doi:10.1016/0012-821X(83)90025-0.
Haxby, W. F., and D. L. Turcotte (1978), On isostatic geoid anomalies, J. Geophys. Res., 83(B11), 5473–5478.
Heiskanen, W. A., and F. A. Vening-Meinesz (1958), The Earth and Its Gravity Field, McGraw-Hill, New York.
Hemingway, D., F. Nimmo, H. Zebker, and L. Iess (2013), A rigid and weathered ice shell on Titan, Nature, 500(7464), 550–552,
doi:10.1038/nature12400.
Hubbard, W. B. (1984), Planetary Interiors, Van Nostrand Reinhold, New York.
Iess, L., et al. (2014), The gravity ﬁeld and interior structure of Enceladus, Science, 344(6179), 78–80, doi:10.1126/science.1250551.
Jeﬀreys, H. (1976), The Earth: Its Origin, History, and Physical Constitution, 6th ed., Cambridge Univ. Press, New York.
Kraus, H. (1967), Thin Elastic Shells, Wiley, New York.
Lambeck, K. (1988), Geophysical Geodesy: The Slow Deformations of the Earth, Oxford Univ. Press, New York.
Lambert, W. D. (1930), The form of the geoid on the hypothesis of complete isostatic compensation, Bull. Geod., 26(1), 98–106.
McGovern, P. J., et al. (2002), Localized gravity/topography admittance and correlation spectra on Mars: Implications for regional and global
evolution, J. Geophys. Res, 107(E12), 5136, doi:10.1029/2002JE001854.
McKenzie, D., F. Nimmo, J. A. Jackson, P. B. Gans, and E. L. Miller (2000), Characteristics and consequences of ﬂow in the lower crust,
J. Geophys. Res., 105(B5), 11,029–11,046, doi:10.1029/1999JB900446.
McKinnon, W. B. (2015), Eﬀect of Enceladus’s rapid synchronous spin on interpretation of Cassini gravity, Geophys. Res. Lett., 42, 2137–2143,
doi:10.1002/2015GL063384.
Nimmo, F., B. G. Bills, and P. C. Thomas (2011), Geophysical implications of the long-wavelength topography of the Saturnian satellites,
J. Geophys. Res., 116, E11001, doi:10.1029/2011JE003835.
Ockendon, J. R., and D. L. Turcotte (1977), On the gravitational potential and ﬁeld anomalies due to thin mass layers, Geophys. J. R. Astron.
Soc., 48, 479–492.
Phillips, R. J., and K. Lambeck (1980), Gravity ﬁelds of the terrestrial planets: Long-wavelength anomalies and tectonics, Rev. Geophys. Space
Phys., 18(1), 27–76.
Rambaux, N., F. Chambat, and J. C. Castillo-Rogez (2015), Third-order development of shape, gravity, and moment of inertia for highly
ﬂattened celestial bodies. Application to Ceres, Astron. Astrophys., 584, A127, doi:10.1051/0004-6361/201527005.
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(1), 3951–3959.
Watts, A. B. (2001), Isostasy and Flexure of the Lithosphere, Cambridge Univ. Press, Cambridge, U. K.
Wieczorek, M. A. (2015), Gravity and topography of the terrestrial planets, in Treatise on Geophysics, 2nd ed., vol. 10, pp. 153–193, Oxford,
doi:10.1016/B978-0-444-53802-4.00169-X.
Wieczorek, M. A., and R. J. Phillips (1997), The structure and compensation of the lunar highland crust, J. Geophys. Res., 102(E5),10,933–10,943, doi:10.1029/97JE00666.
Wieczorek, M. A., and R. J. Phillips (1998), Potential anomalies on a sphere: Applications to the thickness of the lunar crust, J. Geophys. Res.,
103(E1), 1715–1724.
Wieczorek, M. A., and M. T. Zuber (2004), Thickness of the Martian crust: Improved constraints from geoid-to-topography ratios, J. Geophys.
Res., 109, E01009, doi:10.1029/2003JE002153.
Willemann, R. J., and D. L. Turcotte (1982), The role of lithospheric stress in the support of the Tharsis Rise, J. Geophys. Res., 87(B12),
9793–9801.
Zhong, S. (1997), Dynamics of crustal compensation and its inﬂuences on crustal isostasy, J. Geophys. Res., 102, 15,287–15,299.
Zhong, S. (2002), Eﬀects of lithosphere on the long-wavelength gravity anomalies and their implications for the formation of the Tharsis
rise on Mars, J. Geophys. Res., 107(E7), 5054, doi:10.1029/2001JE001589.
Zhong, S., and M. T. Zuber (2000), Long-wavelength topographic relaxation for self-gravitating planets and implications for thetime-dependent compensation of surface topography, J. Geophys. Res., 105(E2), 4153–4164.

ISOSTASY ON A SPHERE11