Icarus 332 (2019) 111–131

Contents lists available at ScienceDirect

Icarusjournal homepage: www.elsevier.com/locate/icarus

Enceladus's ice shell structure as a window on internal heat production
a,b,⁎Douglas J. Hemingwayab, Tushar Mittalb

T

Department of Terrestrial Magnetism, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, D.C. 20015, USA
Department of Earth & Planetary Science, University of California Berkeley, Berkeley, CA 94720, USA

ARTICLE INFOABSTRACT

Keywords:
EnceladusInterior structure
Ice shell
Tidal heating

Its extraordinary level of geologic activity, its potential for habitability, and the prospects of returning samples
from its plume of erupting water ice make Saturn's small (∼500 km diameter) moon Enceladus a high priority
target for future exploration and a key to our developing understanding of icy ocean worlds. The structure of its
outer ice shell is particularly important as it relates to the global heat budget, the global-scale response to tidal
forces, and the nature of the ongoing eruptions. It is also diagnostic of how and where heat is dissipated internally. Here, using the most recent shape model and a new approach to modeling isostasy, we obtain a shellstructure that simultaneously accommodates the shape, gravity, and libration observations and suggests thattidal dissipation near the base of the ice shell is likely an important mode of internal heating. The implied
conductive heat loss is greater than the heat loss associated with the eruptions but is nevertheless compatible
with the condition of steady state.

1. Introduction
Enceladus's ice shell structure is a key input for understanding theongoing eruptions (Porco et al., 2006), what governs their timing, andhow ocean-to-surface pathways are initiated and sustained (Manga andWang, 2007; Rudolph and Manga, 2009; Nimmo et al., 2014; Nakajima
and Ingersoll, 2016; Ingersoll and Nakajima, 2016; Kite and Rubin,2016). The thickness and lateral variability of the ice shell also offers
clues about the evolution of Enceladus and the future of its internalliquid water ocean. Whereas accretion heat and radiogenic heating are
negligible for a body as small as Enceladus, tidal heating due to its
forced eccentric orbit can be significant (Fuller et al., 2016; Lainey
et al., 2012). How this tidal dissipation is partitioned between the outerice shell, the subsurface ocean, and the deeper interior will affect the
spatial pattern of internal heating and, in turn, the lateral shell thickness variations. Hence, the ice shell's structure also provides a window
into the mode of internal heat production. This paper thus has two
objectives: 1) to better characterize the interior of Enceladus and, inparticular, the structure of its outer ice shell; and 2) to compare thisinferred ice shell structure with the expectations from different modes
of internal heat production.
The figure of Enceladus exhibits excess polar flattening, departing
substantially from the expectation for a body that has relaxed to hydrostatic equilibrium (Nimmo et al., 2011; Tajeddine et al., 2017;
Thomas, 2010), but this non-hydrostatic topography is highly⁎

compensated, leading the gravitational field to depart only modestly
from the hydrostatic expectation (Iess et al., 2014) (Fig. 1). The relationship between the non-hydrostatic gravity and non-hydrostatic
topography can be used to estimate the compensation depth, or mean
ice shell thickness (e.g., Hemingway et al., 2013; Iess et al., 2014).
However, because such analyses focus on relatively small departuresfrom hydrostatic equilibrium, the results are sensitive to the way thehydrostatic components are computed—the usual first order approximations are not sufficient (see McKinnon, 2015); sections 2.2, 3.2.2;
Fig. 3).
To determine the shell structure, and interior structure more generally, we carried out a new analysis based on the shape, gravity, andlibration observations, incorporating a numerically accurate approach
to finding the hydrostatic equilibrium figures (section 2.2; Tricarico,2014), a new shape model and a new approach to modeling isostasy.
Crucially, unlike the early gravity-based interior models for Enceladus(Čadek et al., 2016; Hemingway et al., 2013; Iess et al., 2014;
McKinnon, 2015), we do not adopt the common approach of defining
isostasy as the requirement of maintaining equal masses in columns ofequal width (e.g., Lambeck, 1988). Instead, we adopt a new approachthat aims to eliminate lateral pressure gradients at depth (Hemingway
and Matsuyama, 2017) (section 2.3). The resulting difference in thebest fitting interior models is significant (section 3.2.2; Fig. 3) and reduces the preferred shell thickness to a range that is compatible withthe large physical libration amplitude (Fig. 4), resolving the apparent

Corresponding author.
E-mail address: dhemingway@carnegiescience.edu (D.J. Hemingway).

https://doi.org/10.1016/j.icarus.2019.03.011
Received 24 October 2018; Received in revised form 20 February 2019; Accepted 13 March 2019
Available online 26 March 20190019-1035/ © 2019 Elsevier Inc. All rights reserved.

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

4200
4000
3800

In section 3, we explore the parameter space and illustrate how it isconstrained by the shape, gravity, and libration observations, with theresults being summarized in section 3.3. In section 4, we discuss internal heat production, the heat budget, and the results of our tidal
dissipation model. Finally, concluding remarks are given in section 5.

Hydrostatic
Nimmo+2011 shape
Tajeddine+2017 shape
Iess+2014 geoid-H 2,0 (meters)

3600
3400

2. Methods

0.40

3200

The gravitational field of Enceladus was measured through radio
Doppler tracking of the Cassini spacecraft during three dedicatedgravity flybys (Iess et al., 2014). The external gravitational potential
can be expressed as a linear combination of spherical harmonic functions as follows

0.36

2800
0.34

2600
0.32

2400

0.30

2200
600

2.1. Gravitational field

0.38

3000

650

700

750

800

850

900

950

1000

l

1050

H 2,2 (meters)

Rref lClm Ylm ( , )rGMr l=0 m= l

U (r , , ) =(1)

where the potential is evaluated at radius r, colatitude θ, and longitude
ϕ, where Clm represents the degree-l and order-m dimensionless gravitational potential coefficients at some reference radius (Rref), andwhere Ylm(θ, ϕ) are the spherical harmonic functions (e.g., Wieczorek,2015). The point mass term is captured by the fact that C00 ≡ 1. In thecenter of mass frame, the degree-1 coefficients are zero. Note that
whereas Iess et al. (2014) originally reported their coefficients (Clm) fora reference radius of Rref = 254.2 km (McKinnon, 2015), we use
Rref = 252.1 km or Rref = 252.22 km corresponding to the older(Nimmo et al., 2011; Thomas, 2010) and newer (Tajeddine et al., 2017)
shape models, respectively. Note that we sometimes follow the convention of representing the zonal gravity coefficients as Jl = − Cl0.
The geoid is an equipotential surface corresponding to the meanpotential at the reference radius, Uref = − GM/Rref, and is given by

Fig. 1. Comparison of the dominant degree-2 coefficients for published shapemodels (Nimmo et al., 2011; Tajeddine et al., 2017) and the geoid computedfrom the measured gravitational field (Iess et al., 2014), with shaded ellipses
representing the 1σ uncertainties. The dashed line represents the expectation
for a relaxed hydrostatic body (section 2.2; Tricarico, 2014) with moment ofinertia factors ranging from 0.3 to 0.4 (the latter corresponding to a uniform
interior). The thin dotted line corresponds to a slope of 10/3, consistent with afirst-order approximation of the hydrostatic expectation.

discrepancy between some of the early gravity-based results(McKinnon, 2015) and the libration result (Thomas et al., 2016) (seealso Beuthe et al., 2016). Whereas, with the recent exception of Čadek
et al. (2019), all previously published Enceladus interior models are
based on the shape model of Nimmo et al. (2011), we additionally use
the newer shape model of Tajeddine et al. (2017). The newer shapemodel exhibits a smaller departure from the geoid (Fig. 1), leading to alarger preferred compensation depth, or shell thickness, though not solarge as to be incompatible with the physical libration amplitude(Fig. 4). For completeness, we undertook a comprehensive explorationof the parameter space (section 3.2), including examination of both
Airy (lateral shell thickness variations) and Pratt (lateral density variations) compensation mechanisms, and the possible role of elastic
flexure (sections 3.2.3, 3.2.4, and 3.2.5, respectively). We place ourresults in the context of previous studies in order to emphasize how
different observations (including their uncertainties) and different
modeling assumptions can affect the conclusions.
Having determined the likely shell structure and the range of possibilities that are compatible with the observations (section 3), we thendiscuss internal heat production (section 4) and compare the inferredshell structure with the outcomes expected if the internal dissipation is
primarily concentrated in the core, the ocean, or the ice shell. For thelatter case, we use a thin-shell tidal heating model (section 2.8) tocompute the spatial pattern of tidal dissipation within the ice shell(Fig. 11a) and the resulting equilibrium ice shell structure (Fig. 11c).
We show that the inferred ice shell structure (Fig. 11d) is broadly
consistent with the expectations from tidal dissipation in the ice shelland we argue that the inferred shell structure is not as readily explainable if the tidal dissipation is focussed primarily within the oceanor the core (section 4.2). Our model may thus be regarded as an alternative to the Choblet et al. (2017) model, which favors dissipation
within the core; we discuss the advantages and disadvantages of eachmodel in sections 4.2 and 4.3. We also show that maintaining the inferred shell structure near steady state requires internal heat productionin the range of 20 − 40 GW, which is marginally compatible with thedissipation implied by recent astrometric observations (Lainey et al.,2012, 2017), permitting a long-lived internal ocean (section 4.1).
Section 2 provides details on the main methods used in our analysis.

l

N ( , ) = Rref

Clm Ylm ( , )l= 0 m = l

5
6

2R 2
ref

g

Y20 ( , ) +

1
4

2R 2
ref

g

Y22 ( , )(2)

where the latter two terms account for the combined tidal and rotational potentials, where g = GM/Rref2 is the mean surface gravity. Weemphasize that the spherical harmonic functions Y20 and Y22 are notnormalized here.
We model Enceladus as a two- or three-layer body with near spherical symmetry. Departure from spherical symmetry occurs due to tidal/
rotational deformation plus the elastically/isostatically supported nonhydrostatic topography. For the ith layer, the departure from sphericalsymmetry (i.e., the spherically-referenced topography) is described by
Hilm. In general, the ith layer may also exhibit lateral variations in
density, described by δρilm.
The gravitational potential can be written

U (r , , ) =

4 G

1
r

r
0

(r ) r 2dr +Rr(r ) r drl+Ulm (r ) Ylm ( , )(3)

l=1 m= l′where ρ(r ) describes the mean (spherically symmetric) radial density
structure and where the coefficients describing the asymmetric part ofthe potential, Ulm, are given by
Ulm (r ) =4 Gr
2l + 1+i Hilmibelowr
r =0

Ri l + 2+rilm (r )

rr

i Hilm
above

l +2

dr +

Rt
r =r

r
Ri

l 1

lm (r )

+
r
r

l 1

dr

(4)

where Δρi is the density contrast between layer i and the layer above it,
and Ri is the mean outer radius of the ith layer; Rt is the mean outerradius of the top layer. The summation performed over ibelow includes
only the layers whose outer radius is below (or at) the point of interest
112

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal(i.e., r ≥ Ri) and the summation performed over iabove includes only
layers whose outer radius is above the point of interest (i.e., r < Ri).
Layers with laterally variable density will have non-zero density
anomaly expansion coefficients (δρlm) and layers with laterally variabletopography at their upper boundaries will have non-zero topography
expansion coefficients (Hilm). Eq. (4) neglects the small contributions
from terms in δρlmHilm.
Because the amplitude of the long wavelength topography can besubstantial compared with the size of Enceladus, we also include theterrain- or finite-amplitude correction when evaluating the potential
arising from each layer's topography. This correction is denoted in Eq.
(4) by the superscripts + or −, indicating that the topography coefficients have been adjusted (as described below) such that the resulting
potential accounts for the finite amplitude of the topography. This isimportant mainly for the large basal topography describing the isostaticroots in the case of our Airy compensation models. For Enceladus, andassuming ice shell and ocean densities of 925 kg/m3 and 1020 kg/m3,
respectively, the J2 part of the gravity field could be underestimated by∼5% if this finite-amplitude correction is not included. Here we follow
Wieczorek and Phillips, (1998), but with slightly modified notation. We
evaluate the topography of the ith layer, Hi(θ, ϕ), measured with respect
to the sphere Ri, and raise that topography to the nth power. The result
is then expanded in spherical harmonics, with the coefficients denoted
Hnilm. That is,inferred from our isostasy model as described in section 2.3.
2.2. Hydrostatic equilibrium figures
For each of our models, we begin by computing hydrostatic equilibrium figures for each layer. For synchronous satellites, the hydrostatic equilibrium figure is a function of the spatially varying combined
tidal and rotational potential, given by

V( , )=

Hnilm Ylm ( , ).

H hyd ( , ) =For r ≥ Ri, we then compute the finite-amplitude-adjusted topographyasn(l + 4
j=1

l +3

+
Hilm
= Ri

Hnilm
Rin
n=1

j)

n! (l + 3)

(6)

where the first term in the summation is simply the true topography
and the remaining terms constitute the small adjustments required toobtain the correct potential via Eq. (4). The identical result is obtainedby combining Eqs. (9) and (10) of Wieczorek and Phillips (1998).
When r < Ri, we instead compute the finite-amplitude-adjusted topography as

Hilm = Rin(l + j
j=1

Hnilm
Rin
n= 1

n! (l

3)

2)

(7)

which, when combined with Eq. (4) gives the same result as combining
Eqs. (11) and (12) of Wieczorek and Phillips (1998).
In practice, our analysis here always involves evaluating the potential at the reference radius (r = Rref), so Eq. (4) reduces to

Ulm (r ) =4 Gr
2l + 1+i Hilmi

Rtil+ 3 Rbil+ 3
Ri l + 2+ l +ilm
2
r
r
l+3

(8)

where Rti and Rbi are the mean radii corresponding to, respectively, thetop and bottom of the ith layer. The latter term is generally zero except
when we allow for Pratt-type isostatic compensation. Although Eq. (8)is only strictly valid at radii outside the highest standing topography(Wieczorek and Phillips, 1998), we will nevertheless use it to compute a
notional potential at the reference radius (r = Rref) for purposes ofcomparison with the measured gravitational potential, which is often
expressed at the mean outer radius rather than at some arbitrary larger
radius outside all the topography.
We treat each layer's shape (Hilm) as the sum of a hydrostatic
equilibrium figure and a non-hydrostatic component that is presumed
to be supported isostatically and/or elasticallyhydnhHilm = Hilm
+ Hilm.

1
Y22 ( , )4

(10)V( , )hfg(11)

where g is the surface gravity (radial gradient of the potential), and hf isthe fluid (or secular) Love number, which can be related to the momentof inertia—and thus the internal mass distribution—via the RadauDarwin equation (Darwin, 1899; Murray and Dermott, 1999). Thisapproach was followed in obtaining the first Enceladus interior models
based on gravity data (Iess et al., 2014). However, this approach canlead to significant errors when the results depend on small departuresfrom hydrostatic equilibrium, especially in the case of fast-rotating
bodies like Enceladus (McKinnon, 2015; Tricarico, 2014) (see section
3.2.2).
Using a fourth order recursive solution developed by Tricarico(2014) for two-layer bodies, McKinnon (2015) re-examined the analysis
of Iess et al. (2014), obtaining an updated interior structure model. The
Tricarico (2014) formulation adopts the same assumptions about the
tidal/rotational disturbing potential captured in Eq. (10), but retains
higher order terms in computing the resulting equilibrium figure, which
is approximated as a triaxial ellipsoid. We follow McKinnon (2015) inusing the Tricarico (2014) formulation except that, instead of the fourthorder recursive solution, which is limited to two-layer bodies, we use
the numerical solution, which accommodates an arbitrary number oflayers. This approach delivers exact numerical solutions for triaxial
ellipsoidal equipotential surfaces resulting from the tidal/rotational
disturbing potential of Eq. (10). We represent the resulting triaxial ellipsoidal figures as strictly degree-2 shapes with power at only m = 0
and m = 2. The reference hydrostatic equilibrium shapes for the ithhydhydlayer are thus represented by just two coefficients: Hi20
and Hi22. In
reality, these triaxial ellipsoids have some small additional power at allother even harmonics, but for purposes of making comparisons with theobserved gravity, we cannot presently make use of anything beyondl = 3, and hence a degree-2 hydrostatic figure is sufficient.
Because the non-hydrostatic topography, which is the focus of muchof this analysis, is sensitive to small changes in the hydrostatic equilibrium figure, these improvements over first order theory can be important. For example, assuming a moment of inertia of 0.335MR2 forEnceladus (Iess et al., 2014), first order theory underestimates the hyhyd
drostatic H22by ∼4% compared with the numerical solution of
Tricarico (2014), leading to a ∼30% overestimate of the non-hydronh
static component, H22. The fourth order solution yields non-hydrostatic
topography that is within 1% of the numerical results, but an advantageof the numerical approach is that it avoids having to treat the combined
ocean and ice shell as a single layer for purposes of computing hydrostatic equilibrium figures, and that it delivers equilibrium shapes for not
only the core and surface, but also for the ocean/ice interface.

