GEOPHYSICAL RESEARCH LETTERS

Supporting Information for “Isostatic equilibrium inspherical coordinates and implications for crustal
thickness on the Moon, Mars, Enceladus, andelsewhere”
1

Douglas J. Hemingway, and Isamu Matsuyama,

2

Contents of this file
1. Text S1 to S3
2. Figures S1 to S7

Corresponding author: D. J. Hemingway, Department of Earth and Planetary Science, University of California, Berkeley, 94720, USA. (djheming@berkeley.edu)
1

Department of Earth and Planetary

Science, University of California, Berkeley,USA.
2

Lunar and Planetary Laboratory,

University of Arizona, Tucson, Arizona,
USA.

X-2

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

S1. Internal Structure, Pressure, and Gravity
S1.1. Basic interior structure
Consider a simplified model for a planetary body consisting of N concentric sphericalshells of different densities, with the requirement that the densities increase monotonicallyinward. The mean density, ρ̄, of such a body is
ρ̄ =

N
X


∆ρi

i=1

Ri
R

3
(S1)

where R is the body’s full radius, Ri is the radius of the ith layer, and ∆ρi is the density
contrast between layer i and the layer above it.
If we relax the assumption of spherical symmetry, we can express the shape of the ith
layer in terms of a spherical harmonic expansion as
Hi (θ, φ) = Ri +

∞ X
l
X

Hilm Ylm (θ, φ)(S2)

l=1 m=−lwhere θ and φ are the colatitude and longitude, respectively, Ylm (θ, φ) are the sphericalharmonic functions for degree-l and order-m [e.g., Wieczorek , 2015], Ri is the mean radiusof the ith layer, and where the coefficients Hilm describe the departure from sphericalsymmetry for the ith layer.
S1.2. Internal pressure
The hydrostatic (or lithostatic) pressure at an arbitrary location within the interior isgiven byZ ∞p (r, θ, φ) =ρ (r0 , θ, φ) g (r0 ) dr0

(S3)

r

where ρ (r0 , θ, φ) can be determined from the layer shapes, Hi (θ, φ), given the assumption of uniform density within each layer (ρi ), and where the radial component of theHEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

X-3

gravitational acceleration is given by1
g (r) = 4πG 2r

Z r

ρ (r0 ) r02 dr0

(S4)

0

where we neglect the small lateral variation in gravitational acceleration that arises because of the topography (Hilm ). For sufficiently simplified interior models, the internalhydrostatic pressure can be computed analytically by substituting (S4) into (S3) andintegrating, as we have done in order to produce Figure 1 in the main text.
It will become useful to define the mean density below radius r as
M (r)3ρ̄r = 4 3 = 3
r
πr
3

Z r

ρ (r0 ) r02 dr0

(S5)

0

allowing us to express the gravitational acceleration as
4g (r) = πGrρ̄r
3(S6)

S1.3. Gravitational potential
The gravitational potential at an arbitrary location, following the convention of thepotential being everywhere negative, is given by [e.g., Hubbard , 1984]

 Z r
Z R
1
0
0
0
0
02
0
ρ (r ) r dr + ∆U (r, θ, φ)U (r, θ, φ) = −4πG
ρ (r ) r dr +
r 0
r

(S7)

where ∆U (r, θ, φ) represents the asymmetric part of the potential, and is given by∆U (r, θ, φ) = U

rot(r, θ, φ) + U

tid(r, θ, φ) +

∞ X
l
X

Ulm (r) Ylm (θ, φ)(S8)

l=1 m=−lwhere U rot and U tid are the laterally varying rotational and (if applicable) tidal potentials,
respectively, given by
U

rot



1 1
(r, θ, φ) = ω r − + Y20 (θ, φ)
3 3
2 2

X-4

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHEREand (for a tidally locked satellite)
U

tid

2 2



(r, θ, φ) = ω r