(5)

l =1 m= l

5
Y20 ( , )6

where ω and R are the mean rotation rate and radius of Enceladus. Thisapproximation assumes only that the tidal potential can be well represented by a degree-2 function and that the mass of the parent body ismuch larger than the mass of the satellite (a very good assumption inthis case since Saturn is more than 106 times as massive as Enceladus).
The resulting equilibrium figure is typically computed using the firstorder approximationlHin ( , ) =

2R2

(9)

hyd
) as described in section 2.2. The
We compute the hydrostatic term (Hilmnhnon-hydrostatic terms (Hilm) may be constrained by observations or113

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

2.3. Isostatically supported non-hydrostatic topographythat simultaneously accounts for the behavior of everything from thehighly viscous ice to the dynamics of the inviscid ocean is desirable. For
the present study, however, we accept the limitations of the concept ofa static isostasy and we primarily make use of the Hemingway andMatsuyama (2017) model, though we also show results from theLambeck (1988) model as a point of comparison in section 3.2.2.
Although the Hemingway and Matsuyama (2017) approach to Airy
isostasy differs in detail from the minimum stress isostasy approach of
Beuthe et al. (2016), the two share the basic premise that the topography is supported mainly through buoyancy rather than lithosphericstresses, and neither method involves the requirement of equal masses
in columns (or cones) of equal width. Using the Nimmo et al. (2011)
shape model, our preferred ice shell and ocean thicknesses do not differ
significantly from those of Beuthe et al. (2016) (whereas they obtainshell and ocean thicknesses of 19 − 27 km and 34 − 42 km, respectively, we obtain 21 − 25 km and 37 − 44 km, respectively; Table 2).

For our three-layer models, we ensure that the non-hydrostatic
hyd
nh
obs
surface topography (Htlm= Htlm
− Htlm) is supported through a
combination of elastic forces and isostasy, either in the Airy sense oflateral shell thickness variations or in the Pratt sense of lateral density
variations.
For the end member case of Airy isostasy, the topography is supported buoyantly by isostatic roots, i.e., non-hydrostatic topography atnhthe base of the ice shell (Hblm) . For clarity, we use the subscripts t and b
to represent the top and bottom of the ice shell, respectively, rather
than numerical indices which would depend on the numbering scheme.
We set the surface relief to match the observations and compute therelief at the base of the ice shell as (Hemingway and Matsuyama, 2017)
nh
Hblm
=

nh
Htlm

ice

gt
(12)

gb

2.4. Flexural support

where ρice is the density of the ice, Δρ is the density contrast at the ice/
ocean interface, and where gt and gb are the gravitational accelerations
at the top and bottom of the ice shell, respectively.
Under the opposite end member assumption of Pratt isostasy (notexplicitly modeled in previous work on Enceladus), the high standing
non-hydrostatic topography is associated with systematically lower
crustal densities. The lateral density variations required to achieve
isostatic equilibrium, again following the approach of Hemingway andMatsuyama (2017), are given by

lm

ice

nh
Htlm
Rt Rb

Whereas Eq. (12) assumes complete Airy compensation, we also
consider the possibility of additional elastic flexural support. (We do
not consider the scenario in which flexural support is combined withpartial Pratt isostasy.) Treating the lithosphere as a uniform thin elastic
shell (Kraus, 1967; Turcotte et al., 1981; Willemann and Turcotte,1982), it can be shown (Appendix A) that Eq. (12) generalizes to
nh
Hblm
=

nh
Htlm

ice gt1
gb 1 +

2
gb

1+ g

t

(14)

gb

when the loading is entirely from the top of the shell, or

(13)

1

where ρice is the mean density of the icy crust.
Eqs. (12) and (13) account for the spherical geometry of the isostasy
problem in a way that assures approximately equal pressures alongequipotential surfaces at depth (Hemingway and Matsuyama, 2017).
Note that Eq. (12) differs from the corresponding expressions in previous studies of Enceladus's ice shell (e.g., Iess et al., 2014; McKinnon,
2015; Čadek et al., 2016) and other similar work. Many of these studies
follow the convention of defining isostasy as requiring equal masses incolumns (or cones) of equal width (e.g., Lambeck, 1988). However,Hemingway and Matsuyama (2017) showed, via integration of thehydrostatic equation, that this “equal masses” requirement does notproduce an equilibrium condition and can lead to a significantly largernh) . Conestimate of the amplitude of the required isostatic roots (Hblmsider, for example, the approach of Čadek et al. (2016), which is similar
to the Lambeck (1988) model except that they also allow for the radial
variation in gravity. In comparison to the Hemingway and Matsuyama(2017) model, Eq. (1) of Čadek et al. (2016), assuming complete
compensation, leads to a basal topography amplitude that is larger by afactor of (Rt/Rb)2 where Rt and Rb are the mean radii at the top andbottom of the ice shell, respectively. This factor can be significant in thecase of Enceladus since the mean shell thickness could be a substantial
fraction of the total radius. For example, if the mean ice shell thickness
is 25 km, then compared with our approach, Eq. (1) of Čadek et al.,2016) implies a basal topography amplitude that is larger by ∼23%.
As discussed by Hemingway and Matsuyama (2017), however, it
should be noted that the principle of isostasy as a static equilibrium
condition may be flawed in that it assumes a lithosphere that is too
weak to support any vertical shear stresses, and yet is simultaneously
strong enough to resist viscous relaxation. Ideally, the full dynamic
nature of the system should be taken into account, as has been done insome recent work (e.g., Čadek et al., 2017, 2019). These dynamic
models, however, entail a number of additional assumptions and simplifications of their own. In particular, the Čadek et al., (2017); Čadek
et al. (2019) models do not explicitly address the effects of the lateralpressure gradients within the ocean. A self-consistent dynamic model

nh
Hblm
=

nh
Htlm

ice gtgb 1 +

1(15)

ice gtwhen the loading is entirely from the base of the shell. The quantity ζ
has units of force per unit volume and serves as a shorthand for the
flexural rigidity at a particular wavelength (spherical harmonic degreel), and is given by (see Appendix A)(l, R, E , Te , ) =

ETe
R2

Te2
R212(13
32) (l (l + 1)

4l2 (l + 1)2 + 4l (l + 1)) + l (l + 1)
l (l + 1)

(1

2

)

(16)
where Te is the effective elastic thickness, and E and ν are the Young's
modulus and Poisson's ratio for the icy lithosphere. R is the radius of thethin elastic shell, which we take to be the body's outer radius. The
quantities in parentheses in Eqs. (14) and (15) may be regarded ascompensation factors. In both cases, decreasing rigidity (Te → 0, ζ → 0)
leads the compensation factor to approach unity, at which point both
Eqs. (14) and (15) reduce to Eq. (12). For non-zero rigidity, the compensation factors mean different things for the cases of bottom and top
loading. Increasing shell rigidity (increasing Te and therefore ζ) reduces
the compensation factor in either case, but this implies reduced basal
topography in the case of top loading, or reduced surface topography inthe case of bottom loading (see section 3.2.5).
2.5. Libration amplitude
Due to its eccentric orbit, Enceladus moves relatively faster near
pericenter and slower near apocenter. This results in small misalignments between the long (tidally bulged) axis of Enceladus and the line
connecting Enceladus with Saturn. This misalignment means that
Saturn exerts gravitational torques on Enceladus's tidal bulge, forcing
Enceladus to experience physical librations—small oscillations in its
otherwise steady rotation rate. The amplitude of these forced physical114

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

Table 1
Observed gravity coefficients and their uncertainties (Iess et al., 2014), along with the corresponding uncertainties from the shape models, their effect on the model
gravity uncertainties, computed via Eq. (25), and the combined effect (σlm). The gravity coefficients in the upper and lower parts of the table differ slightly because
each is expressed at a reference radius corresponding to the mean radius for each shape model: Rref = 252.1 km and Rref = 252.22 km for the Nimmo et al. (2011)
and Tajeddine et al. (2017) shape models, respectively. The assumed shell density is ρshell = 925 kg/m3, and all coefficients are unnormalized.
Shape modell, mNimmo et al. (2011)
Tajeddine et al. (2017)

2, 0
2, 2
3, 0
2, 0
2, 2
3, 0

Clmgrav

σlmgrav−6−5526 × 10
1576 × 10−6
118 × 10−6
−5521 × 10−61574 × 10−6
118 × 10−6

Hlmtopo−6

35.4 × 10
15.9 × 10−6
23.4 × 10−6
35.4 × 10−6
15.8 × 10−6
23.4 × 10−6

librations is a function of the orbital eccentricity and moments of inertia
of Enceladus. In the case of a multi-layer body, with all the layers
physically coupled to one another, the physical libration amplitude isgiven by (e.g., Murray and Dermott, 1999)

=

1

2e
(n / 0 ) 2

3(B

(17)

A)

L=

(18)

C

A
C

=

i
i