1
1
Y20 (θ, φ) − Y22 (θ, φ)2
4

where ω is the spin period. The coefficients Ulm in (S8) account for the gravitation associated with the topography and thus depend on the layer shapes and densities (assumed
uniform within each layer), and are given by
Ulm (r) = −

N
X

4πGr
2l + 1 i=1

  Ri l+2r∆ρi Hilm  r l−1r ≥ Ri(S9)

r < RiRi

The gravitational potential exterior to the body (r ≥ R) can also be expressed as∞l

GM X XU (r, θ, φ) = −r l=0 m=−l



Rref
r

l
Clm Ylm (θ, φ)(S10)

where M is the total mass of the body, and where Clm are dimensionless gravitational
potential coefficients, with C00 = 1 by definition, and with the other coefficients being
obtained by evaluating (S9) at some reference radius, Rref , and dividing by −GM/Rref ,yielding
3Clm =2l + 1



2
Rref
ρ̄R3

X
N


∆ρi Hilm

i=1

Ri
Rref

l+2
(S11)

where again, R is the body’s mean radius, which in general may differ from the reference
radius, Rref . If Rref = R, then3Clm =2l + 1



1
ρ̄R

X
N
i=1


∆ρi Hilm

Ri
R

l+2
(S12)

S1.4. Gravitational acceleration
In spite of the lateral variations in gravitational acceleration being small enough toneglect for purposes of computing the internal pressure, we are sometimes concernedwith the magnitude of these small lateral variations, or gravity anomalies. These gravity

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

X-5

anomalies, evaluated at or above the surface (r ≥ R), can be obtained by differentiating(S9) with respect to r to get the (positive downward) gravitational acceleration
 l+2NX
l+1
Ri
glm (r) =4πG∆ρi Hilm
2l + 1
r
i=1

(S13)

When the potential is instead expressed in terms of the dimensionless potential coefficients, Clm , the equivalent result is obtained by differentiating (S10) to getGMglm (r) = (l + 1) 2
r



Rref
r

l
Clm

(S14)

S1.5. Radial variation in gravity
As we showed in the main text, the radial variation in gravity is relevant to the isostatic equilibrium problem and has implications for spectral admittance and gravity-totopography ratio (GTR) models. From (S6), it is clear that the ratio of the mean surface
gravity (gt ) to the mean gravity at the base of the crust (gb ) is
ρ̄Rt
gt
=
gb
ρ̄b Rb

(S15)

where ρ̄ is the body’s mean density and ρ̄b is the mean density of all the material below
the base of the crust. Using (S1), it can be further shown that, if the crustal density isρc , such thatρ̄ = ρc + (ρ̄b − ρc )

Rb
Rt

3
(S16)

then
(Rb /Rt )2
gt

=
gb
1 + (Rb /Rt )3 − 1 ρρ̄c

(S17)

The ratio gt /gb thus depends only on the relative compensation depth (d/Rt , whered = Rt − Rb ) and the internal density structure, the relevant part of which is effectively
captured in the ratio ρc /ρ̄ (Figure S1).

X-6

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

S2. Airy vs Pratt Isostasy
In the main text, we focused on the case of Airy isostasy, wherein the topography issupported by buoyant roots at the base of the crust (Figure S2a). In Pratt isostasy,
the topography is supported instead by lateral variations in crustal density (Figure S2b).
The problem of modeling Pratt-type isostatic compensation can be framed as the need tocompute the lateral variations in crustal density required to eliminate lateral pressure gradients at the depth of compensation (d = Rt − Rb ). Given the known surface topography(ht ) and assuming a Cartesian geometry with constant gravity, we haveρc gd = (ρc + δρ) g (d + ht )(S18)

from which it follows thatδρ = −ρcht
d + ht(S19)

Provided that |ht |  d, we can linearize (S19) to obtain the more convenient expressionδρ = −ρc

ht
d


(S20)

Whereas this Cartesian formulation is adequate in most cases, we wish to examine
the effect of taking the spherical geometry into account. Here again we examine thetwo distinct conceptions of the condition of isostatic equilibrium in spherical coordinates:
1) the requirement of equal masses in columns of equal solid angle; and 2) the requirementof the absence of lateral pressure gradients at depth, where pressure is assumed to behydrostatic. Figure S2b illustrates a cross section of a two-layer body having a crust withmean density ρc . The left hand side of the figure shows a part of the body with zero
topography with respect to the mean radius Rt . The right hand side of the figure shows a

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

X-7

part of the body with some positive topographic anomaly at the top of the crust (ht > 0)
and a corresponding crustal density anomaly, δρ.
S2.1. Equal Masses Pratt Isostasy
Equating the wedge mass in the absence of the topographic anomaly (left side of Figure S2b) with the wedge mass in the presence of the topographic anomaly (right side of
Figure S2b), yields
Z Rt

2

Z Rt +ht

r dr = (ρc + δρ)

ρc
Rb

r2 dr

Rb

After integrating and some manipulation, and making use of a substitution for thecrustal thickness (d = Rt − Rb ), we obtain 2 !

 !


1
Rb
h2t
h2t
Rb
ht
ht
δρ
+
+
+
1+
d+ 1+
ht = −ρc 1 +
ht
3
Rt
Rt
Rt 3Rt2
Rt 3Rt2
As long as |ht |  Rt , we can say that






ht



δρ ≈ −ρc 
 2


Rb
Rb
1
+
1
+
d
+
h
t
3
Rt
Rt

(S21)

As expected, this expression converges to the Cartesian equivalent, (S19), when thecompensation depth is shallow (as Rb → Rt ).
Provided that |ht |  d, we can again linearize (S21) to obtain


 
3
ht 

δρ ≈ −ρc 2 
d
b
b
1+ R
+ R
Rt
Rt

(S22)

S2.2. Equal Pressures Pratt Isostasy
Equating the hydrostatic pressure in the absence of the topographic anomaly (left sideof Figure S2b) with the hydrostatic pressure in the presence of the topographic anomaly

X-8

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE(right side of Figure S2b), in both cases evaluated at r = Rd , we obtain
Z Rt +ht

Z Rt +ht
g (r) dr = −ρcδρg (r) dr
Rt

Rb

Once again, if |ht |  Rt , then over the small radial distance between Rt and Rt +ht , theintegrand on the right hand side has a nearly constant value of gt , so the right hand side ofthis equation becomes ≈ −ρc gt ht . The integral on the left hand side is less straightforward.
However, if we assume that the gravity varies linearly from a value of gb at the base of thecrust to gt at and near the top of the crust, then we can approximate this integral by atrapezoidal region of width d + ht and starting and ending heights gb and gt , making theleft hand side of this equation ≈ δρ

gt +gb2

(d + ht ) (this approximation is good to within

2% as long as ρc /ρ̄ > 0.3 and Rb /Rt > 0.8). The result isδρ ≈ −ρcht
d + ht

2
1 + ggbt!
(S23)

As expected, this expression likewise converges to the Cartesian equivalent, (S19), when
the gravity at the compensation depth is similar to the surface gravity (i.e., as gb → gt ).
Provided that |ht |  d, we can linearize (S23) to obtain
!
 
ht
2
δρ ≈ −ρcd
1 + ggbt(S24)

Compared with the “equal pressures” condition, the “equal masses” condition always
overestimates the magnitude of the lateral density variations (Figure S3b).
S2.3. Implications under Pratt modelS2.3.1. Spectral Admittance under Pratt isostasy
Under the model of Pratt compensation, we must abandon the assumption of uniform
density within each layer, meaning that equation (S9) is no longer sufficient to describe

X-9

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