2
bi2 ) ai bi cii (a i22
i (ai + bi ) ai bi ci

(19)
th

where Δρi is the density contrast between the i layer and the layer
above it, and where ai, bi and ci are, respectively, the semi-major, semiintermediate, and semi-minor axes corresponding to the ith layer, which
is taken to be a triaxial ellipsoid.
In the special case of a three layer body in which the intermediate
layer is taken to be a global liquid layer that physically decouples thecore from the outer ice shell, the physical libration amplitude is instead
given by (e.g., Richard et al., 2014; Thomas et al., 2016)=

2e [K s (K c + 2Kint n2Cc) + 2Kint K c ]2
2
2
Cs Cc (n21 )(n2)

178.9 m19.4 m
4.8 m
3.9 m
1.3 m
4.5 m

σlm

∗

−6

244.8 × 10
26.5 × 10−6
4.7 × 10−6
5.3 × 10−6
1.7 × 10−6
4.4 × 10−6

247.3 × 10−6
30.9 × 10−6
23.9 × 10−6
35.8 × 10−6
15.9 × 10−6
23.8 × 10−6

(X

Y)T

1 (X

(21)

Y)

where X is a vector containing the model results (gravitational potential
coefficients and libration amplitude) and where Y is a vector containing
the observations (the observed gravitational potential coefficients (Iesset al., 2014) and libration amplitude (Thomas et al., 2016)), and where
it is understood that only the measurably non-zero gravitational potential coefficients are included (i.e., C20, C22, C30, in the present case).
The covariance matrix, Σ, captures the uncertainties in both the modeland the observations, and is given by

where A < B < C are the principle moments of inertia. These are captured in the moon's dynamical triaxiality, given byB−3846 m
917 m
384 m
−3510 m
857 m
420 m

σlmtopodimensional parameter space exhaustively (with steps of 1 km in layer
thicknesses, 25 kg/m3 in layer densities, and 10 m in elastic thicknesses) in order to find the best fitting parameters. We use the techniques described above to develop interior models over the full grid, andthen compute the resulting gravitational field using Eq. (8) and the libration amplitude using Eqs. (17) or (20), depending on whether adecoupling ocean is present. At each point in the parameter space, wethen quantify the goodness of fit using the Mahalanobis distance, given
bywhere e is Enceladus's orbital eccentricity, n is its mean motion, and ω0
is its natural frequency, given by
0 = n

σlmtopo

=

obs +

(22)

model

obs

where Σ is the covariance matrix corresponding to the observations,consisting of the estimated gravitational potential coefficients covariances (D. Durante & L. Iess, personal communication) and the varianceof the estimated libration amplitude (Thomas et al., 2016); the libration
amplitude and gravitational potential coefficients are assumed to be
uncorrelated. Likewise, Σmodel is a covariance matrix comprising model
uncertainties. Uncertainties in the model gravitational potential arise
from the propagation of uncertainties from the shape models (Nimmo
et al., 2011; Tajeddine et al., 2017) through Eq. (8). Similarly, uncertainties in the model libration amplitude arise from propagating theshape model uncertainties through the libration calculation describedin section 2.5. Because it is impractical to directly propagate these
uncertainties every time the libration amplitude is computed from themodel, we instead used a Monte-Carlo approach to determine how thecomputed libration amplitude uncertainty is related to the shape model
uncertainties, and then used a conservative fixed value for each shapemodel. The Nimmo et al. (2011) shape model uncertainties (seeTable 1) translate to a 1σ uncertainty of <0.004∘ in the computed libration amplitude while the Tajeddine et al. (2017) uncertainties
translate to a 1σ uncertainty of <0.00025∘. Squaring this value gives
the last term in Σmodel. Since the shape model coefficients were assumed
to be independent (i.e., covariances were not published by either
Nimmo et al. (2011) or Tajeddine et al. (2017)), we must accept theapproximation that Σmodel is diagonal.
Assuming the errors are normally distributed, we compute theprobability density at every point in the parameter space as(20)

where Cs and Cc are the separate polar moments of inertia for the shelland core, respectively. In this case, there are multiple separate torques
to consider: Ks and Kc are the amplitudes of the effective torques applied by Saturn on the shell and core, respectively, and are mainly
functions of moment of inertia differences (B-A) in those layers.
Meanwhile, Kint is the amplitude of the effective torque applied by thecore on the shell, and is a function of the shapes and density contrasts atthe core/ocean boundary and at the ocean/shell boundary (Richard
et al., 2014). Finally, ω1 and ω2 are the system's natural frequencies. For
further details, see Richard et al. (2014); Thomas et al. (2016); Van
Hoolst et al. (2009, 2013). Whereas the core is assumed to be hydrostatic (though see 3.2.6), the outer shape of the ice shell comes directlyfrom the observations (Nimmo et al., 2011; Tajeddine et al., 2017) andthe shape of the ice/ocean interface is obtained as described in section
2.3 or 2.4. Uncertainties in the outer shape model are handled as described in section 2.6 below. Consistent with Thomas et al. (2016), we
do not consider the effects of diurnal tidal variations in the figures,
which have been shown to reduce the libration amplitude, at least in
the case of large ocean worlds (Van Hoolst et al., 2013); the effect isprobably small for Enceladus (Van Hoolst et al., 2016).

f (X , Y , ) =2.6. Goodness of fit

1| |(2 ) k2

e L /2(23)

where k is the number of degrees of freedom, which we take to be thenumber of measurable gravitational potential coefficients plus the

For each of our models, we explore the two-, three-, or four115

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

60

10 9

0
200

9 8

400

600

800

Shell Thickness (km)

60

10
9
180

40

0
200

9
7

0
200

400

600

800

1000

Shell Density (kg/m3)

600

800

1000

Shell Density (kg/m3)

2400
2300
2200

60
40

2000
170

180

190

200

210

220

230

Core Radius (km)

h

2700

20

400

2500

2800

g

80

0
200

2600

2100

100

80

20

40

120

f

100

Shell Thickness (km)1000

Shell Density (kg/m3)

0
18
89

60

20

120

e

7 8
190

80

d

2700

Core Density (kg/m 3)

20

10
7
710
98

2800

c

100

Shell Thickness (km)8010Shell Thickness (km)

100

40

120

b

Core Density (kg/m 3)

120

a

2600
2500
2400
2300
2200
2100

400

600

800

1000

Shell Density (kg/m3)

2000
170

180

190

200

210

220

230

Core Radius (km)Fig. 2. Misfits and probability densities across the parameter space assuming a two-layer model, and comparing results derived from the shape models of (a–d)
Nimmo et al. (2011) and (e–h) Tajeddine et al. (2017). (a,e) Lines indicate combinations of shell thickness and density that satisfy the observed J2 (blue), C22 (red), orJ3 (gold) dimensionless gravitational potential coefficients (Iess et al., 2014), with shaded bands indicating 1σ (dark) and 2σ (pale) bounds from combined shape andgravity model uncertainties (section 2.6). (b,f) Misfit between the model and the observations, given by Eq. (21). (c,g) 68% (dark), 95% (intermediate), and 99.7%(pale) confidence contours showing shell thicknesses and densities that best account for the observed gravitational potential, and (d,h) corresponding core radius anddensity (section 2.6). The minimum misfits are L = 6.1 and L = 6.7 for the Nimmo et al. (2011) and Tajeddine et al. (2017) shape models, respectively, but require
unreasonably small ice shell densities. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

a

b

c

d

e

f

g

h

Fig. 3. Sensitivity to methodological choices, adopting the shape models of (a–d) Nimmo et al. (2011) or (e–h) Tajeddine et al. (2017). The curves illustrate the partsof the parameter space that satisfy the observed J2 (blue), C22 (red), or J3 (gold) dimensionless gravitational potential coefficients (Iess et al., 2014), and the observed
libration amplitude (purple) (Thomas et al., 2016), with shaded bands indicating 1σ (dark) and 2σ (pale) uncertainties (section 2.6). (a,e) Results obtained following
the approach of Iess et al. (2014). (b,f) Results obtained following the approach of McKinnon (2015), wherein the equilibrium figure is computed following theapproach of Tricarico (2014). (c,g) Results obtained with our three-layer model, which includes the finite amplitude correction (Wieczorek and Phillips, 1998), andretaining the “equal masses” conception of isostasy used in earlier studies. (d,h) Results obtained when the “equal masses” isostatic model is replaced with the “equalpressures” model of Hemingway and Matsuyama (2017). For all panels, we assume ρshell = 925 kg/m3 and ρocean = 1020 kg/m3. The minimum misfits when using
the Nimmo et al. (2011) shape model are: (a) L = 0.4, (b) L = 1.5, (c) L = 3.3, (d) L = 1.6. When using the Tajeddine et al. (2017) shape model, they are: (e) L = 1.3,(f) L = 2.9, (g) L = 6.0, (h) L = 3.7. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
116

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

9

60
8
7

9

7

6

8
7

5

8
2

8

3

6

4
81095 6
7

10

20

30

40

7

8109

50

60

c

50
40
30
20

0

70

f

0

10

20

30

40

50

60

70

Ocean Thickness (km)

g

60

30

4
5

6
18
0

10
0

1980
7
5

97

20

0

10

20

30

40

50
40
30
20
10

50

60

70

0

2400
2300
2200

2000
170

180

190

200

210

220

230

Core Radius (km)

h

2700

Core Density (kg/m 3)

Shell Thickness (km)

6

80
1

Shell Thickness (km)1
7 89 0

9

40

2500

2800

10

50

2600

2100

70

60

d

2700

10

Ocean Thickness (km)

70

2800

60

6
5

0

3

4

10 9 7

20
10

e

6

5

9

4

30

0

10

9
10

Shell Thickness (km)

50
40

b

10

8

70

Core Density (kg/m 3)a9Shell Thickness (km)

70

2600
2500
2400
2300
2200
2100

0

10

20

Ocean Thickness (km)

30

40

50

60

70

2000
170

180

Ocean Thickness (km)

190

200

210

220

230

Core Radius (km)Fig. 4. Misfits and probability densities across the parameter space assuming Airy compensation, with ρshell = 925 kg/m3 and ρocean = 1020 kg/m3, and comparing
results derived from the shape models of (a–d) Nimmo et al. (2011) and (e–h) Tajeddine et al. (2017). (a,e) Lines indicate combinations of shell and ocean thicknesses
that satisfy the observed J2 (blue), C22 (red), or J3 (gold) dimensionless gravitational potential coefficients (Iess et al., 2014), and the observed libration amplitude(purple) (Thomas et al., 2016), with shaded bands indicating 1σ (dark) and 2σ (pale) uncertainties on the estimated values (section 2.6); same as Fig. 3d, h. (b,f)
Misfit between the model and the observations, given by Eq. (21). (c,g) 68% (dark), 95% (intermediate), and 99.7% (pale) confidence contours showing ice shell andocean thicknesses that best account for the observed gravitational potential and libration amplitude, and (d,h) corresponding core radius and density. The sharp
lower boundary on mean shell thickness in (c) is the result of excluding parts of the parameter space where the south polar shell thickness is less than zero; the same
criterion is applied in (g) but the effect is subtle because the south polar topographic depression is less pronounced in the Tajeddine et al. (2017) shape model. The
minimum misfits are L = 1.6 and L = 3.7 for the Nimmo et al. (2011) and Tajeddine et al. (2017) shape models, respectively. (For interpretation of the references tocolor in this figure legend, the reader is referred to the web version of this article.)libration amplitude (γ). In the present case, k = 4 (i.e., C20, C22, C30, γ).
We do not attempt to evaluate the likelihood of each model in absolute terms, which would require some assumptions about priorprobabilities across the parameter space. Instead, the contours in
Figs. 2, 4, 6, S1, S2, S3, and S4, are chosen such that they enclose
68.27% (dark), 95.45% (intermediate), and 99.73% (pale) of theprobability density-weighted volume of the explored parameter space,with the implicit assumption that the cumulative probability within thisvolume is equal to unity. Hence, our confidence contours show only the
relative likelihood of the model parameters given each particular set ofmodel assumptions (i.e., the chosen shape and gravity models, thechosen isostasy model, the ranges of layer densities considered, and therange of allowed elastic support). For our three-layer Airy models,
however, we require that the ice shell and ocean thickness be everywhere greater than zero. This is accomplished by forcing the probability
density to zero in parts of the parameter space where this condition isviolated.
The shaded bands in Figs. 2, 3, 4, 6, and 8 illustrate the combined 1σ
or 2σ uncertainties, neglecting any correlation between the terms. Thatis, the width of each dark band is the standard deviation, given by thesquare root of the diagonal elements of the combined covariance matrix
Σ. Equivalently, the uncertainties relating to the gravitational potential
coefficients are given bylm =grav 2
topo 2( lm) + ( lm
)

topo
lm

=

3 shell
R (2l + 1)

topo
lm .

(25)

Similarly, for the libration amplitude, the shaded purple bands in
Figs. 3 and 4 are the combined uncertainties, given by=2,obs +

2,model(26)

where σγ, obs = 0.014∘ (Thomas et al., 2016) and where σγ, model arises
from the shape model uncertainties described above, with values ofσγ, model = 0.004∘ and σγ, model = 0.00025∘ for the shape models of
Nimmo et al. (2011) and Tajeddine et al. (2017), respectively.
2.7. Radiative equilibrium surface temperature
To compute the equilibrium surface temperature of Enceladus as afunction of latitude, we take an approach similar to that of Ojakangas
and Stevenson (1989) for Europa. Some differences include the fact that
we cannot assume a small angle for the inclination (solar obliquity)
because, whereas the Jupiter system is inclined by just ∼3.13° withrespect to its orbital plane, the Saturn system's inclination is ∼26.73°.
Also, rather than dividing the problem into its Arctic and non-Arctic
parts, making appropriate approximations, and solving analytically, we
simply carry out the full calculation numerically.
Consider a point, P, on the surface of Enceladus at colatitude θ. In an
inertial frame centered on Enceladus with the z-axis coinciding with thespin vector, we have a unit vector pointed towards P and another unit
vector pointed towards the sun. The angle, ψ, between the two vectors
is given by

(24)

where
are the uncertainties on the measured coefficients (Iesstopo
et al., 2014) and where lm
are the uncertainties in the model potential field arising from the shape model uncertainties (Nimmo et al.,
topo
2011; Tajeddine et al., 2017), lm. The two are related, neglecting theeffects of compensation, via Eqs. (8) and (1), as follows
grav
lm

117

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittalcos= sin ·cos i·cos

t ·cos nt + sin ·sin

+ sin i·cos ·cos

t·sin nt

deq ( , ) =

(27)

t

where i is the inclination of the Saturn system with respect to its orbital
plane (and thus the solar obliquity of Enceladus given that its obliquity
with respect to Saturn is not measurably different from zero), where Ωis the mean motion of Saturn's orbit around the sun, and where n is themean motion of Enceladus in its orbit around Saturn.
The time- and latitude-dependent solar flux, F(θ, t), at a point on thesurface of Enceladus is the product of the solar constant at Saturn(Fs ≈ 14.9 W/m2) and cos ψ, and is obviously zero when cos ψ ≤ 0,
corresponding to the night side. That is

Fs cosF ( , t) =
0

cos
cos

Tb
T=0

2 /

2

t=0

Tb

(29)

Finally, we compute the equilibrium surface temperature asT( ) =(1A) FA ( )

1/4

(30)
−8

−2 −4

where A = 0.8 is the surface albedo and σ = 5.67 × 10 Wm K is
the Stefan-Boltzmann constant. The approximation proposed by
Nadeau and Mcgehee (2017) delivers a nearly identical result.

=d c dT.
dz T dz(33)

q (T ) dT
2µ
=
( , )
T
nl
2

tan 1 n 0µ

(34)

T( , )cln sd( , )Tb

1d( , )
R(35)

where the term in parentheses is required to account for the sphericalgeometry. Although this equation is only strictly valid when there are
no heat sources within the ice, this approximation is adequate because
the dissipation is concentrated very close to the base of the shell, suchthat the heat flux is constant across nearly all of the ice shell.
The total tidal dissipation within Enceladus can be related to theproperties of its orbit and its interior via (Peale and Cassen, 1978)

The diurnally variable tidal stresses cause strains that dissipate heatwithin the ice shell. The volumetric heat dissipation, q, is a function of
the assumed ice rheology and varies spatially across the body as well aswith depth (due to the dependence on viscosity and thus temperature).
In an equilibrium configuration, the thickness and temperature structure of the ice shell will be such that all of the heat dissipated at a given
latitude and longitude will be removed radially outward by conduction(due to the much smaller temperature gradients, lateral heat flow will
be negligible by comparison). That is, the surface heat flux will equal
the total heat dissipated within the ice shell at a given latitude andlongitude (plus any additional heat that is supplied to the shell frombelow). Given some surface temperature (Ts; section 2.7) and basal
temperature (Tb), we wish to solve for the equilibrium shell thickness,deq.
Assuming radial outward heat conduction and using a coordinate
system with z positive downward, the volumetric heat dissipation at
any given point within the ice shell is also the radial gradient of thelocal heat flux, and can thus be related to the temperature structure as
follows

dT
k
dz

(32)

q (T ) dTT

Fcond ( , ) =2.8. Tidal dissipation, equilibrium shell thickness, and heat lossdq (T ) =
dz

1/2
H 2
c

where μ is the rigidity of the ice, n is the mean motion of the satellite inits orbit (the same as its rate of spin), and η0 is a reference viscosity
corresponding to the base of the ice shell. The term ( , ) is the timeaveraged value of the sum of the squares of the elements of the strain
rate tensor (see Ojakangas and Stevenson, 1989, Appendix B, noting thecorrections by Nimmo et al. (2007)). Note that this approach assumes
the dissipation takes place in a uniform thin shell. If the low viscosity
part of the shell has significant variability in thickness, dissipation may
be enhanced in the thinner regions (e.g., Běhounková et al., 2017;
Beuthe, 2018).
The laterally varying conductive heat loss may be approximated by(28)

F ( , t ) dt .

( ))is a function of the assumed ice rheology. Assuming a Maxwell
rheology, for example, yields (Ojakangas and Stevenson, 1989)

After integrating (numerically, with dt = 1 hr) and averaging overan entire Saturnian year (2π/Ω), we obtain the latitude-dependent
annual mean flux

FA ( ) =Tb q (T ) dT
+
T =0
T

where θ and ϕ are the colatitude and longitude, respectively, and where
the dissipation integral

T=0

>0
.
0