the asymmetries in the potential. If the non-hydrostatic topography is limited to thesurface (Htlm ), and if the lateral density variations (ρlm ) are confined to the outermost
layer (the crust), and do not vary radially within that layer, then the asymmetries in thegravitational potential, evaluated at or above the surface (r ≥ Rt ), are given by l+2
Z Rt +Htlm  0 l+2 !
r
4πGr
Rt
Ulm (r) = −(ρc + ρlm ) Htlm
+ ρlm
dr0
2l + 1rr
r0 =Rb
The first term represents the effect of the extra mass associated with the surface topography (Htlm ) and the second term represents the effect of the lateral density variations
within the crust (ρlm ). When the density variations are small compared with the meancrustal density (ρlm  ρc ) and when the surface relief is small compared with the meancrustal thickness (Htlm  Rt − Rb ), we can neglect the terms in (ρlm Htlm ), and thisexpression simplifies to
4πGr
Ulm (r) = −2l + 1ρc Htlm

Rt
r

l+2

Z Rt  0 l+2 !
r
+ ρlm
dr0
r
0
r =Rb

(S25)

Differentiating (S25) with respect to r, we obtain the (positive downward) gravitational
acceleration, which we evaluate at r = Rt to get
l+1Rt
glm (Rt ) =4πG ρc Htlm +2l + 1(l + 3)


1−

Rb
Rt

!

l+3 !
ρlm

Setting ρlm to the linearized expression in equation (S20), we then obtain!!
 l+3 !
l+1
1Rbρc Htlmglm (Rt ) =4πG ρc Htlm −1−b
2l + 1(l + 3)Rt
1− R
Rt
from which it follows that the spectral admittance is


Zl =



Rb
Rt

l+3 

l+1
1 1 −

4πGρc 1 −b
2l + 1(l + 3)
1− R
R
t




(S26)

X - 10HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

Using the “equal masses” equation (S22) instead yields l+3  Rb
1  1 − Rt
3
l+1



4πGρc 1 −
Zl =


  2 Rb
2l + 1(l + 3)
1 − Rt
1 + Rb + Rb
Rt

(S27)

Rt

Finally, the “equal pressures” equation (S24) leads to l+3 Rb
1  1 − Rt
l+1
4πGρc 1 −
Zl =


Rb
2l + 1(l + 3)
1 − Rt


!
2


gb
1 + gt(S28)

As with the case of Airy compensation (Figure 2), the “equal masses” model for Prattcompensation underestimates the admittance (Figure S4) compared with the “equal pressures” model. The discrepancy is especially pronounced at the lowest spherical harmonic
degrees (longest wavelengths) and when compensation depths are large. Compared with
the Airy model, admittances computed assuming a Pratt model increase more gradually,both with increasing compensation depth (Figure S4a) and increasing spherical harmonic
degree (Figure S4b). That is, for a given spherical harmonic degree, a Pratt model will
require a larger compensation depth to produce the same admittance.
S2.3.2. Geoid-to-topography ratio (GTR) under Pratt isostasy
Evaluating (S25) at r = Rt , we obtain
4πGRt
Ulm (Rt ) = −2l + 1Rt
ρc Htlm + ρlm(l + 3)


1−

Rb
Rt

l+3 !!

from which it follows that the dimensionless gravitational potential coefficients of (S10)become3Clm =2l + 1



ρc
Rt ρ̄



ρlm RtHtlm +ρc (l + 3)


1−

Rb
Rt

l+3 !!

Assuming complete Pratt compensation, with the lateral density variations computed
according to equation (S20), the geoid-to-topography ratio, given by equations (14) andX - 11HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE(15), becomes
lmax
X

GTR =



Wl

l=lmin32l + 1



ρc
ρ̄





1 −



Rb
Rt

l+3 

1 1 −
b(l + 3)
1− R
R




t

When the lateral density variations are computed according to the “equal masses”
equation (S22), we have

GTR =



lmax
X


Wl

l=lmin32l + 1