(

ln (Tb/Ts ( , ))2cE=

21 k2 GMS2 5 2RE ne
2 Q a6(36)

where k2 is the degree-2 tidal potential Love number for Enceladus,
describing the magnitude of the response to tidal forcing, Q is thequality factor related to how much energy is dissipated on each tidal
cycle, MS is the mass of Saturn, RE is the radius of Enceladus, a is its
semi-major axis, n is its mean motion, and e is its orbital eccentricity.
See Nimmo et al. (2018) for a more thorough discussion of tidal dissipation as it relates to Enceladus.
3. Ice shell structure
3.1. Observations
The shape of Enceladus has been estimated based on analysis oflimb profile measurements (Nimmo et al., 2011; Thomas, 2010) and,more recently, by additionally incorporating data from a control point
network (Tajeddine et al., 2017). The latter approach differs in detailfrom the earlier work, in particular in the way uncertainties are calculated. Because the newer analysis incorporates additional data, one
would expect somewhat smaller uncertainties. However, the degree-2
coefficient uncertainties reported by Tajeddine et al. (2017) are smaller
than those of Nimmo et al. (2011) by a factor of 16–46, depending on
which coefficients are compared. We regard this as an indication thatthe uncertainties reported by Nimmo et al. (2011) and Tajeddine et al.
(2017) could have somewhat distinct physical meanings, and may nottherefore be directly comparable. We proceed here using the as-published uncertainties in both shape models, but we urge caution incomparing the results, as we discuss further below.
In spite of the differences between the published shape models, it is(31)

Here, we have expressed the thermal conductivity as a function oftemperature, k = c/T, with c = 651 W/m being an experimentally derived constant (Petrenko and Whitworth, 1999, p.43). Following theapproach of Ojakangas and Stevenson (1989), we can integrate thisheat dissipation across the shell subject to the boundary conditions thatthe surface temperature is given by Eq. (30) and that the heat flux into
the shell from below is H. Assuming that dissipation is a strong functionof viscosity and thus temperature, such that nearly all of the dissipationis concentrated near the relatively low viscosity base of the shell, andthat the gradient in ln T is approximately constant across almost theentire shell, then the laterally variable equilibrium shell thickness isgiven by (Ojakangas and Stevenson, 1989)
118

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

clear that Enceladus's shape departs substantially from the expectation
for a body in hydrostatic equilibrium (Fig. 1). Whereas, for a hydrostatic Enceladus, the ratio between the two major degree-2 terms
should be −H20/H22 ≈ 3.25 (e.g., McKinnon, 2015; Tricarico, 2014),
the observed ratio is 4.20 ± 0.22 or 4.09 ± 0.01, according to theshape models of Nimmo et al. (2011) and Tajeddine et al. (2017), respectively. The corresponding ratio for the geoid, computed using Eq.
(2) from the measured gravitational field (Iess et al., 2014), is just
3.42 ± 0.02. While this still represents a statistically significant departure from hydrostatic equilibrium, it is much closer to the hydrostatic expectation (dashed line in Fig. 1), immediately suggesting thatthe non-hydrostatic topography is largely compensated.
Because of the departure from hydrostatic equilibrium, we cannot
obtain the moment of inertia directly. Instead, we must treat Enceladus
as a mostly hydrostatic body (corresponding to some point on the hydrostatic line in Fig. 1) superimposed with some non-hydrostatic topography and its corresponding non-hydrostatic gravity. Finding a solution that satisfies the observational constraints, while handlingcompensation self-consistently, requires some exploration of parameter
space.
To illustrate the way in which the different observations constrain
the interior structure, we use shaded bands (in Figs. 2, 3, 4, 6, and 8) toshow the parts of the parameter space that satisfy each observation
within their uncertainties (see section 2.6). Table 1 shows, for the important gravity coefficients (C20, C22, C30), the relative contributions to
the overall uncertainty from the separate uncertainties in the gravity
and shape models.

range, computing the resulting gravity coefficients for each model using
Eq. (8). Fig. 2 illustrates the ranges of model parameters that yieldgravitational fields that match the observed J2 (blue), C22 (red), or J3(gold) dimensionless gravitational potential coefficients (Iess et al.,2014). We compute the libration amplitude, via Eq. (17), across thisparameter space but its maximum value is ∼0.04°, more than 5σ below
the observed value of 0.120 ± 0.014°. Only models with ice shell densities below ∼600 kg/m3 can simultaneously satisfy the observed J2
and C22 within their 1σ bounds, though there is considerable agreement
within their 2σ bounds. To satisfy the observed J3, however, requires
unreasonably small ice shell densities. The magnitude of the J3 gravity
term, obtained via Eq. (8), and resulting mainly from the large topographic depression at the south pole, is generally far too large (e.g., forρshell = 900 kg/m3, J3model ≈ − 375 × 10−6, nearly 10σ from the observed value). Hence, for reasonable ice shell densities, and independent of any other assumptions regarding the interior structure, it
is clear that substantial compensation is required. Next, we explore
three-layer models that allow for compensation in either the Airy sense(lateral thickness variations in a shell overlying a weaker layer) or the
Pratt sense (lateral density variations in the outermost portion of thebody).
3.2.2. Three-layer model with Airy compensation—sensitivity to methods
For our three-layer Airy isostatic model, there are four free parameters: the density and mean thickness of the ice shell, and the density
and mean thickness of the underlying lower viscosity layer (which wecall the ocean); the core radius and density are derived from the known
total radius and bulk density. Although this parameter space has four
dimensions, we will begin here by fixing the ice shell and ocean densities, reducing the number of free parameters to two.
The results of this type of analysis are sensitive to certain methodological choices including the way the hydrostatic equilibrium figure
is computed, the way isostatic support is modeled, and whether or not
the finite amplitude of the topography is taken into account. We illustrate the effects of these choices by comparing our results against results
obtained under different assumptions and in previous work. Fig. 3a, forinstance, follows the assumptions of Iess et al. (2014) in which theequilibrium figure is computed according to Eq. (11), placing it on thedotted line in Fig. 1. Although Fig. 3a does not resemble any figure in
Iess et al. (2014), the information is essentially the same, with the exceptions that we have re-referenced the gravitational potential coefficients to a radius of 252.1 km and, for the sake of uniformity across thepanels of Fig. 3, we assume here that ρshell = 925 kg/m3 andρocean = 1020 kg/m3. As discussed by McKinnon (2015), Enceladus's
fast rotation causes its equilibrium figure to depart slightly from the

3.2. Parameter space exploration
Here, starting from the observed shape, we consider a range ofpossible interior models in terms of the number of layers, layer thicknesses, layer densities, and compensation mechanisms. For each model,
we compute the corresponding gravitational potential using Eq. (8) andcompare the result to the observed values (Iess et al., 2014), as described in section 2.6.
3.2.1. Two-layer model
For our two-layer model (core and shell), there are only two freeparameters: the density and mean thickness of the ice shell. Having
fixed these two values, the density and mean radius of the core areconstrained by the known total radius and bulk density. We set the coreshape to a hydrostatic equilibrium figure (section 2.2) and we set theexterior to match the observed shape (Nimmo et al., 2011 or Tajeddine
et al., 2017). We then vary the two free parameters over a considerable

Table 2Preferred interior model parameters assuming Airy compensation, with the ice shell supported buoyantly in a subsurface ocean. The table shows the results for thenominal densities (corresponding to Fig. 4) as well as the results emerging from the four-dimensional parameter space exploration (corresponding to Fig. S1), in
which the layer densities were allowed to vary. Since the ice shell and ocean densities could not be constrained in that case, we show the full ranges considered anddesignate them with an asterisk (*). For all other parameters, the reported ranges correspond to the 68% confidence ranges emerging from the parameter space
exploration. For the ice shell thicknesses, we also show the thickness found at the south pole (SP). The minimum misfit (Lmin) across the parameter space is also given
for each case (see section 2.6).
2D grid with nominal densities

4D grid with density ranges

Shape model

Layer

Mean thickness

Density

Mean thickness

Density

Nimmo et al. (2011)

Ice shell

21 − 25 km(0 − 3 km at SP)
37 − 44 km
185 − 192 km
Lmin = 1.6
19 − 24 km(7 − 11 km at SP)
35 − 39 km
192 − 195 km
Lmin = 3.7

925 kg/m3

19 − 24 km(4 − 12 km at SP)
37 − 48 km
184 − 193 km
Lmin = 0.6
21 − 26 km(9 − 17 km at SP)
30 − 39 km
191 − 198 km
Lmin = 3.2*850 − 950 kg/m3

Ocean
Core
Tajeddine et al. (2017)

Ice shell
Ocean
Core

1020 kg/m3
2420 − 2590 kg/m3
925 kg/m3
1020 kg/m3
2340 − 2410 kg/m3

119*1000 − 1100 kg/m3
2390 − 2570 kg/m3*850 − 950 kg/m3*1000 − 1100 kg/m3
2320 − 2410 kg/m3

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

90°N
60°N

a

b

5
10

15

15

30°N
20
25

30

25

25

20

20

15
10
5

60°S
90°S
90°N

20

10

c

20

d

15

15
10

15

60°N

25

15

20

20

25

25

30
25

30°S

0

25

30

30

0°

5
25

25

25

30°N

Ice Shell Thickness (km)

30

35

25

35

30

25

25

30°S

30

25

30

30

30

20

25
35

2
35 30 5

0°

60°S

25

25

20

20

15
10

90°S
0°

15

5

45°E

90°E

135°E

180°E

225°E

270°E

315°E

0°

45°E

90°E

135°E

180°E

225°E

270°E

315°E

360°E

Fig. 5. Lateral shell thickness variations implied by Airy compensation assuming the topography is supported isostatically. Representative examples are shown basedon the shape models of (a,b) Nimmo et al. (2011) and (c,d) Tajeddine et al. (2017), where both expansions are limited to spherical harmonic degree 8. Left panels(a,c) correspond to the nominal densities of ρshell = 925 kg/m3 and ρocean = 1020 kg/m3 whereas right panels (b,d) correspond to a larger density contrast, withρshell = 900 kg/m3 and ρocean = 1100 kg/m3, thus yielding more subtle shell thickness variations. The mean layer thicknesses are chosen to give the best fit for eachcombination of shape model and layer densities. The mean layer thicknesses (shell/ocean/core), in kilometers, and misfits [in square brackets], are: (a) 22/42/188[1.6], (b) 21/46/185 [0.6], (c) 21/37/194 [3.7], (d) 23/38/191 [3.2].

usual approximation—whereas it is usually assumed that
hyd
hyd
H20
/H22
= − 10/3, the appropriate value for Enceladus is ≈ − 3.25.
The latter value is obtained following the approach of Tricarico (2014).
Although the difference in this ratio is only ∼2.5%, the effect is asignificant change to the point at which the J2 (blue) and C22 (red)curves intersect—at a shell thickness of ∼49 km rather than ∼33 km,
when the Nimmo et al. (2011) shape model is adopted (Fig. 3b). Onceagain, our illustration does not resemble any figure in McKinnon(2015), but the calculations are the same, except for the assumed layer
densities, the effects of which are subtle. Figs. 3e and f show the corresponding results when the Tajeddine et al. (2017) shape model isadopted instead.
Unlike the work of Iess et al. (2014) and McKinnon, (2015), ourcurrent approach involves computing the hydrostatic equilibrium figures for all three layers (section 2.2). The hydrostatic equilibrium figureis then maintained for the core (a topic we return to in section 3.2.6)
while the shape of the outermost surface is set to match the observations (Nimmo et al., 2011 or Tajeddine et al., 2017), and the shape of
the ice/ocean interface is set to ensure that the topography is supported
via Airy isostatic compensation. Also unlike Iess et al., (2014) and
McKinnon, (2015), which effectively treat the topography as surface
density anomalies for purposes of computing the resulting gravitational
field, our current approach includes a correction accounting for thefinite amplitude of the topography (Wieczorek and Phillips, 1998).
Fig. 3c illustrates the results obtained when we compute the equilibrium figures using the numerical approach of Tricarico (2014) andwhen we apply the finite amplitude correction (Wieczorek and Phillips,
1998).
Both Iess et al. (2014) and McKinnon (2015) followed what we callthe “equal masses” approach to isostasy (e.g., Lambeck, 1988) in whichcolumns of equal width are required to contain equal masses. As discussed by Hemingway and Matsuyama, (2017), however, this does not
eliminate lateral pressure gradients at depth, and therefore may notcorrespond to an equilibrium condition (see also Beuthe et al., 2016).

Whereas the “equal masses” approach to isostasy was maintained forFig. 3c, we instead advocate the “equal pressures” approach to isostasy(Hemingway and Matsuyama, 2017), wherein the goal is to eliminate
lateral pressure gradients at depth (see discussion in section 2.3). The
effect is dramatic (Fig. 3d), with the J2 (blue) and C22 (red) curves now
intersecting at a shell thickness of ∼30 km, when the Nimmo et al.
(2011) shape model is adopted. Moreover, whereas the J3 (gold) curve
was not sensitive to the way in which the equilibrium figure wascomputed (compare Fig. 3a and b), it is strongly affected by the change
in the approach to isostasy, shifting from a preferred shell thickness ofmore than 40 km (Fig. 3c, “equal masses”) to just ∼20 km (Fig. 3d,“equal pressures”). Figs. 3g and h illustrate the corresponding results
when the Tajeddine et al. (2017) shape model is adopted instead.
Fig. 3c, d, g, and h additionally illustrate the parts of the parameter
space that best satisfy the observed libration amplitude (purple line),demonstrating that the libration constraint agrees best with the gravity
constraints when the “equal pressures” model is adopted, consistent
with the results of Beuthe et al. (2016). This is also evident from theminimum misfit values provided in the caption. From here on, our results employ the “equal pressures” isostasy model exclusively.
3.2.3. Three-layer model with Airy compensation—results
Fig. 4 shows, for one example combination of ice shell and ocean
densities (ρshell = 925 kg/m3 and ρocean = 1020 kg/m3), the ranges ofshell and ocean thicknesses that yield, via Eq. (8), gravitational fields
that match the observed J2 (blue), C22 (red), or J3 (gold) (Iess et al.,2014), as well as those thicknesses that yield, via Eq. (20), a libration
amplitude that matches the observed value (purple). This two-dimensional slice of the parameter space serves to illustrate how the model
gravitational potential coefficients and libration amplitude vary as afunction of shell and ocean thickness. All three gravity coefficients are
sensitive to the assumed shell thickness because this determines how
effectively the isostatic compensation reduces the magnitude of the
non-hydrostatic gravity arising from the non-hydrostatic topography.
120

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

9

7

10

40 8
9

9

10

30

7

8

8
9

20

10
9

10

10

7

20

30

40

50

60

9
8

40
8
190

30

10

20
10

10

0

10

20

30

40

20

2600
2500
2400
2300
2200
2100

0

10

20

30

40

50

60

70

Mantle Thickness (km)

g

50
40
30
20

50

60

70

180

190

200

210

220

230

Core Radius (km)

h

2700

10

Mantle Thickness (km)

2000
170
2800

60

Compensation Depth (km)109Compensation Depth (km)

8
190

30

70

f

60
50

40

70

Mantle Thickness (km)

70

50

8

10

d

2700

10

8
9

0

e

8

7

c

60

Core Density (kg/m 3)

50

7Compensation Depth (km)

60

2800

70

b

10

Compensation Depth (km)

a

8

7

Core Density (kg/m 3)

70

2600
2500
2400
2300
2200
2100

0

10

20

30

40

50

Mantle Thickness (km)

60

70

2000
170

180

190

200

210

220

230

Core Radius (km)Fig. 6. Misfits and probability densities across the parameter space assuming a three-layer model with Pratt compensation, assuming mean densities ρshell = 925 kg/m3
and ρmantle = 1030 kg/m3, and comparing results derived from the shape models of (a–d) Nimmo et al. (2011) and (e–h) Tajeddine et al. (2017). (a,e) Lines indicate
combinations of compensation depth and mantle thickness that satisfy the observed J2 (blue), C22 (red), or J3 (gold) dimensionless gravitational potential coefficients(Iess et al., 2014), with shaded bands indicating 1σ (dark) and 2σ (pale) bounds from combined shape and gravity model uncertainties (section 2.6). (b,f) Misfit between
the model and the observations, given by Eq. (21). (c,g) 68% (dark), 95% (intermediate), and 99.7% (pale) confidence contours showing compensation depths andmantle thicknesses that best account for the observed gravitational potential, and (d,h) corresponding core radius and density (section 2.6). The minimum misfits are
L = 6.0 and L = 7.0 for the Nimmo et al. (2011) and Tajeddine et al. (2017) shape models, respectively. (For interpretation of the references to color in this figurelegend, the reader is referred to the web version of this article.)

Only the two degree-2 coefficients, however, are also sensitive to theassumed ocean thickness. This is because the magnitude of the hydrostatic component of those coefficients depends on the overall radial
density structure of the body. Since there is no hydrostatic component
to the J3 term, it is not sensitive to ocean thickness, except insofar as thehydrostatic components affect the total amplitude of the topography,
thereby affecting the potential via the finite amplitude correction(Wieczorek and Phillips, 1998). Similarly, the libration amplitude ismainly sensitive to the mean shell thickness because this directly determines the shell's moment of inertia (Cs). The libration amplitude isalso weakly sensitive to the mean ocean thickness because this slightly
alters the moment of inertia of the core and also affects the decoupledsystem's natural frequencies as well as the torques acting on the different parts of the system (section 2.5). Because the libration amplitude
is an inverse function of shell thickness, it increases rapidly at small
shell thicknesses, effectively providing a strong lower bound on shell
thickness—this is evident in the shapes of the misfit contours in Fig. 4b
and f.
The misfits are larger when the Tajeddine et al. (2017) shape modelis adopted instead of that of Nimmo et al. (2011), primarily because of
the much smaller uncertainties. As we discussed in section 3.1, however, the two studies took very different approaches in their analyses of
uncertainties such that the two results may not be directly comparable.
In principle, the minimum misfit of L = 3.7 (Fig. 4f) can be interpreted
as an indication that there are no points in the parameter space thatsimultaneously satisfy all of the observations within their 1σ uncertainties and that additional model parameters may be required toaccommodate the observations. Although we do explore additionalmodel parameters below (sections 3.2.5 and 3.2.6), we emphasize thatthe meaning of the uncertainties reported by Tajeddine et al. (2017) is
not entirely clear and should perhaps be interpreted with caution. As
discussed in section 2.6, since we cannot evaluate the likelihood of any

particular model in absolute terms, we instead focus on the shapes ofthe probability density functions to indicate, within the context of eachmodel, which model parameters fit best (e.g., Fig. 4c, d, g, and h).
Whereas Fig. 4 illustrates only a two dimensional slice of the parameter space, with a fixed combination of shell and ocean densities, weconsider a wide range of possible densities for both the shell and theocean, leading to a four dimensional parameter space. We compute themisfit and use it to compute the probability density everywhere acrossthis four dimensional parameter space (section 2.6), subject to the constraint that the shell and ocean thicknesses must be everywhere greaterthan zero. Fig. S1 shows confidence contours illustrating the ranges ofmodel parameters that best account for the observed gravitational potential and libration amplitude. The preferred shell thickness is notsensitive to the assumed densities for the ice shell and ocean, thoughslightly better fits tend to occur for smaller ice shell densities and larger
ocean densities. The preferred ocean thickness also increases slightly
with both ice shell density and ocean density, but neither of these densities are well constrained. Though there is considerable overlap between
the confidence contours arising from the two different shape models, thenewer shape model (Tajeddine et al., 2017) tends to favor slightly larger
ice shell thicknesses and smaller ocean thicknesses (Table 2).
Although Figs. 4 and S1 refer only to the mean shell thickness, Airycompensation entails considerable lateral shell thickness variations.
Fig. 5 illustrates, for four example interior models, the lateral shell
thickness variations assuming that all of the topography is supported by
Airy isostasy. The four cases correspond to the best fitting interior
models for each shape model given two different sets of layer densities(a nominal reference set of densities and a high density contrast case).
The amplitude of the shell thickness variations is a function, via Eq.
(12), of the density contrast between the ice shell and the ocean, withlarger density contrasts corresponding to subtler shell thickness variations (Fig. 5b and d). The shell is always thinnest at the poles and121

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

b

960

0

940

900

93

0

0

0

920

0

93

0

940

1050

93

0

94

0

930

920

920

95

940
950
960

60°S

930
0
92

93

30°S

910

900 0
91

91

0

91

9

910 0
90

0

910

89

0°

93

960

30

0

0

960
970
980
990

95

1000

970

c

960

d

930

940

950
960

940

60°N

97

920
910

30°N

90°S
90°N

970

980

Crust Density (kg/m 3 )

60°N

950

a

95 0
94
0

90°N

30°N
0
0

60°S

0

95

0

93
0

94

960

940

970

950

90°S
0°

45°E

0

0

92

0

900

94

920

95

93

940

0

94

910

910

30°S

94
0

92

950

0

93

0°

90°E

135°E

180°E

225°E

270°E

315°E

0°

45°E

90°E

135°E

180°E

225°E

270°E

315°E

360°E

Fig. 7. Lateral shell density variations implied by Pratt compensation assuming all of the model topography is supported isostatically. Representative examples areshown based on the shape models of (a,b) Nimmo et al. (2011) and (c,d) Tajeddine et al. (2017). Left panels (a,c) correspond to the nominal densities of ρshell = 925
kg/m3 and ρmantle = 1030 kg/m3 whereas right panels (b,d) correspond to larger densities, with ρshell = 950 kg/m3 and ρmantle = 1300 kg/m3. The approximate mean
layer thicknesses (shell/mantle/core), in kilometers, and misfits [in square brackets], are: (a) 39/21/192 [6.0], (b) 38/35/179 [6.0], (c) 49/7/196 [7.0], (d) 49/13/
190 [7.0].

a

b

c

d

122

Fig. 8. Effect of elastic flexural support on thebest fitting shell thickness under the assumptionof either (a,c) top loading; or (b,d) bottom
loading, and adopting the shape models of (a,b)
Nimmo et al. (2011) or (c,d) Tajeddine et al.
(2017). Lines indicate combinations of ice shell
thickness (i.e., compensation depth) and elastic
layer thicknesses (Te) that satisfy the observed
J2 (blue), C22 (red), or J3 (gold) dimensionlessgravitational potential coefficients (Iess et al.,2014), and the observed libration amplitude(purple) (Thomas et al., 2016), with shaded
bands indicating 1σ (dark) and 2σ (pale) uncertainties on the estimated values (section 2.6).
The minimum misfits are: (a) L = 1.9, (b)L = 1.9, (c) L = 3.4, and (d) L = 3.7. (For interpretation of the references to color in this figurelegend, the reader is referred to the web versionof this article.)Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

thickest at the equator, especially near the prime- and anti-meridians.
The thinnest parts of the shell correspond to the topographic depressions near the poles, with the thinnest part of the shell coinciding with
the South Polar Terrain (SPT). We emphasize that, although we show
only four examples in Fig. 5, a range of interior models and shell
thicknesses are compatible with the observational constraints (Fig. S1,Table 2).

shell and mantle, though slightly better fits tend to occur for larger shelland mantle densities. Owing to the relatively smaller uncertainties in
the Tajeddine et al. (2017) shape model, the confidence contours enclose a much smaller part of the parameter space, with a strong preference for large shell thicknesses and small mantle thicknesses.
So far, we have discussed only the mean density of the outer ice
shell, but the Pratt compensation mechanism entails lateral variationsin the ice shell density. Fig. 7 illustrates, for four example interior
models, the required lateral density variations under the assumption of
Pratt isostatic compensation. The four cases correspond to the best
fitting interior models for each shape model given two different sets ofmean layer densities (a nominal reference set of densities and the
maximum density case). When the density variations occur in a 40 km
thick compensating layer, Pratt compensation implies 4 − 5% higher
density ice in the South Polar Terrain (SPT) and in the heavily crateredregion on the eastern hemisphere at high northern latitudes. When thecompensating layer is only 20 km thick, the required density anomaly is
8 − 10%.

3.2.4. Three-layer model with Pratt compensation
For our three-layer Pratt isostatic model, there are four free parameters: the mean densities and thicknesses of the two outer layers; thecore radius and density are derived from the known total radius andbulk density. We assume the outermost layer is composed of water ice
with a nominal mean density of 925 kg/m3 (though we consider a range
between 850 and 950 kg/m3). Pratt compensation is accomplishedthrough relatively smaller lateral density variations in this outer layer.
We refer to the mean thickness of the outer layer as the compensation
depth, since the compensating density anomalies occur throughout thislayer. We refer to the intermediate layer between the ice shell and thecore as the mantle and we assume it is part of the primarily water ice
envelope surrounding the core, and hence its density is only slightly
higher than that of the overlying ice shell (we consider mantle densities
in the range of 1000 − 1300 kg/m3). As with the Airy model, we begin
by computing the hydrostatic equilibrium figures for all layers (section
2.2). The hydrostatic equilibrium figure is then maintained for the core(see section 3.2.6) while the shape of the outermost surface is set tomatch the observations (section 3.1). Finally, the required lateral
density anomalies, computed via Eq. (13), are applied to the outer layer
in order to achieve Pratt isostatic compensation (section 2.3;
Hemingway and Matsuyama, 2017).
Fig. 6 shows, for one example combination of mean shell and mantle
densities (ρshell = 925 kg/m3 and ρmantle = 1030 kg/m3), the ranges ofshell and mantle thicknesses that yield, via Eq. (8), gravity fields thatmatch the observed J2 (blue), C22 (red), or J3 (gold) (Iess et al., 2014).
This two-dimensional slice of the parameter space serves to illustrate
how the model gravitational potential coefficients vary as a function ofcompensation depth and mantle thickness. As with the Airy compensation case, all three coefficients are sensitive to the assumed compensation depth because this determines how effectively the isostaticcompensation reduces the magnitude of the non-hydrostatic gravity
arising from the non-hydrostatic topography. Only the two degree-2
coefficients, however, are also sensitive to the assumed mantle thickness. This is because the magnitude of the hydrostatic component of
those coefficients depends on the overall radial density structure of the
body. As with the two-layer model, for the three-layer Pratt model,
where it assumed that the layers are physically coupled (cannot easily
slide against each other), the maximum libration amplitude, computed
via Eq. (17), is ∼0.04°, more than 5σ below the observed value of0.120 ± 0.014°, meaning that no part of this parameter space is compatible with the observed libration amplitude (Thomas et al., 2016).
Consideration of the Pratt model nevertheless remains useful for determining the compensation depth that would be required to accommodate the gravity observations; in principle, a thin decoupling liquid
layer could be inserted between the mantle and core, permitting alarger libration amplitude.
Whereas Fig. 6 illustrates only a two dimensional slice of theparameter space, with a fixed combination of shell and mantle densities, we consider a wide range of possible densities for both the shelland the mantle, leading to a four dimensional parameter space. We
compute the misfit and use it to compute the probability density everywhere across this four dimensional parameter space (section 2.6),
subject to the constraint that the shell and mantle thicknesses must be
everywhere greater than zero. Fig. S2 shows 68%, 95%, and 99.7%
confidence contours illustrating the ranges of model parameters thatbest account for the observed gravitational potential. The preferredshell thickness is not sensitive to the assumed mean densities for the ice

3.2.5. Flexural support
If part of the ice shell behaves elastically, and if the shell thickness
variations were generated after the formation of this elastic layer, then,in addition to buoyant support, the topography will be supported inpart by bending and membrane stresses (e.g., Turcotte et al., 1981). In
such a situation, the basal topography cannot be computed with Eq.
(12) and must instead be computed using Eqs. (14) or (15), depending
on how the thickness variations were generated (see Appendix A). The
distinction between Eqs. (14) and (15) is important. For a given observed surface relief, in the case of loading (thickening/thinning) at thesurface, additional rigidity (increasing Te and thus ζ) means relatively
smaller topography at the base of the shell whereas, for basal loading,
additional rigidity means relatively greater basal topography. The
choice of top versus bottom loading thus drives the resulting basal topography amplitude in opposite senses, with the difference becoming
increasingly significant for increasingly stiff elastic shells (Fig. 8).
The assumption of top loading was adopted in a study by Čadek
et al. (2016), leading the best fitting ice shell thickness to decrease withincreasing shell rigidity until it came into agreement with the results
implied by the libration observations (Thomas et al., 2016). A similar
situation is illustrated in Fig. 8a and c, except that, as discussed above,
we are using a different approach to computing the equilibrium figures
and a different approach to isostasy. However, the assumption of top
loading is appropriate only if the topography arises due to surface
processes such as impacts or the emplacement of erupted materials. By
contrast, we regard tidal heating (and the associated asymmetries inmelting and freezing at the base of the ice shell) as the dominant sourceof long wavelength shell thickness variations. That is, we regard bottom
loading as the more appropriate model for the generation of long wavelength shell thickness variations on Enceladus. Accordingly, increasing elastic thickness translates to relatively larger basal topography and therefore requires larger mean shell thicknesses in order tosatisfy the observed gravity (Fig. 8b and d) (see also Beuthe et al., 2016,
Text S7).
For Fig. 8, in order to clearly illustrate how the non-hydrostaticparts of the gravitational field depend on the compensation depth (i.e.,shell thickness) and the elastic layer thickness, we fixed the moment ofinertia factor to an arbitrary value of 0.335, therefore maintaining afixed hydrostatic component across this parameter space. Because lines
of constant moment of inertia factor are nearly parallel to lines of
constant C22 in Fig. 4a and e, C22 varies only slightly across the parameter space illustrated in Fig. 8. Hence, the C22 gravity observation
does not provide a strong constraint on the elastic layer thickness, as
illustrated by the broad red shaded bands in Fig. 8. Due to the large
uncertainties in the Nimmo et al. (2011) shape model, the J2 term issimilarly of limited usefulness. The J3 term is more useful, as illustrated
by the relatively smaller gold bands in Fig. 8. Finally, the observed123

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

libration amplitude provides an independent constraint that agrees bestwith the gravity constraints when the elastic thickness is near zero (atthe left edge of the parameter space). The libration amplitude depends
weakly on the elastic layer thickness, mainly because of the effect on
the shape of the ice/ocean interface. We do not take into consideration
any dynamical effects that the elasticity of the ice shell might have on
the libration amplitude, and we do not expect such effects to be significant for Enceladus (see section 2.5; Van Hoolst et al., 2016).
We also carried out a three dimensional parameter space exploration, varying shell and ocean thickness as well as elastic layer thickness,while keeping the shell and ocean densities fixed at their nominal values of ρshell = 925 kg/m3 and ρocean = 1020 kg/m3. We did this forboth the top- and bottom-loading cases (Figs. S3 and S4). The three
dimensional parameter space exploration demonstrates that the best
fitting ice shell thickness increases with increasing elastic thickness
when bottom loading is assumed, and decreases with increasing elastic
thickness when top loading is assumed, consistent with the discussion
above and in Appendix A. In any case, the best fits always correspond toan elastic layer thickness ≤20 m—that is, simultaneously accommodating the gravity and libration observations becomes increasingly
difficult with increasing elastic layer thickness (Figs. S3 and S4). We
conclude from this that elastic support must not be effective for the
longest wavelength topography on Enceladus. This conclusion is perhaps not surprising given the heavily fractured south polar region,
where long wavelength bending and membrane stresses may be interrupted (e.g., Soucek et al., 2016; Beuthe et al., 2016). This is not to say
the elastic thickness must be everywhere negligible, as the cold upper
part of the ice shell should behave elastically in regions where it is not
too heavily fractured. But even on local scales, for example, the effective elastic layer thickness has been estimated at just ∼0.3 km based on
apparent plate flexure at a rift zone (Giese et al., 2008) or 0.4 − 1.4 km
based on a series of ridges and troughs under an unstable extension
model (Bland et al., 2007).

Nevertheless, in order to provide a sense of what core topography
would be required to achieve a perfect match to the observed gravitational field, we solved Eq. (8) for the required core topography given
the observed gravity (Iess et al., 2014), setting aside the finite amplitude correction for simplicity. When there is no isostatic compensation,some 1500 m of non-hydrostatic core topography is required (Fig. S5a)
to accommodate the observed gravity. When the topography is supported by Airy compensation as in section 3.2.3, where the smallness ofthe non-hydrostatic gravity is explained mainly by isostatic compensation, then just ∼350 m of non-hydrostatic core topography is sufficient (Fig. S5b) to accommodate the small remaining non-hydrostatic
gravity.
3.3. Summary
The above parameter space exploration demonstrates that, in orderto accommodate the observed shape and gravity, significant compensation is required (section 3.2.1). In principle, the observations could beaccommodated by either Airy or Pratt compensation (sections 3.2.3 and3.2.4, respectively), but the Airy model is preferred for several reasons.
First, the large amplitude of Enceladus's forced physical librations(Thomas et al., 2016) implies a global subsurface ocean, which implies
a floating ice shell. A floating ice shell is likely to experience laterally
variable tidal dissipation (see section 4), naturally leading to lateralshell thickness variations and therefore at least some degree of Airycompensation. The presence of a decoupling ocean in the Airy models
also permits large libration amplitudes and is one of the main reasonsthe misfits are much smaller in the Airy models (Fig. 4) than in theothers (Figs. 2 and 6).
Second, even if Pratt isostasy were to dominate within the floating
ice shell, it would require the systematic lateral density variations to be
distributed through a layer that is some 40 − 60 km thick (Fig. 6g)
when the Tajeddine et al. (2017) shape model is adopted, or 20 − 60
km thick (Fig. 6c) when the Nimmo et al. (2011) shape model isadopted. The libration amplitude, however, implies an ice shell thickness of roughly 16 − 22 km, depending on the ocean thickness andassumed layer densities (Fig. 4a, e; cf. Thomas et al., 2016; Van Hoolst
et al., 2016). Whereas this range is marginally compatible with thegravity results when the older shape model (Nimmo et al., 2011) is
adopted, the newer shape model (Tajeddine et al., 2017) leads to a
result that requires an ice shell roughly twice as thick as that implied by
the libration amplitude.
Third, although the amplitude of the Pratt-required systematic lateral density variations is modest, and could in principle be accounted
for by variations in porosity (Besserer et al., 2013), the spatial pattern(Fig. 7) is only partially consistent with the geological constraints (e.g.,
Crow-Willard and Pappalardo, 2015). Specifically, the Pratt mechanism
requires relatively higher densities in the south polar region as well as
at high northern latitudes on the trailing hemisphere. Whereas thegeologically active south polar region is characterized by a youthful
surface and may well have experienced thermally driven compaction(Besserer et al., 2013), the high northern latitudes are instead characterized by an ancient and heavily cratered surface with no signs ofanomalous heating or porosity reduction that could lead to a positive
density anomaly. This is not to say that lateral density variations do notexist. Indeed, compaction may be the best way to account for some of
the large (but not global) scale topographic basins (e.g., those discussed
by Nimmo et al., 2011; Schenk and McKinnon, 2009; Tajeddine et al.,2017). Some hybrid of Airy and Pratt support is of course possible at thelongest length scales as well. However, it is clear that the Pratt mechanism cannot dominate.
For our Airy compensation model, the primary parameters are theice shell and ocean thicknesses while the core radius and density arederived parameters. We assume a hydrostatic figure for the core because, although contributions from core topography cannot be ruled
out, they are not expected to bias the results one way or the other

3.2.6. Core topography
Our analysis has assumed that the shape of the core corresponds tothe expectation for hydrostatic equilibrium. If the core is weak (e.g.,
Roberts, 2015), as expected if it consists of fluid-saturated unconsolidated materials (i.e., a mud-like core), or even if it is presently
strong but was weak for an extended period early in its history, it may
be reasonable to expect the core to conform to a hydrostatic equilibrium
figure. If, however, the core exhibits sufficient strength to maintain adeparture from the expected hydrostatic equilibrium figure, then thisnon-hydrostatic core topography will contribute to the observed gravitational field (e.g., McKinnon, 2013; Tajeddine et al., 2014). And because of the relatively large density contrast between the core and theoverlying water/ice, a significant contribution can be made with relatively little core topography.
While it is always possible to arrange such unseen core topography
so as to aid in reducing the observed non-hydrostatic gravity coefficients (i.e., compensating for the non-hydrostatic topography), there isno physical reason to expect it to be so arranged—effectively mirroring
the long wavelength non-hydrostatic topography in the same way thatisostatic roots would, resulting in a reduction in the non-hydrostaticpart of J2. If anything, core topography might be expected to make apositive contribution to the non-hydrostatic part of J2 since Enceladus
was likely spinning even faster (experiencing greater centrifugal flattening) earlier in its history when the core was forming. Since both themagnitude and sign of any non-hydrostatic core topography are unknown, we adopt the simplest assumption, which is that the core shape
matches the expectation for hydrostatic equilibrium. In principle, theunknown core topography could be taken into account using additional
dimensions in the parameter space, with some a priori assumptions
about the reasonable range of possibilities. Doing so would have theeffect of widening the probability density functions shown above, butshould not bias the best fitting results one way or another.
124

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal(section 3.2.6). We additionally take the ice shell and ocean densities tobe free parameters, but their effects are so subtle that they cannot bemeaningfully constrained by shape, gravity, and libration observations
alone (Fig. S1). The dominant terms of the quadrupole gravity field, J2
and C22, effectively constrain the thickness of the H2O envelope, andtherefore the core radius and density, but allow for trade-offs in how
this envelope is partitioned into its solid and liquid components(Fig. 4a, e). By contrast, the zonal part of the octopole field, J3, is insensitive to ocean thickness and therefore better constrains ice shell
thickness (Fig. 4a, e), underscoring the importance of obtaining good
observational constraints on the higher order terms of the gravitational
field. The libration amplitude is likewise strongly sensitive to shell
thickness but is only subtly and indirectly affected by ocean thickness(Fig. 4a, e; section 2.5). Expanding the parameter space exploration toallow for a range of possible ice and ocean densities does not significantly alter the best fitting parameters. This additional parameter
space does, however, permit larger density contrasts and therefore
subtler shell thickness variations, allowing for larger south polar shell
thicknesses (Table 2).
Assuming the long wavelength topography is primarily the result ofasymmetric ocean freezing, the inclusion of an elastic layer would mutethe surface expression of the lateral shell thickness variations, implying
larger isostatic roots (section 3.2.5). In that case, a larger mean shell
thickness would be required to accommodate the observed gravitational potential (cf. Čadek et al., 2016). Such models, however, lead tolarger misfits compared with the observations and, in some cases,
grounding of the ice shell onto the core (as the mean ocean thickness
becomes too small to accommodate the considerable topography at theice-ocean interface)—a situation that is precluded by the large physicallibration amplitude. Furthermore, significant elastic support requiresthat large membrane stresses be transmitted across the heavily fractured South Polar Terrain (SPT), which may be difficult (Soucek et al.,2016) given that the Tiger Stripe fractures likely penetrate the ice shellcompletely.
Adopting a pure Airy isostasy model then, and focusing on the results based on the new shape model (Tajeddine et al., 2017), the preferred mean shell and ocean thicknesses are 19 − 24 km and 35 − 39
km (68% confidence ranges), respectively (Fig. 4g) and the corresponding core radius and density are 192 − 195 km and 2340 − 2410
kg/m3 (68% confidence ranges) (Fig. 4h), a density suggestive of hydrated silicates and/or some degree of porosity in the core (Table 2).
Importantly, the shell thickness ranges from typically ∼35 km near theprime- and anti-meridians to as little as ∼5 km at the south pole(Fig. 5c). We emphasize that a range of different mean shell thicknesses
are permitted by the observations and that this range should be taken
into consideration in any attempt to understand the evolution of the ice
shell, its response to tidal stresses, and the nature of ocean to surface
pathways. This is particularly important for the south polar region
where the ice shell thickness could be anywhere from 4 km to 14 km(our 95% confidence range), assuming nominal ice and ocean densities
and adopting the newer shape model (Tajeddine et al., 2017).

scale ice shell convection to be unlikely and we focus instead on conductive heat loss. Conductive heat loss is an inverse function of shell
thickness and, focusing on the result based on the new shape model(Tajeddine et al., 2017), ranges from ∼20 to 35 GW, where we haveconsidered ocean temperatures in the range of 220 K to 273 K (Fig. 10).
This calculation does not include the additional advective heat loss
associated with the eruptions in the South Polar Terrain (SPT), which
may be in the vicinity of ∼5 GW (Howett et al., 2011; Spencer et al.,2013). Hence, within our 68% confidence range of likely mean shell
thicknesses, the implied total heat loss is in the range 25 − 40 GW.
By comparison, assuming a chondritic radiogenic heating rate of
4.5 × 10−12 W/kg (Spohn and Schubert, 2003) and a chondritic massfraction of ∼0.5 (obtained by assuming the chondritic component
density to be 3500 kg/m3 while the remainder has a density of 1000 kg/m3), radiogenic heating within the rocky core is <0.3 GW. Tidal dissipation is potentially greater as Enceladus's eccentric orbit is maintained by its 2:1 mean motion resonance with Dione, and Saturn'sabundant rotational energy is gradually transferred into the satellites.
Making traditional assumptions about Saturn's dissipation quality factor(Q ≈ 18,000), the long term average tidal dissipation within Enceladus
is, however, limited to just ∼1.1 GW (Meyer and Wisdom, 2007, 2008).
This significant shortage of heat supply would imply a rapidly
freezing ocean (e.g., Roberts and Nimmo, 2008). However, recent astrometric observations (Lainey et al., 2012, 2017) indicate that Saturn
may be effectively much more dissipative than previously thought, witha Q as small as perhaps ∼2000, permitting the tidal dissipation within
Enceladus to be as much as ∼25 GW; if Q < 2000, the dissipation
would be even greater still. Note that this low value of Q for Saturn may
be time- and frequency-dependent, and may vary from satellite to satellite as each satellite's tidal interactions with Saturn occur on different
timescales. Importantly, the resonance locking mechanism proposed by
Fuller et al. (2016) permits these low effective present-day values for Q,
and the corresponding large tidal dissipation rates, without requiring
the satellites to be young (see Nimmo et al., 2018 for a more thoroughdiscussion). Hence, the long-term average tidal dissipation rate within
Enceladus may be comparable to the inferred conductive heat loss,
permitting a balanced heat budget and an interior that is close to steadystate.
4.2. Spatial pattern
If the tidal dissipation occurs mainly near the base of the ice shell,
where the ice behaves viscously on a timescale comparable to the orbital period of 1.37 days, we can use a thin-shell tidal heating model(section 2.8) to compute the spatial pattern of tidal dissipation within
the ice shell (Fig. 11a). The corresponding equilibrium ice shell structure is characterized by being thinnest at the poles and thickest at theequator, especially at the prime- and anti-meridians (Fig. 11c), inagreement with the inferred shell structure (Fig. 11d). Note that thispredicted shell structure differs somewhat from that expected forEuropa (e.g., Ojakangas and Stevenson, 1989) where the shell thickness
is greater at the poles due to the low surface temperatures resulting
from the small solar obliquity. Assuming an ocean temperature of
273 K, our model gives a mean shell thickness of 21 km (matching theinferred mean thickness) when the total dissipation is ∼35 GW(Fig. 11b), corresponding to an interior characterized by k2/Q ≈ 0.02
(see section 2.8). Note that this total dissipation is a free parameter in
our ice shell tidal heating model, which we tune to achieve the desired
mean shell thickness (see discussion in section 4.3, below). That is, ourtidal heating model does not make a prediction about the mean shell
thickness, only the spatial pattern in the shell thickness variations.
Nevertheless, the broad agreement between the inferred shell structure
and the lateral thickness variations predicted by the ice shell tidal
heating model suggests that tidal dissipation near the base of the ice
shell may be an important source of internal heat production.
By contrast, assuming eccentricity tides, tidal dissipation within the

4. Internal heat production
Here we discuss what the inferred shell structure implies about theheat budget of Enceladus and how and where heat is dissipated within
its interior.
4.1. Heat budget
Given the range of possible mean shell thicknesses (Table 2), we canuse Eq. (35) to approximate the likely range of total conductive heatloss. If solid state convection were taking place in the ice shell, the heatloss would be much greater and maintaining significant lateral shell
thickness variations would be very difficult. Given that we infer significant lateral shell thickness variations (Fig. 5), we consider global
125Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

45
220 K
250 K
273 K

40

Shell Thickness (km)

35
30
25
20
15
10
5
0
Fig. 9. Depiction of the approximate interior structure of Enceladus, to scale,assuming the topography (Tajeddine et al., 2017) is supported by complete Airycompensation (Hemingway and Matsuyama, 2017), assuming nominal densitiesof ρshell = 925 kg/m3 and ρocean = 1020 kg/m3, and mean shell and ocean
thicknesses of 21 km and 37 km, respectively. The ice shell and ocean are cut
away between the prime meridian and 90°E. The corresponding core radius anddensity are approximately 194 km and 2400 kg/m3. Surface texture is a global
map of Enceladus produced by Paul Schenk (Lunar and Planetary Institute)from Cassini ISS data (NASA, JPL).

0

20

40

60

80

100

Heat Loss (GW)
Fig. 10. Mean shell thickness versus total heat loss assuming conductive
cooling and assuming a representative surface temperature of 60 K, and shown
for a series of different temperatures at the ice/ocean interface. The horizontal
purple shaded bands illustrate the 1σ (dark) and 2σ (pale) ranges of likely ice
shell thicknesses (Fig. 4g). The vertical gray band represents the internal dissipation permitted assuming the present-day effective dissipation within Saturn
corresponds to 1000 < Q < 4000, with the dashed line representing Q = 2000
(Fuller et al., 2016; Lainey et al., 2012). (For interpretation of the references tocolor in this figure legend, the reader is referred to the web version of thisarticle.)ocean produces an altogether distinct pattern, with maximum dissipation
at the equator, especially at the centers of the leading and trailinghemispheres, and minimum dissipation at the poles, leading to a shellthat is thinnest at the equator and thickest at the poles (Matsuyama et al.,2018). Moreover, such dissipation may only be a few ×10−6 W/m2(∼2 MW in total), some four orders of magnitude smaller than the value
corresponding to the inferred shell thickness (Fig. 11b). Obliquity tides
produce heating that is maximal at the poles but the magnitude is at least
four orders of magnitude too small, when adopting the parametersconsidered by Matsuyama et al. (2018). More significant dissipation
within the ocean may be possible, especially if turbulence is important at
the ocean boundaries where the topography may be rough (Wilson andKerswell, 2018), or due to resonant waves that may occur in the ocean(Lemasquerier et al., 2017), though further studies on these phenomena
are needed.
Tidal dissipation within the rocky core may occur in a pattern similar to what is expected for heating within the ice shell (Choblet et al.,2017), but that heat must then be transported through the ocean, which
could significantly alter or remove lateral variations in the heat supplied to the base of the ice shell. If the ocean is not stratified, narrowplumes may occur over the most strongly heated parts of the core,leading to localized ice shell thinning (e.g., Goodman et al., 2004;
Goodman and Lenferink, 2012; Vance and Goodman, 2009). If the moststrongly heated parts of the core are at the poles (Choblet et al., 2017),
this could help explain the thinner ice in the polar regions. Such polar
upwellings would not, however, account for the longitudinal variationsin the shell's structure. Upwellings could also be disrupted by any
convection that may be taking place within the ocean. And if convection in the ocean is sufficiently vigorous, the heat supplied to the ice
shell could in fact be concentrated at low latitudes (Soderlund et al.,2014), which would be inconsistent with the ice shell being thickest
near the equator. Although it is not known whether or not the ocean isstratified, the presence of salts and nanosilica grains in the ocean (Hsuet al., 2015; Postberg et al., 2011) along with the likely ongoingfreezing and melting at the ocean-ice interface (e.g., Čadek et al., 2017,2019) could lead to brine rejection at low latitudes and ponding of fresh

water at the poles. Such processes may lead to ocean stratification
which would inhibit direct radial heat transport between the core andthe ice shell. If the ocean is indeed stratified, the heat emerging fromthe core may be efficiently redistributed within the ocean, resulting in anearly uniform distribution of heat supplied to the base of the ice shell,
at least on the large scales we consider here, and therefore not contributing to global scale shell thickness variations. Hence, althoughsome degree of tidal dissipation probably occurs within the ocean andthe core, it is not clear that either of these modes of internal heatproduction could account for the inferred pattern of shell thickness
variations (Fig. 11d). We therefore argue that tidal dissipation within
the ice shell itself is likely an important component of internal heatproduction.
4.3. Discussion
Although our ice shell tidal heating model and our new interior
model exhibit broadly consistent ice shell structures, a few caveats
should be considered. Because our tidal heating model computes dissipation in a thin, symmetric shell, with dissipation being related to thesquares of the strain rates arising from the time varying degree-2 tidal
stresses (Ojakangas and Stevenson, 1989), the model predicts a symmetric pattern with power only up to spherical harmonic degree 4.
Consequently, the model does not predict many of the shorter wavelength features, nor the south polar anomaly. Even neglecting theseanomalies, our tidal heating model slightly underestimates the amplitude of the lateral shell thickness variations.
Our calculation of the equilibrium ice shell structure also neglects
the effect of any ongoing viscous relaxation. The lateral shell thickness
variations create lateral pressure gradients, which will cause secularflow in the low viscosity region near the base of the ice shell (e.g.,Čadek et al., 2017, 2019; Kamata and Nimmo, 2017), which will tend to
remove the lateral thickness variations. The resulting imbalance between conductive heat loss and local tidal dissipation will cause ocean

126

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

d

30

35

14

14

30

18

18

18

25

22

26

30°N

22

26

60°S

18

10

14

5

10
6

45°E

90°E

135°E

180°E

225°E

270°E

315°E

0°

45°E

90°E

135°E

180°E

15

22

26

22

18

26

26

26

30°S

20
22

30

30

26

22

26

30

30

22

30

22

26

30

0°

90°S
0°

5
45 0

50

0
35

225°E

270°E

315°E

360°E

Ice Shell Thickness (km)

30
35

c

10

26

60°N

20

60
70
80
100
120

60
65

90°S
90°N

30

45

30
35
40

55

60°S

40

35

50

40

35

30

40

30

45

50

40

35

40

30°S

60

40

35

40

40

0°

70
60

30

45

70

35

50

80

80

55

Heat Flux (mW/m2 )

60

60

30°N

80

b

65

60

60°N

a

70

90°N

0

Fig. 11. Comparison between predicted and inferred shell structure and heat flux. (a) Predicted pattern of tidal dissipation and (c) the corresponding equilibrium ice
shell structure when the total dissipation is ∼35 GW such that the mean shell thickness is 21 km. (b) Conductive heat loss corresponding to (d) the nominal inferredshell structure, with mean shell and ocean thicknesses of 21 km and 37 km, respectively (Fig. 4g), as derived from the newer shape model (Tajeddine et al. 2017),
truncated to spherical harmonic degree 8.

freezing beneath the thickest parts of the shell and basal melting at thethinnest parts of the shell. Taking this effect into account should further
mute the expected lateral shell thickness variations, making the large
amplitude of the inferred thickness variations yet more difficult to explain. On the other hand, our model also assumes that tidal dissipation
occurs in a uniform thin shell whereas, in reality, dissipation should beenhanced in the thinner parts of the shell near the poles (Běhounková
et al., 2017; Beuthe, 2018; Soucek et al., 2016)—a feedback that willlead to larger amplitude shell thickness variations. Similarly, because
our tidal dissipation model assumes no ongoing freezing and melting, it
is an approximation that is only valid if the rate of freezing and melting
is sufficiently slow—as it must be in order to maintain significant lateralshell thickness variations (see also Čadek et al., 2017, 2019; Kvorka
et al., 2018).
A related problem is that the ice must have a rheology that permits
the required ∼35 GW of heat dissipation without having a basal viscosity that is too low to maintain shell thickness variations. For example,
to dissipate 35 GW using the Efroimsky (2018) model (which, forsimplicity, assumes a uniform interior), and using an Andrade rheology(e.g., Shoji et al., 2013), would require a viscosity of ∼1013 Pa · s, depending on the assumed shear modulus. While such a viscosity is notunreasonable for water ice near its melting temperature, Čadek et al.
(2019) argue that if the basal viscosity is less than ∼3 × 1014 Pa · s,
viscous relaxation would be too rapid to be balanced by ongoingfreezing and melting. We note, however, that viscous relaxation andtidal dissipation occur on very different timescales (millions of years
versus 1.37 days) and given that the two processes may involve different mechanisms (e.g., dislocation creep versus grain sliding) andvery different strain rates, the relevant effective viscosities could besomewhat different (see also McCarthy and Cooper, 2016). An advantage of core dissipation models (e.g., Choblet et al., 2017), of course,is that high viscosities in the ice shell are not problematic because
dissipation in the ice shell is not required, as long as the core itself canbe sufficiently dissipative at tidal frequencies—a matter that still requires experimental validation (see Methods of Choblet et al., 2017). In
any case, further work is needed to identify realistic rheologies that candeliver the required dissipation, whether in the ice shell or the core,
within the context of a realistic three-dimensional model for Enceladus.
Finally, even when we assume that the tidal dissipation in the ice
shell is sufficient to produce a mean ice shell thickness of 21 km, theexpected dissipation pattern (Fig. 11a) slightly underestimates the heat
flux corresponding to the inferred shell structure (Fig. 11b) in certain
areas, such as at mid northern latitudes and particularly in the southpolar region, where the inferred heat flux is approximately twice what
is predicted by the model (Fig. S6). This could be explained in part by
the fact that our tidal dissipation model assumes a uniform shell
thickness whereas dissipation should be enhanced in the thinner partsof the shell (Beuthe, 2018; Soucek et al., 2016). Because we have set themagnitude of the dissipation to produce the desired mean shell thickness, there are naturally also regions where the model overestimates theheat flux slightly. It should be noted that the most significant excess
heating is confined to a region poleward of roughly 70°S, accounting forjust 3% of the surface area, such that the total heat anomaly illustrated
in Fig. S6 is only ∼3 GW. In reality, the conductive heat loss in thesouth polar region may be even less than 3 GW due to the elevatedsurface temperatures (Howett et al., 2011; Spencer et al., 2013) andcorrespondingly shallower geothermal gradients. Of course, the totalheat loss is likely somewhat greater, however, due to the advective heatloss associated with the eruptions.
5. Conclusions
Determining the structure of Enceladus's ice shell is challenging inpart because the analysis is sensitive to factors such as the approach tocomputing the reference hydrostatic equilibrium figure and the ways inwhich flexural and isostatic support are modeled (see section 3.2.2).
After an exhaustive parameter space exploration, and employing a
numerically accurate approach to calculating the equilibrium figures(section 2.2; Tricarico, 2014), we have obtained a new interior model(Figs. 4 and 9) that satisfies the observational constraints using the
Hemingway and Matsuyama (2017) approach to Airy isostasy (section2.3), which is meant to eliminate lateral pressure gradients at depth,127

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

and without the need for elastic support (sections 2.4 and 3.2.5; cf.
Čadek et al., 2016).
The ice shell is sufficiently thin that conductive heat loss significantly exceeds the already large ∼5 GW associated with the eruptions along the Tiger Stripe fractures (Howett et al., 2011; Spencer
et al., 2013). We find that global conductive heat loss alone is in therange of 20 − 35 GW, far in excess of radiogenic heating and tidal
dissipation when traditional assumptions about tidal interactions with
Saturn and Dione are adopted (Meyer and Wisdom, 2007, 2008).
Whereas this would imply a rapidly freezing ocean, recent astrometricobservations (Lainey et al., 2012, 2017) suggest that much greater tidal
dissipation is possible (∼25 GW assuming QSaturn = 2000) and, allowing for the resonance locking mechanism of Fuller et al. (2016), is
compatible with long-lived Saturnian satellites. Hence, although a balanced heat budget is not assured, it is at least possible, meaning thatthe interior of Enceladus could be close to a steady state, permitting along-lived internal liquid water ocean.
The implied shell structure (Fig. 11d) resembles the spatial pattern
expected if tidal dissipation occurs primarily near the base of the ice
shell (Fig. 11c). By contrast, tidal heating within the subsurface ocean(Matsuyama et al., 2018) may not be able to account for the observed
shell thickness variations both because the spatial pattern does notmatch and because the amplitude is orders of magnitude too small(section 4.2). Whereas tidal dissipation within the core could—if theocean is not stratified—contribute to ice shell thinning at the poles(Choblet et al., 2017), it cannot easily account for the inferred longitudinal variations in shell thickness, unless the core heat can be effectively transported radially outward through the ocean without

inducing ocean currents that would alter the spatial pattern. We
therefore conclude that tidal dissipation within the ice shell itself must
be an important mode of internal heat production (section 4.2).
Although our tidal heating model and our new interior model exhibit broadly consistent ice shell structures, future work is required todetermine the likely equilibrium between the opposing effects of viscous relaxation (e.g., Čadek et al., 2017, 2019) and enhanced dissipation in the thinner parts of the shell (Běhounková et al., 2017; Beuthe,2018) (section 4.3). Future measurements of heat flux outside the SPT,
observations of higher order terms in the gravitational field, and measurement of gravitational field variations over orbital timescales, could
help to further constrain the ice shell's structure and its implications forthe past and future evolution of Enceladus.
Acknowledgments
This research was funded primarily by the Miller Institute for Basic
Research in Science at the University of California, Berkeley, and wasmade possible by the NASA/ESA Cassini mission to Saturn and, inparticular, the work of the Radio Science and Imaging Science
Subsystem teams. We thank Francis Nimmo, Ondřej Čadek, and Anton
Ermakov for their thorough and constructive reviews. This work also
benefited from discussions with Michael Beuthe, Bruce Buffett, Michael
Manga, Bill McKinnon, Naor Movshovitz, Radwan Tajeddine, andGabriel Tobie. We thank Daniele Durante and Luciano Iess for supplying the covariance matrix corresponding to the estimated gravitational potential coefficients. All other data are publicly available in the
published literature.

Appendix A. Flexural support
We obtain Eqs. (14) and (15) following previous work (Kraus, 1967; Turcotte et al., 1981), but with several assumptions relaxed or altered tobetter suit the present problem. Treating the lithosphere as a uniform thin elastic shell (Kraus, 1967; Turcotte et al., 1981; Willemann and Turcotte,1982), we can write the pressure (pw) associated with the elastic restoring stresses resulting from the deflection of the shell (w) asD

6 w + 4D 4 w + 4D

2 w + ET R2e2 w + 2ET R2w
e

R4 (

2 + 1

(37)

) pw = 0

where D represents flexural rigidity and is given by

D=

ETe3
2)
12(1

(38)
2

where Te is the effective elastic thickness, E is Young's modulus, and ν is Poisson's ratio, and where ∇ is the Laplacian operator. R is the radius of thethin elastic shell, which we take to be the body's outer radius. Here we have taken w as positive upward deflection of the elastic shell, and pw as
positive downward pressure. Hence, pw and w must have the same sign in order to satisfy Eq. (37); positive upward deflection results in a positive
downward restoring pressure. Note that Turcotte et al. (1981) use a different convention in their Eq. (1) with the deflection, w, being positive
downward and with their pressure term referring to the applied pressure rather than the elastic response of the shell. Note also that, as pointed out byWillemann and Turcotte (1982), Turcotte et al. (1981) neglected to include the third term in Eq. (37), which has a small effect on our Eq. (43).
To balance the pressures at depth beneath unloaded and loaded portions of the shell, we use Eq. (3) from Hemingway and Matsuyama (2017),
except that we additionally include the effects of the elastic restoring stress. The result is

pw =

(

gb hb +

(39)c gt ht )where, ht and hb are, respectively, the (positive upward) topographic relief at the top and bottom of the ice shell, measured with respect to the
hypothetical equilibrium figure (i.e., ht and hb are the non-hydrostatic topography), where ρc is the density of the shell, Δρ is the density contrast
between the shell and the ocean, and gt and gb are the mean gravitational accelerations at the top and bottom of the shell, respectively.
Eq. (39) differs from Eq. (3) of Turcotte et al. (1981) in three ways. First, there are differences of sign convention, as discussed. Second, we
account for the possibility that the gravitational acceleration may vary between the top and bottom of the shell (this can be important for Enceladus;
see Hemingway and Matsuyama, 2017). Finally, we do not include the effects of self-gravitation. While self-gravitation of the relief does affect the
local equipotential surfaces, it does not make sense to alter Eq. (39) by the elevation of the geoid at the surface. As pointed out by Hemingway andMatsuyama (2017), the goal is not to balance the weight of surface topography against the buoyancy of basal topography. Instead, what is importantis the accumulation of pressure at depth computed via integration of the hydrostatic equation. A datum equipotential surface at depth will experienceelevation changes due to the surface relief and the compensating basal relief (only the former of which is addressed by adding the geoid elevation to
Eq. (39)). But at long wavelengths, even with significant elastic support, the compensation is large enough that the effects of the surface and basal
relief largely cancel, meaning that the datum equipotential surface elevation can be taken as fixed, and can therefore be neglected in Eq. (39).
In this model, the topographic relief at the top and bottom of the ice shell (ht and hb, respectively) are a function of the thickening or thinning of
the shell and the resulting lithospheric deflection (w). Thickening or thinning can occur due to: excavation of surface material due to impacts;
sediment transport; emplacement of surface loads due to volcanism; viscous relaxation; tectonic spreading or compression; or freezing or melting at
128

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

the base of the ice shell. If the load thicknesses at the top and bottom of the shell are dt and db, respectively, with positive values indicating thickening
and negative values indicating thinning, then (neglecting radial compression of the loaded shell) the topographic relief at the top and bottom of theshell can be written

h t = w + dt

(40)

hb = w(41)db.

For clarity and simplicity, we have not included the subscripts lm in the above expressions but these should be taken as the amplitudes of periodic
functions and the solutions here apply to individual spherical harmonic degrees and orders. We also assume that the periodic loads, dt and db, and theresulting lithospheric deflection, w, and topographic relief, ht and hb, are all in-phase—a good assumption when the loading is either purely from the
top or purely from the bottom.
Making use of the geometric constraints of Eqs. (40) and (41), and the requirement of equal pressures at depth (39), we can solve Eq. (37) for the
deflection, w, obtaining

w=

gb db

c gt dt

gb +

(42)

c gt +

where we have introduced an additional parameter, ζ, which has units of force per unit volume and will serve as a shorthand for the flexural rigidity
at a particular wavelength, and is given by
Te2ET(l, R, E , Te , ) = 2eR2 ) (lR212(13 (l + 1)3

4l2 (l + 1) 2 + 4l (l + 1)) + l (l + 1)
l (l + 1)

(1

2

)

(43)

2

where we have replaced each ∇ with −l(l + 1) (Turcotte et al., 1981; Willemann and Turcotte, 1982). Here again, R is the body's outer radius.
For finite rigidity (ζ < ∞), and if both ht and db are specified, then combining Eqs. (40), (41), and (42), it can be shown that
d

hb =

b
gt 1 + ht c gt

c

ht

gb

1+

.
(44)

gb

Similarly, if ζ < ∞, and both ht and dt are specified, then it can be shown that

hb =

gt

c

ht

gb

1+

c gt

1

dt
ht

.

(45)

This expression corresponds to Eq. (S23) of Hemingway et al. (2013) except that here we allow for the radial variation in gravity and we do notassume equal masses in columns of equal width (see Hemingway and Matsuyama, 2017).
For non-zero rigidity, and if the loading is entirely at the top of the shell, then db = 0 and Eq. (44) reduces to

hb =

gt

c

ht

gb

(Ct )

(46)

or, when the loading is entirely at the base of the shell, and dt = 0, Eq. (45) reduces to

hb =

gt 1
.
gb Cb

c

ht

(47)

In these expressions, Ct and Cb are compensation factors corresponding to the cases of top and bottom loading, respectively, and are given by

Ct =

1
1+

(48)

gb

and

Cb =

1
1+

.
(49)

c gt

In both cases, decreasing rigidity (Te → 0, ζ → 0) leads the compensation factor to approach unity. In the case of zero rigidity, both Eqs. (44) and(45) reduce to Eq. (12). For non-zero rigidity, the compensation factors mean different things for the cases of bottom and top loading. An infinitely
rigid shell (Te → ∞, ζ → ∞) leads to zero compensation (C = 0) in either case, but that means zero root thickness in the case of top loading or zerosurface relief in the case of bottom loading. For a finite elastic thickness, 0 < Te < ∞, in general Ct ≠ Cb. When gt/gb is assumed to be unity, our Eq.
(46) is identical to Eq. (21) of Turcotte et al. (1981), except that our ζ incorporates the term that was missing from Eq. (1) of Turcotte et al. (1981)(Willemann and Turcotte, 1982).
Appendix B. Supplementary data
Supplementary data to this article can be found online at https://doi.org/10.1016/j.icarus.2019.03.011.

129

Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal

References

Kraus, H., 1967. Thin Elastic Shells. Wiley, New York.
Kvorka, J., Čadek, O., Tobie, G., Choblet, G., 2018. Does Titan's long-wavelength topography contain information about subsurface ocean dynamics? Icarus 310, 149–164.
https://doi.org/10.1016/j.icarus.2017.12.010.
Lainey, V., Karatekin, Ö., Desmars, J., Charnoz, S., Arlot, J.-E., Emelyanov, N., Le PoncinLafitte, C., Mathis, S., Remus, F., Tobie, G., Zahn, J.-P., 2012. Strong tidal dissipationin Saturn and constraints on Enceladus' thermal state from astrometry. Astrophys. J.
https://doi.org/10.1088/0004-637X/752/1/14. http://stacks.iop.org/0004-637X/
752/i=1/a=14?key=crossref.d68a903d62213eb4f281cc2f4669349f.
Lainey, V., Jacobson, R.A., Tajeddine, R., Cooper, N.J., Murray, C., Robert, V., Tobie, G.,
Guillot, T., Mathis, S., Remus, F., Desmars, J., Arlot, J.-E., De Cuyper, J.-P., Dehant,
V., Pascu, D., Thuillot, W., Poncin-Lafitte, C.L., Zahn, J.-P., 2017. New constraints on
Saturn's interior from Cassini astrometric data. Icarus 281, 286–296. http://linkinghub.elsevier.com/retrieve/pii/S0019103516304183 https://doi.org/10.
1016/j.icarus.2016.07.014.
Lambeck, K., 1988. Geophysical Geodesy: The Slow Deformations of the Earth. Oxford
University Press, New York.
Lemasquerier, D., Grannan, A.M., Vidal, J., Cébron, D., Favier, B., Le Bars, M., Aurnou,
J.M., 2017. Libration-driven flows in ellipsoidal shells. Journal of GeophysicalResearch: Planets 122, 1926–1950. URL: http://doi.wiley.com/10.1002/
2017JE005340. https://doi.org/10.1002/2017JE005340.
Manga, M., Wang, C.Y., 2007. Pressurized oceans and the eruption of liquid water on
Europa and Enceladus. Geophys. Res. Lett. 1–5. https://doi.org/10.1029/
2007GL029297.
Matsuyama, I., Beuthe, M., Hay, H.C., Nimmo, F., Kamata, S., 2018. Ocean tidal heatingin icy satellites with solid shells. Icarus 312, 208–230. https://doi.org/10.1016/j.
icarus.2018.04.013. doi:10.1016/j.icarus.2018.04.013. arXiv:1804.07727.
McCarthy, C., Cooper, R.F., 2016. Tidal dissipation in creeping ice and the thermal
evolution of Europa. Earth Planet. Sci. Lett. 443, 185–194.
McKinnon, W.B., 2013. The shape of Enceladus as explained by an irregular core: implications for gravity, libration, and survival of its subsurface ocean. Journal ofGeophysical Research: Planets 118https://doi.org/10.1002/jgre.20122. http://doi.
wiley.com/10.1002/jgre.20122.
McKinnon, W.B., 2015. Effect of Enceladus's rapid synchronous spin on interpretation of
Cassini gravity. Geophys. Res. Lett. https://doi.org/10.1002/2015GL063384.
Meyer, J., Wisdom, J., 2007. Tidal heating in Enceladus. Icarus 188, 535–539. http://linkinghub.elsevier.com/retrieve/pii/S0019103507001054 https://doi.org/10.
1016/j.icarus.2007.03.001.
Meyer, J., Wisdom, J., 2008. Tidal evolution of Mimas, Enceladus, and Dione. Icarus 193,
213–223. https://doi.org/10.1016/j.icarus.2007.09.008.
Murray, C.D., Dermott, S., 1999. Solar System Dynamics. Cambridge University Press.
Nadeau, A., Mcgehee, R., 2017. A simple formula for a planet's mean annual insolation bylatitude. Icarus 291, 46–50. https://doi.org/10.1016/j.icarus.2017.01.040.
Nakajima, M., Ingersoll, A.P., 2016. Controlled boiling on Enceladus. 1. Model of thevapor-driven jets. Icarus 272, 309–318. https://doi.org/10.1016/j.icarus.2016.02.
027.
Nimmo, F., Thomas, P.C., Pappalardo, R.T., Moore, W.B., 2007. The global shape ofEuropa: constraints on lateral shell thickness variations. Icarus 191, 183–192.
https://doi.org/10.1016/j.icarus.2007.04.021.
Nimmo, F., Bills, B.G., Thomas, P.C., 2011. Geophysical implications of the long-wavelength topography of the Saturnian satellites. J. Geophys. Res. E11001. http://doi.
wiley.com/10.1029/2011JE003835 https://doi.org/10.1029/2011JE003835.
Nimmo, F., Porco, C., Mitchell, C., 2014. Tidally modulated eruptions on Enceladus:
Cassini ISS observations and models. Astron. J. 46. http://stacks.iop.org/1538-3881/
148/i=3/a=46?key=crossref.aa4e6bbc50ecc4ba5aec7e3316ec9995 https://doi.
org/10.1088/0004-6256/148/3/46.
Nimmo, F., Barr, A. C., Běhounková, M., & Mckinnon, W. B. (2018). The thermal andorbital evolution of Enceladus: observational constraints and models. In P. M. Schenk,
R. N. Clark, C. J. A. Howett, & A. J. Verbiscer (Eds.), Enceladus and the Icy Moons ofSaturn (pp. 79–94). Tucson, Arizona: University of Arizona Press. https://websites.
pmc.ucsc.edu/ fnimmo/website/Enceladus_in_press.pdf.
Ojakangas, G.W., Stevenson, D.J., 1989. Thermal state of an ice shell on Europa. Icarus
81, 220–241. https://doi.org/10.1016/0019-1035(89)90053-5.
Peale, S., Cassen, P., 1978. Contribution of tidal dissipation to lunar thermal history.
Icarus 36, 245–269. http://www.sciencedirect.com/science/article/pii/
0019103578901094 https://doi.org/10.1016/0019-1035(78)90109-4.
Petrenko, V.F., Whitworth, R.W., 1999. Physics of Ice. Oxford University Press.
Porco, C.C., Helfenstein, P., Thomas, P.C., Ingersoll, A.P., Wisdom, J., West, R., Neukum,G., Denk, T., Wagner, R., Roatsch, T., Kieffer, S., Turtle, E., McEwen, A., Johnson,
T.V., Rathbun, J., Veverka, J., Wilson, D., Perry, J., Spitale, J., Brahic, A., Burns, J.A.,
Delgenio, A.D., Dones, L., Murray, C.D., Squyres, S., 2006. Cassini observes the active
south pole of Enceladus. Science 1393–1401. (New York, N.Y.). http://www.ncbi.
nlm.nih.gov/pubmed/16527964. https://doi.org/10.1126/science.1123013.
Postberg, F., Schmidt, J., Hillier, J., Kempf, S., Srama, R., 2011. A salt-water reservoir as
the source of a compositionally stratified plume on Enceladus. Nature 474, 620–622.
http://www.nature.com/doifinder/10.1038/nature10175 https://doi.org/10.1038/
nature10175.
Richard, A., Rambaux, N., Charnay, B., 2014. Librational response of a deformed 3-layer
Titan perturbed by non-Keplerian orbit and atmospheric couplings. Planetary andSpace Science 94, 22–34. URL: https://doi.org/10.1016/j.pss.2014.02.006(doi:10.1016/j.pss.2014.02.006. arXiv:arXiv:1402.2493v1).
Roberts, J.H., 2015. The fluffy core of Enceladus. Icarus 258, 54–66. http://linkinghub.
elsevier.com/retrieve/pii/S001910351500250X https://doi.org/10.1016/j.icarus.
2015.05.033.
Roberts, J.H., Nimmo, F., 2008. Tidal heating and the long-term stability of a subsurface
ocean on Enceladus. Icarus 194, 675–689. http://linkinghub.elsevier.com/retrieve/Běhounková, M., Souček, O., Hron, J., Čadek, O., 2017. Plume activity and tidal deformation on Enceladus influenced by faults and variable ice shell thickness.
Astrobiology. https://doi.org/10.1089/ast.2016.1629.(ast.2016.1629. URL: http://
online.liebertpub.com/doi/10.1089/ast.2016.1629).
Besserer, J., Nimmo, F., Roberts, J.H., Pappalardo, R.T., 2013. Convection-driven compaction as a possible origin of Enceladus's long wavelength topography. Journal ofGeophysical Research E: Planets 118, 908–915. https://doi.org/10.1002/jgre.20079.
Beuthe, M., 2018. Enceladus's crust as a non-uniform thin shell: I tidal deformations.
Icarus 302, 145–174. https://doi.org/10.1016/j.icarus.2017.11.009. doi:10.1016/j.
icarus.2017.11.009. arXiv:1711.08236.
Beuthe, M., Rivoldini, A., & Trinh, A. (2016). Enceladus' and Dione's floating ice shells
supported by minimum stress isostasy. Geophys. Res. Lett., URL: http://doi.wiley.
com/10.1002/2016GL070650. doi:https://doi.org/10.1002/2016GL070650.
Bland, M.T., Beyer, R.A., Showman, A.P., 2007. Unstable extension of Enceladus' lithosphere. Icarus 192, 92–105. https://doi.org/10.1016/j.icarus.2007.06.011.
Čadek, O., Tobie, G., Hoolst, T.V., Mass, M., Choblet, Gael, Lefevre, A., Mitri, G., Baland,
R.-m., Běhounková, M., Bourgeois, O., Trinh, A., 2016. Enceladus's internal ocean
and ice shell constrained from Cassini gravity, shape and libration data. Geophys.
Res. Lett. https://doi.org/10.1002/2016GL068634.
Čadek, O., Běhounková, M., Tobie, G., Choblet, G., 2017. Viscoelastic relaxation of
Enceladus's ice shell. Icarus 291, 31–35. https://doi.org/10.1016/j.icarus.2017.03.
011.
Čadek, O., Souček, O., Běhounková, M., Choblet, G., Tobie, G., Hron, J., 2019. Long-term
stability of Enceladus' uneven ice shell. Icarus 319, 476–484.. URL: https://www.
sciencedirect.com/science/article/pii/S0019103518303762. https://doi.org/10.
1016/J.ICARUS.2018.10.003.
Choblet, G., Tobie, G., Sotin, C., Běhounková, M., Čadek, O., Postberg, F., Souček, O.,
2017. Powering prolonged hydrothermal activity inside Enceladus. Nature
Astronomy. https://doi.org/10.1038/s41550-017-0289-8. URL: http://www.nature.
com/articles/s41550-017-0289-8.
Crow-Willard, E.N., Pappalardo, R.T., 2015. Structural mapping of Enceladus and implications for formation of tectonized regions. Journal of Geophysical Research:
Planets 120, 928–950. https://doi.org/10.1002/2014JE004627.Received.
Darwin, G.H., 1899. The theory of the figure of the earth carried to the second order ofsmall quantities. Mon. Not. R. Astron. Soc. 60, 82–124.
Efroimsky, M., 2018. Dissipation in a tidally perturbed body librating in longitude. Icarus
306, 328–354. URL: https://doi.org/10.1016/j.icarus.2017.10.020 doi:10.1016/
j.icarus.2017.10.020. arXiv:1706.08999.
Fuller, J., Luan, J., Quataert, E., 2016. Resonance locking as the source of rapid tidal
migration in the Jupiter and Saturn moon systems. Mon. Not. R. Astron. Soc. 458,
3867–3879. http://arxiv.org/abs/1601.05804 https://doi.org/10.1093/mnras/
stw609 arXiv:1601.05804.
Giese, B., Wagner, R., Hussmann, H., Neukum, G., Perry, J., Helfenstein, P., Thomas, P.C.,
2008. Enceladus: an estimate of heat flux and lithospheric thickness from flexurally
supported topography. Geophys. Res. Lett. 1–5. https://doi.org/10.1029/
2008GL036149.
Goodman, J.C., Lenferink, E., 2012. Numerical simulations of marine hydrothermal
plumes for Europa and other icy worlds. Icarus 221, 970–983. https://www.sciencedirect.com/science/article/pii/S0019103512003466. https://doi.org/10.1016/J.
ICARUS.2012.08.027.
Goodman, J.C., Collins, G.C., Marshall, J., Pierrehumbert, R.T., 2004. Hydrothermal
plume dynamics on Europa: implications for chaos formation. J. Geophys. Res.
https://doi.org/10.1029/2003JE002073. URL: http://doi.wiley.com/10.1029/
2003JE002073.
Hemingway, D.J., Matsuyama, I., 2017. Isostatic equilibrium in spherical coordinates andimplications for crustal thickness on the Moon, Mars, Enceladus, and elsewhere.
Geophys. Res. Lett. 7695–7705. https://doi.org/10.1002/2017GL073334.
Hemingway, D.J., Nimmo, F., Iess, L., 2013. Enceladus' internal structure inferred fromanalysis of Cassini-derived gravity and topography. In: AGU Fall Meeting 2013
P53E-03.
Howett, C.J.A., Spencer, J.R., Pearl, J., Segura, M., 2011. High heat flow from Enceladus'
south polar region measured using 10-600 cm-1 Cassini/CIRS data. Journal ofGeophysical Research E: Planets 116, 1–15. https://doi.org/10.1029/2010JE003718.
Hsu, H.-W., Postberg, F., Sekine, Y., Shibuya, T., Kempf, S., Horányi, M., Juhász, A.,
Altobelli, N., Suzuki, K., Masaki, Y., Kuwatani, T., Tachibana, S., Sirono, S.-i.,Moragas-Klostermeyer, G., Srama, R., 2015. Ongoing hydrothermal activities within
Enceladus. Nature 519, 207–210. http://www.nature.com/doifinder/10.1038/
nature14262 https://doi.org/10.1038/nature14262.
Iess, L., Stevenson, D.J., Parisi, M., Hemingway, D.J., Jacobson, R.A., Lunine, J.I., Nimmo,
F., Armstrong, J.W., Asmar, S.W., Ducci, M., Tortora, P., 2014. The gravity field andinterior structure of Enceladus. Science (New York, N.Y.) 78–80. URL: http://
www.sciencemag.org/cgi/doi/10.1126/science.1250551. https://doi.org/10.1126/
science.1250551.
Ingersoll, A.P., Nakajima, M., 2016. Controlled boiling on Enceladus. 2. Model of theliquid-filled cracks. Icarus 272, 319–326. https://doi.org/10.1016/j.icarus.2015.12.
040.
Kamata, S., Nimmo, F., 2017. Interior thermal state of Enceladus inferred from the viscoelastic state of the ice shell. Icarus 284, 387–393. URL. https://doi.org/10.1016/j.
icarus.2016.11.034 (doi:10.1029/2011JE003835).
Kite, E.S., Rubin, A.M., 2016. Sustained eruptions on Enceladus explained by turbulent
dissipation in tiger stripes. Proc. Natl. Acad. Sci. 113, 201520507. URL: http://www.
pnas.org/lookup/doi/10.1073/pnas.1520507113. https://doi.org/10.1073/pnas.
1520507113.

130Icarus 332 (2019) 111–131

D.J. Hemingway and T. Mittal
pii/S0019103507005593 https://doi.org/10.1016/j.icarus.2007.11.010.
Rudolph, M.L., Manga, M., 2009. Fracture penetration in planetary ice shells. Icarus 199,
536–541. URL: https://doi.org/10.1016/j.icarus.2008.10.010 (doi:10.1016/
j.icarus.2008.10.010).
Schenk, P.M., McKinnon, W.B., 2009. One-hundred-km-scale basins on Enceladus: evidence for an active ice shell. Geophys. Res. Lett. https://doi.org/10.1029/
2009GL039916. URL: http://www.agu.org/pubs/crossref/2009/2009GL039916.
shtml.
Shoji, D., Hussmann, H., Kurita, K., Sohl, F., 2013. Ice rheology and tidal heating of
Enceladus. Icarus 226, 10–19. doi:10.1016/j.icarus.2013.05.004. https://doi.org/
10.1016/j.icarus.2013.05.004.
Soderlund, K.M., Schmidt, B.E., Wicht, J., Blankenship, D.D., 2014. Ocean-driven heatingof Europa's icy shell at low latitudes. Nat. Geosci. 16–19. https://doi.org/10.1038/
ngeo2021.
Soucek, O., Hron, J., Behounkova, M., Cadek, O., 2016. Effect of the tiger stripes on the
deformation of Saturn's moon Enceladus. Geophys. Res. Lett. 43, 7417–7423. https://doi.org/10.1002/2016GL069415.
Spencer, J.R., Howett, C.J.A., Verbiscer, A., Hurford, T.A., Segura, M., Spencer, D.C.,
2013. Enceladus heat flow from high spatial resolution thermal emission observations. In: European Planetary Science Congress. 8. pp. 840.
Spohn, T., Schubert, G., 2003. Oceans in the icy Galilean satellites of Jupiter? Icarus 161,
456–467. https://doi.org/10.1016/S0019-1035(02)00048-9.
Tajeddine, R., Rambaux, N., Lainey, V., Charnoz, S., Richard, A., Rivoldini, A., Noyelles,
B., 2014. Constraints on Mimas' interior from Cassini ISS libration measurements.
Science 346, 322–324.
Tajeddine, R., Soderlund, K.M., Thomas, P.C., Helfenstein, P., Hedman, M.M., Burns, J.A.,
Schenk, P.M., 2017. True polar wander of Enceladus from topographic data. Icarus
295, 46–60. http://linkinghub.elsevier.com/retrieve/pii/S001910351630584X
https://doi.org/10.1016/j.icarus.2017.04.019.
Thomas, P., 2010. Sizes, shapes, and derived properties of the Saturnian satellites after
the Cassini nominal mission. Icarus 208, 395–401. http://linkinghub.elsevier.com/