ρc
ρ̄



1 −





Rb
Rt

l+3  1 1 −
b(l + 3)
1− R
R




t


3
b
1+ R
+
Rt



Rb
Rt


 2 
(S29)

Finally, using instead the “equal pressures” equation (S24), we have
 l+3 
!
Rb
 

lX
max
1  1 − Rt
ρc 
2
3


Wl
GTR =
1 −



gb
Rb
2l + 1ρ̄(l + 3)
1 + gt
1− R
l=l

(S30)

t

min

For reference, the dipole moment approximation [Ockendon and Turcotte, 1977; Haxby
and Turcotte, 1978] gives  

3
ρc
Rb
GTR =
1−
4
ρ̄
Rt

(S31)

As with the case of Airy compensation, the “equal masses” model always underestimates
the GTR, or equivalently, overestimates the compensation depth corresponding to a given
measured GTR (Figure S5). The relative differences between the “equal masses”, “equalpressures”, and “dipole moment” models are similar under either Airy or Pratt isostasy.
The most significant difference between the Airy and Pratt models is that the Prattmodel gives much smaller (roughly half the magnitude) GTRs for any given compensation
depth. Equivalently, for a given GTR, the Pratt model requires a significantly larger
compensation depth.

X - 12HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

S3. Applications
In the main text section 3.3, we discussed several examples of how crustal thickness
estimates are affected when the “equal pressures” model for Airy isostatic equilibrium is
adopted in place of the “equal masses” model.
S3.1. The Moon and Mars
Here we reproduce Figure 3a from Wieczorek and Phillips [1997], but with one additional curve, corresponding to the “equal pressures” model described by equation (19)(Figure S6a). We likewise reproduce Figure 1 from Wieczorek and Zuber [2004], again
with the additional “equal pressures” curve (Figure S6b). In both cases, we computethe power spectra anew using more recent topography models—namely degree and order2600 spherical harmonic models for the lunar and Martian topography (both available atmarkwieczorek.github.io, and ultimately sourced from LOLA and MOLA gridded data
products available on NASA’s planetary data system). While we also repeated the GTR
computations from these and recent gravity models—namely the degree and order 420
GRAIL model for the Moon and the degree and order 120 MRO model for Mars (bothavailable on NASA’s planetary data system)—we did this only as a sanity check to ensure that our GTR values would be similar to those of Wieczorek and Phillips [1997] andWieczorek and Zuber [2004]. Rather than attempting to reproduce the masks required toisolate the highland regions of interest in both cases, we simply used the highlands-only
GTR values obtained by Wieczorek and Phillips [1997] and Wieczorek and Zuber [2004].
Whereas the GTRs for the nearside lunar highlands were given explicitly in Wieczorek

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

X - 13and Phillips [1997], we worked backwards from the reported crustal thickness estimatesof Wieczorek and Zuber [2004] to obtain the GTRs for the Martian highlands.
S3.2. Enceladus
We can also compute admittance as a function of compensation depth for the caseof Enceladus and compare it to the observed value (Figure S7). The observed degree3 admittance [Iess et al., 2014] of 14 ± 2.8 mGal/km suggests a compensation depth
of of 30 ± 6 km when the “equal masses” model is used, but just 17 ± 4 km when the“equal pressures” model is used instead. This example serves to illustrate that, when theavailable gravity/topography data are limited to the lowest spherical harmonic degrees,
the compensation depth estimate is extremely sensitive to the difference between “equalmasses” and “equal pressures “ isostasy.

References
Belleguic, V., P. Lognonné, and M. Wieczorek (2005), Constraints on the Martian lithosphere from gravity and topography data, Journal of Geophysical Research, 110 (E11),
E11,005, doi:10.1029/2005JE002437.
Beuthe, M. (2008), Thin elastic shells with variable thickness for lithospheric flexure of one-plate planets, Geophysical Journal International, 172 (2), 817–841, doi:
10.1111/j.1365-246X.2007.03671.x.
Beuthe, M., A. Rivoldini, and A. Trinh (2016), Enceladus’ and Dione’s floating ice
shells supported by minimum stress isostasy, Geophysical Research Letters, 43, doi:
10.1002/2016GL070650.

X - 14HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

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, doi:10.1002/2016GL068634.
Čadek, O., M. Běhounková, 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, Journal of Geophysical Research, 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 and Planetary Science Letters, 63 (1), 97–109, doi:
10.1016/0012-821X(83)90025-0.
Haxby, W. F., and D. L. Turcotte (1978), On Isostatic Geoid Anomalies, Journal ofGeophysical Research, 83 (B11), 5473–5478.
Heiskanen, W., and F. Vening-Meinesz (1958), Historical development of the idea ofisostasy, in The Earth and Its Gravity Field, pp. 124–146, McGraw-Hill, New York.
Hemingway, D., F. Nimmo, H. Zebker, and L. Iess (2013), A rigid and weathered ice shellon Titan, Nature, 500 (7464), 550–552, doi:10.1038/nature12400.
Hemingway, D. J., M. Zannoni, P. Tortora, F. Nimmo, and S. W. Asmar (2016), Dione’s
Internal Structure Inferred from Cassini Gravity and Topography, in Lunar and Planetary Science Conference, vol. 47, p. 1314.
Hubbard, W. B. (1984), Planetary Interiors, Van Nostrand Reinhold Company, New York.
Iess, L., N. J. Rappaport, R. A. Jacobson, P. Racioppa, D. J. Stevenson, P. Tortora, J. W.
Armstrong, and S. W. Asmar (2010), Gravity field, shape, and moment of inertia ofTitan., Science (New York, N.Y.), 327 (5971), 1367–9, doi:10.1126/science.1182583.

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

X - 15

Iess, L., et al. (2014), The Gravity Field and Interior Structure of Enceladus, Science(New York, N.Y.), 344 (6179), 78–80, doi:10.1126/science.1250551.
Jeffreys, H. (1976), The Earth: its origin, history, and physical constitution, 6th ed.,Cambridge University Press, New York.
Kraus, H. (1967), Thin Elastic Shells, Wiley, New York.
Lambeck, K. (1988), Geophysical Geodesy: The Slow Deformations of the Earth, Oxford
University Press, New York.
Lambert, W. D. (1930), The form of the geoid on the hypothesis of complete isostaticcompensation, Bulletin géodésique, 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, Journal of GeophysicalResearch, 107 (E12), 1–25, doi:10.1029/2002JE001854.
McKenzie, D., F. Nimmo, J. A. Jackson, P. B. Gans, and E. L. Miller (2000), Characteristics and consequences of flow in the lower crust, Journal of Geophysical Research,
105 (B5), 11,029–11,046, doi:10.1029/1999JB900446.
McKinnon, W. B. (2015), Effect of Enceladus’s rapid synchronous spin on interpretationof Cassini gravity, Geophysical Research Letters, 42, doi:10.1002/2015GL063384.
Nimmo, F., B. G. Bills, and P. C. Thomas (2011), Geophysical implications of the longwavelength topography of the Saturnian satellites, Journal of Geophysical Research,
116 (E11), E11,001, doi:10.1029/2011JE003835.
Ockendon, J. R., and D. L. Turcotte (1977), On the gravitational potential and fieldanomalies due to thin mass layers, Geophysical Journal of the Royal Astronomical So-

X - 16HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