retrieve/pii/S0019103510000448 https://doi.org/10.1016/j.icarus.2010.01.025.
Thomas, P., Tajeddine, R., Tiscareno, M., Burns, J., Joseph, J., Loredo, T., Helfenstein, P.,
Porco, C., 2016. Enceladus's measured physical libration requires a global subsurface
ocean. Icarus 264, 37–47. http://linkinghub.elsevier.com/retrieve/pii/
S0019103515003899 https://doi.org/10.1016/j.icarus.2015.08.037.
Tricarico, P., 2014. Multi-layer hydrostatic equilibrium of planets and synchronous
moons: theory and application to Ceres and to solar system moons. Astrophys. J. 99.
http://stacks.iop.org/0004-637X/782/i=2/a=99?key=crossref.
51a63b8cfe5c9568438c8375874b01f7 https://doi.org/10.1088/0004-637X/782/
2/99.
Turcotte, D.L., Willemann, R.J., Haxby, W.F., Norberry, J., 1981. Role of membranestresses in the support of planetary topography. J. Geophys. Res. 3951–3959.
Van Hoolst, T., Rambaux, N., Karatekin, Ö., Baland, R.M., 2009. The effect of gravitational and pressure torques on Titan's length-of-day variations. Icarus 200, 256–264.
https://doi.org/10.1016/j.icarus.2008.11.009.
Van Hoolst, T., Baland, R.M., Trinh, A., 2013. On the librations and tides of large icy
satellites. Icarus 226, 299–315. https://doi.org/10.1016/j.icarus.2013.05.036.
Van Hoolst, T., Baland, R.-M., Trinh, A., 2016. The diurnal libration and interior structureof Enceladus. Icarus 277, 311–318. http://linkinghub.elsevier.com/retrieve/pii/
S0019103516302020 https://doi.org/10.1016/j.icarus.2016.05.025.
Vance, S., Goodman, J., 2009. Oceanography of an Ice-Covered Moon. https://www.
researchgate.net/publication/241525803.
Wieczorek, M., 2015. Gravity and topography of the terrestrial planets. In: Treatise on
Geophysics, 2nd ed. Elsevier B.V, pp. 153–193. https://doi.org/10.1016/B978-0444-53802-4.00169-X.
Wieczorek, M.A., Phillips, R.J., 1998. Potential anomalies on a sphere: applications to thethickness of the lunar crust. J. Geophys. Res. 1715–1724.
Willemann, R.J., Turcotte, D.L., 1982. The role of lithospheric stress in the support of theTharsis rise. J. Geophys. Res. 9793–9801.
Wilson, A., Kerswell, R.R., 2018. Can libration maintain Enceladus's ocean? Earth Planet.
Sci. Lett. 41–46. https://doi.org/10.1016/j.epsl.2018.08.012.

131