ciety, 48, 479–492.
Phillips, R. J., and K. Lambeck (1980), Gravity fields of the terrestrial planets: Longwavelength anomalies and tectonics, Review of Geophysics and Space Physics, 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 flattened celestial bodies. Application
to Ceres, Astronomy & Astrophysics, 584, A127, doi:10.1051/0004-6361/201527005.
Turcotte, D. L., R. J. Willemann, W. F. Haxby, and J. Norberry (1981), Role of MembraneStresses in the Support of Planetary Topography, Journal of Geophysical Research,
86 (1), 3951–3959.
Watts, A. B. (2001), Isostasy and flexure of the lithosphere, Cambridge University Press.
Wieczorek, M. (2015), Gravity and Topography of the Terrestrial Planets, in Treatise on
Geophysics, 2nd ed., pp. 153–193, Elsevier B.V., doi:10.1016/B978-0-444-53802-4.00169X.
Wieczorek, M. A., and R. J. Phillips (1997), The structure and compensation of
the lunar highland crust, Journal of Geophysical Research, 102 (E5), 10,933, 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, Journal of Geophysical Research, 103 (E1), 1715–
1724.
Wieczorek, M. A., and M. T. Zuber (2004), Thickness of the Martian crust: Improved
constraints from geoid-to-topography ratios, Journal of Geophysical Research, 109 (E1),HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

1–16, doi:10.1029/2003JE002153.
Willemann, R. J., and D. L. Turcotte (1982), The role of lithospheric stress in the supportof the Tharsis Rise, Journal of Geophysical Research, 87 (B12), 9793–9801.
Zhong, S. (1997), Dynamics of crustal compensation and its influences on crustal isostasy,
Journal of Geophysical Research, 102.
Zhong, S. (2002), Effects of lithosphere on the long-wavelength gravity anomalies andtheir implications for the formation of the Tharsis rise on Mars, Journal of GeophysicalResearch, 107, doi:10.1029/2001JE001589.
Zhong, S., and M. T. Zuber (2000), Long-wavelength topographic relaxation for selfgravitating planets and implications for the time-dependent compensation of surface
topography, Journal of Geophysical Research, 105 (E2), 4153–4164.

X - 17

X - 18

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

1.1
1
0.9

g t /g b

0.8
0.7
0.6
0.5

0.3
0.5
0.7
0.8

0.4
0.3

0

0.1

0.2

0.3

0.4

0.5

Relative Compensation Depth (d/R)

Figure S1.

Relationship between surface gravity (gt ) and the gravity at the base of

the crust (gb ) as a function of relative compensation depth (d/Rt , where d = Rt − Rb ),for various internal density structures (described by the ratio of the crustal density to thebody’s bulk density, ρc /ρ̄).

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

a

b
ρc

ρm

ρc+δρ

ρc

ρm

Figure S2. Comparison of two hydrostatically supported wedges under the assumptionsof Airy (a) and Pratt (b) isostatic equilibrium. In both cases, the left hand side of the
figure shows a region having no topography with respect to the reference surface Rt ,while on the right, there is a topographic anomaly of height ht . In the Airy case (a), thetopography is compensated by an accompanying isostatic root at the base of the crust.
In the Pratt case (b), the nominal crustal density is ρc but the crust coinciding with the
topographic anomaly has density ρc + δρ (when ht is positive, δρ is negative).

X - 19

X - 20HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

a2.2

Basal Topography (EM/EP)

2

b

1.45
0.3
0.5
0.7
0.8

1.4

Lateral Density Anomaly (EM/EP)

0.3
0.5
0.7
0.8

1.8

1.6

1.4

1.2

1.35
1.3
1.25
1.2
1.15
1.1
1.05

1

0

0.05

0.1

0.15

Relative Compensation Depth (d/R)

Figure S3.

0.2

1

0

0.05

0.1

0.15

0.2

Relative Compensation Depth (d/R)

Comparison between “equal masses” (EM) and “equal pressures” (EP)isostasy models as a function of compensation depth and the body’s internal density
structure (characterized by the ratio ρc /ρ̄). For Airy isostasy (a), this is the factor bywhich the amplitude of the basal topography is overestimated by the “equal masses”
model. For Pratt isostasy (b), this is the factor by which the magnitude of the lateral
density variations are overestimated by the “equal masses” model. The discrepancy becomes more pronounced with increasing compensation depth and with decreasing surface
density to bulk density ratios (ρc /ρ̄).

X - 21HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

a

40
35
30
25
20
15
10

l=3
l = 10l = 30l = 100

5
0

0

0.05

0.1

0.15

Relative Compensation Depth (d/R)

Figure S4.

b

45

Admittance (mGal/km) (per c = 1000 kg/m 3)

Admittance (mGal/km) (per c = 1000 kg/m 3)

45

0.2

40
35
30
25
20
15

d/R = 0.01
d/R = 0.02
d/R = 0.05
d/R = 0.10
d/R = 0.20

10
5
0

0

10

20

30

40

50

60

70

Spherical Harmonic Degree

Admittance under Pratt isostasy. (a) Admittance as a function of relativecompensation 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 (S27), which assumes equal masses in equal columns. Solidlines show admittance as computed via (S28), which eliminates lateral pressure gradientsat depth. The “equal masses” conception of isostasy always leads to underestimating
the admittance, especially at low spherical harmonic degrees (long wavelengths). In bothpanels, 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 (S28) also depends
weakly on the internal density structure, which is here arbitrarily defined by ρc /ρ̄ = 0.6.
(cf. main text Figure 2.)

X - 22HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE
45
40
35

GTR (m/km)

30
25
20
15
10
Dipole Moment
Equal Masses
Equal Pressures

5
0

0

0.02

0.04

0.06

0.08

0.1

Relative Compensation Depth (d/R)

Figure S5.

Geoid-to-Topography Ratio (GTR) as a function of relative compensation

depth (d/R). Dotted blue line shows GTR computed via (S31), using the linear dipole
moment approximation. Dashed red line shows GTR computed via (S29), which assumes
equal masses in equal columns. Solid gold line shows GTR computed via (S30), which
avoids lateral pressure gradients at depth. The internal density structure is again arbitrarily defined by ρc /ρ̄ = 0.6. The sum in (14) is taken from lmin = 3 to lmax = 70. The
weighting coefficients are obtained from (16) by assuming a synthetic power spectrum
defined by Shh = Al−1.5 , where A is an arbitrary constant.

X - 23HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHEREaFigure S6.

bGeoid-to-Topography Ratio (GTR) as a function of crustal thickness.

Dotted blue line shows GTR computed via (20), using the linear dipole moment approximation. Dashed red line shows GTR computed via (18), which assumes equal masses inequal columns. Solid gold line shows GTR computed via (19), which avoids lateral pressure
gradients at depth. The grey bands represent the observed GTRs. Panel (a) represents
the case of the nearside lunar highlands with ρc = 2900 kg/m3 such that ρc /ρ̄ ≈ 0.87,
where the sum is taken from lmin = 3 to lmax = 70, and where the weighting coefficientsare obtained using a power spectrum derived from a degree and order 2600 LOLA-based
spherical harmonic model (cf. Figure 3a in Wieczorek and Phillips [1997]). Panel (b)represents the case of the Martian highlands with ρc = 2900 kg/m3 such that ρc /ρ̄ ≈ 0.74,
where the sum is taken from lmin = 11 to lmax = 70 (cf. Figure 1 in Wieczorek and Zuber[2004]), and where we use a degree and order 2600 MOLA-based topography model.

X - 24

Figure S7.

HEMINGWAY & MATSUYAMA: ISOSTASY ON A SPHERE

Degree-3 admittance as a function of compensation depth for Enceladus.

Dashed line shows admittance as computed via (12), which assumes equal masses in equalcolumns. Solid line shows admittance as computed via (13), which eliminates lateralpressure gradients at depth. Grey shaded region indicates degree-3 admittance obtained
from the observed gravity and topography [Iess et al., 2014; Nimmo et al., 2011].