Hemingway D., Iess L., Tajeddine R., and Tobie G. (2018) The interior of Enceladus. In Enceladus and the Icy Moons of Saturn(P. M. Schenk et al., eds.), pp. 57–77. Univ. of Arizona, Tucson, DOI: 10.2458/azu_uapress_9780816537075-ch004.

The Interior of Enceladus
Douglas Hemingway
University of California, Berkeley

Luciano Iess

Università La Sapienza, Rome

Radwan Tajeddine
Cornell University

Gabriel Tobie

Université de Nantes

Knowledge of internal structure places constraints on the origin, evolution, and present
behavior of Enceladus. Although the interior cannot be observed directly, much can be deduced from the magnitude of tidal and rotational asymmetries in the shape and gravity field
and from how these asymmetries affect the rotational dynamics. The data and models suggest
a differentiated interior consisting of a ~2400 kg m–3 rocky core, a 20–50-km mean thickness
global subsurface liquid water ocean, and a 20–30-km mean thickness icy shell. The presenceof an internal liquid water reservoir had been suggested by a series of observations early in the
Cassini mission and was ultimately confirmed by two crucial and independent analyses related
to (1) determination of Enceladus’ low-order gravity field based on Doppler tracking of the
Cassini spacecraft during dedicated flybys and (2) measurement of forced physical librations
based on control point tracking through seven years of Cassini images. Although Cassini-based
observations have yielded a good, basic picture of the interior of Enceladus, many unresolved
issues remain — in particular, relating to the structure and stability of the icy shell and related
implications for the energy budget and thermal evolution of Enceladus — and will need to beaddressed through further modeling and, eventually, additional observations.

1. INTRODUCTION

polar region (see chapters in this volume by Postberg et al.,Goldstein et al., and Spencer et al.). Are these eruptions
sourced from a potentially habitable subsurface liquid water
reservoir? If so, how deep is the reservoir beneath the icy
surface? How thick is the icy crust? What is the structure of
the crust and how does it support the surface topography?
How is the liquid water, which is not buoyant with respect
to the ice, transported to the surface? How does the interior respond to tidal stresses? What modulates the eruptive
activity at the south pole? How large is the internal sea orocean, and how does it keep from freezing solid? What is thecomposition and nature of the core, and what is its capacity
for contributing to the heating of the interior?
Although the internal structure cannot be observed directly, a variety of techniques are available for building a
picture of the interior that can help to address many of thesequestions. Seismic imaging has proven to be the definitivetechnique for constraining the structure of Earth’s interior,and seismic stations placed on the Moon during the Apollo
program have likewise helped to improve our understanding

1.1. Motivation
The internal structure of a planetary body provides essential clues about its origin and evolution. Based on its
internal density structure, what is the likely composition of
the materials from which the body formed (see the chapter
in this volume by McKinnon et al.)? Did the interior become
sufficiently warm to enable differentiation, with heavier
components sinking toward the center? How quickly or
slowly did the body form? What does the structure of thenear surface and deeper interior tell us about the nature of
the topography and the distribution of geologic provinces
across the surface?
The interior of Enceladus is of particular interest becauseof the complex tectonic structures that cover much of its
unusually youthful surface (see chapter in this volume byPatterson et al.) and, most strikingly, because of the ongoing eruptions of water vapor and ice grains from its south
57

58

Enceladus and the Icy Moons of Saturnof the Moon’s interior. However, although seismic stations
may soon be deployed on Mars and other worlds, the studyof planetary interiors aside from Earth and the Moon has so
far had to rely on other methods. In particular, as we will
discuss in detail in this chapter, studies of shape, gravitational fields, and rotational dynamics can place constraints
on planetary interiors, and have enabled significant advancements in our emerging understanding of Enceladus.
1.2. Key Observations and Basic Picture
Prior to the Cassini spacecraft’s arrival in the Saturnsystem in 2004, our knowledge of Enceladus was limited
to what could be inferred from telescopic observations
and the data collected during the Saturn system flybys of
Pioneer 11 in 1979 and the two Voyager spacecraft in 1980
and 1981. However, even these early observations led toconsiderable interest in Enceladus due to its unusually highalbedo (Buratti and Veverka, 1984) and apparently youthfulicy surface (Smith et al., 1982; Clark et al., 1983; Squyres
et al., 1983), as well as its relationship to Saturn’s E ring
(Baum et al., 1981; Haff et al., 1983; Pang et al., 1984), allof which pointed to the incredible possibility of recent, oreven ongoing, geologic activity on this tiny moon. Given
the apparent history of geologic activity, and incorporating
mass estimates based on astrometric observations (Kozai,
1976; Jacobson, 2004), basic interior models were developed for a differentiated Enceladus (Zharkov et al., 1985),
but were limited by the large uncertainties regarding its
size and shape.
Beginning in 2005, the Cassini mission, which included
numerous flybys of Enceladus, led to a series of remarkable
observations that reinforced Enceladus’ status as one of themost compelling exploration targets in the solar system, andprovided clues about the nature of its interior. In particular,
several of Cassini’s instruments detected water vapor andorganic molecules in the vicinity of Enceladus’ south pole
(Brown et al., 2006; Dougherty et al., 2006; Hansen et al.,2006; Spahn et al., 2006; Spencer et al., 2006; Waite et al.,2006), and spectacular back-lit images revealed the source of
this material to be a series of ongoing eruptions concentrated
along four major fissures spanning the highly tectonized southpolar terrain (SPT) (Porco et al., 2006). Subsequent observations identified high heat flow emanating from the southpolar region and especially at the sites of the most prominent
eruptions (Spitale and Porco, 2007; Howett et al., 2011).
The ongoing activity led to the suggestion that a subsurface
sea might be feeding the eruptions (Collins and Goodman,2007) and producing a topographic depression at the southpole, helping to explain the anomalous shape (Thomas et al.,2007). The finding that the erupted ice grains contained salts(Postberg et al., 2009, 2011) further supported the notionof an internal liquid water reservoir. On the other hand, thelimited presumed available tidal heating energy (Meyer andWisdom, 2007), together with the rapid heat loss, suggested
that an internal sea should have frozen solid in just tens of
millions of years (Roberts and Nimmo, 2008).

A major advancement in our knowledge of the interiorcame with Doppler tracking of the Cassini spacecraft during close encounters with Enceladus (Rappaport et al.,2007), and especially with the determination of Enceladus’
quadrupole gravity field and hemispherical asymmetry, once
enough flyby results had been accumulated (Iess et al.,2014). In combination with improved models of the topography (Thomas et al., 2007; Nimmo et al., 2011), the gravity
data provided the first opportunity to effectively probe theinterior, yielding an estimate of the degree of differentiation,and identifying mass anomalies suggestive of a (possibly
global) subsurface liquid ocean. Subsequent measurement of
the forced physical librations (Thomas et al., 2016) independently confirmed that the icy crust must be fully decoupledfrom the deep interior, indicating that the internal ocean mustindeed be global. Although detailed interpretation of theseobservations is complicated and remains an active area of
research (e.g., Čadek et al., 2016, 2017; van Hoolst et al.,2016; Beuthe et al., 2016; Hemingway and Mittal, 2017),
the general conclusion is that Enceladus is differentiated,with an icy shell covering a global subsurface liquid water
ocean, overlying a low-density rocky core (Fig. 1).
In this chapter, we examine the relevant theory anddiscuss the constraints that can be placed on the internal
structure of Enceladus based mainly on the gravity andlibration observations (sections 2 and 3, respectively). We
also discuss the implications for the thermal state of Enceladus and, in particular, the structure and dynamics of its
icy shell (section 4). Finally, we discuss the remaining open

Fig. 1. Approximate interior structure of Enceladus (toscale). Its 252-km radius comprises a ~190-km-radius corewith a density of ~2400 kg m–3, surrounded by layers ofliquid and solid H2O making up the remaining ~60 km. The
ice shell and ocean layer thicknesses vary laterally, with thethinnest part of the ice shell (and thickest part of the ocean)being centered on the south pole (see Table 3).

Hemingway et al.: The Interior of Enceladus   59questions and how future observations and modeling might
help to address them (section 5).
2. IMPLICATIONS FROM GRAVITY
Gravity is a powerful tool for studying planetary interiors
because a body’s gravitational field is ultimately a functionof its internal mass distribution. While the problem of inverting gravity for interior structure suffers from inherent nonuniqueness, with a few reasonable assumptions, it is nevertheless possible to draw useful conclusions. For example,
it is often reasonable to assume that large planetary bodies
have interiors that are weak on long timescales, meaningthat they tend to relax to near spherical symmetry, and that
their internal densities do not generally increase with radial
position — a situation that would be gravitationally unstable.
A straightforward approach to interior modeling is therefore to treat the body as a series of concentric sphericalshells of different densities, with the requirement that thedensities increase monotonically inward. The mean density,−, of such a body is given by
r
ρ=

R 
∆ ρi  i 
 R

∑
i

3

(1)

where R is the body’s full radius, Ri is the outer radius of
the ith layer, and Dri is the density contrast between layer iand the layer above it. This formulation implicitly assumes
uniform density within each layer, which is appropriate
only when internal pressures are small enough not to causesignificant compression. For larger bodies, where compression may be important, equations of state are required tomodel the radial dependence of density on temperature andpressure. We proceed here assuming that compression is notimportant in the interior of Enceladus.
The body’s normalized mean moment of inertia (or moment of inertia factor) is given by
I
MR

=
2

2
5

∆ρi  R i ∑ ρ  R 

5

(2)

i

Although there are generally more than two unknowns,
necessitating additional assumptions, these two equations
provide fundamental constraints on the internal structure.
The mean density can be determined directly from the mass
and radius, and immediately restricts the range of possible
bulk compositions. But the moment of inertia, which reflects
density stratification, must be determined via other means.
The three principal moments of inertia can be defined as)∫ (B = ∫ ρ ( x + z ) dV
C = ∫ ρ ( x + y ) dV
A = ρ y 2 + z 2 dV
2

2

2

2

where r is the density of a volume element dV, at position x,y, z within the body. Here, the principal frame is a coordinate
frame centered on the body with the z-axis coincident withthe spin pole, the x-axis pointing toward the prime meridian, and the y-axis completing the righthanded system (forsynchronous satellites in circular orbits, the x-axis points
along the long axis, toward the parent body). Differencesin these principal moments of inertia are important because
they can, in some cases, be related to observable rotational
dynamical effects such as spin pole precession, obliquity,
nutation, and librations (e.g., van Hoolst, 2015) — thesubject of section 3 of this chapter; or to asymmetries in
the gravitational field — the subject of the rest of this section. In the limit of perfect spherical symmetry, a body’s
gravitational field would be equivalent to that of a point
mass located at its center and, apart from the mean density,
little could be concluded about the interior. Real planetary
bodies, however, exhibit asymmetries that effectively carryinformation about their internal structures.
The case of Enceladus is interesting and instructive forseveral reasons. Because of the short period of its orbit
around Saturn (1.37 days), the asymmetries in its shape andgravity (primarily due to tidal and rotational deformation)
are large, permitting unprecedented relative precision in
their measurement. This precision presents both challenges(it forces us to abandon the usual simplifying assumptionof hydrostatic equilibrium) and opportunities (it allows us
to estimate the thickness of the icy crust via analysis of the
non-hydrostatic gravity and topography), forcing the development of new methods for interior modeling from shapeand gravity (Hemingway et al., 2013b; Iess et al., 2014;
McKinnon, 2015; Čadek et al., 2016; Beuthe et al., 2016).
Finally, the special circumstance of the substantial northsouth polar asymmetry associated with the SPT provides
an opportunity for independent analysis of long-wavelength
gravity and topography.
In section 2.1, we discuss how rotational and tidal effects
lead to characteristic asymmetries in a body’s shape and,
consequently, its gravitational field, and how the magnitudeof those asymmetries is related to the internal mass distribution. To do this, we introduce several concepts, including
gravitational potential and the ways it is affected by tidal androtational forces; hydrostatic equilibrium figure theory; andLove numbers and their connection to the principal moments
of inertia and, ultimately, the density stratification within the
body. In section 2.2, we discuss how the gravitational fieldof Enceladus has been measured based on radio trackingof the Cassini spacecraft during three close flybys. Finally,in section 2.3, we discuss the interpretation of those measurements, including some of the subtleties and challenges
associated with interior modeling.
2.1. Theory(3)

2.1.1. Hydrostatic equilibrium. For planetary bodies
that are sufficiently large, the combination of high internalpressures and low internal viscosities (especially true when60   Enceladus and the Icy Moons of Saturn
there is appreciable internal heating) results in relaxation toa figure that approaches the expectation for a strengthless
fluid body, in which the inward acceleration due to gravity is everywhere balanced by the gradient in fluid pressure — a condition referred to as hydrostatic equilibriumand described by the equation dp = –rgdr, where p, r, g,
and r are the pressure, local density, local gravity, and radial
position, respectively.
Under the influence of self-gravitation alone, the hydrostatic equilibrium figure would be a sphere. Because of theirrotation, however, planetary bodies also experience centrifugal flattening. In addition, satellites that are in synchronous
rotation with their parent bodies (i.e., tidally locked) also
experience permanent elongation along the static tidal axis.
This applies to most of the large natural satellites in the solar
system because of the short timescale associated with tidal
locking (e.g., Gladman et al., 1996; Murray and Dermott,1999). The tidal and rotational deformation results in an
asymmetry in the equilibrium figure — and consequently thegravitational field — of synchronous satellites like Enceladus.
2.1.2. Gravitational potential. A body’s Newtonian gravita‑
tional potential, U, resolved at an arbitrary position r, is afunction of its internal mass distribution, and is given by
U (r )

ρ (r′ )–GdV ′V r – r′∫(4)

where r(r′) is the density at position r′, where G is theuniversal gravitational constant, and where the integral is
performed over the body’s entire volume, V. The potential
represents the work per unit mass done by the gravitational
field and gives rise to the gravitational acceleration g = –∇U.
Everywhere outside the body, the gravitational potential
satisfies Laplace’s equation, ∇2U(r) = 0, and can be expressed
as a linear combination of spherical harmonic functions as
U ( r, θ, φ ) = –

GM
r

∞

l

∑∑

l = 0 m = –ll

 R ref  r  Clm Ylm ( θ, φ )(5)where M is the total mass of the body, r is the radius atwhich the potential is to be expressed, q is colatitude, f
is longitude, Clm are the degree-l and order-m sphericalharmonic expansion coefficients representing the dimensionless gravitational potential at the reference radius Rref,
and where Ylm are the spherical harmonic functions — thenatural set of orthogonal basis functions on a sphere (e.g.,
Wieczorek, 2015).
2.1.3. Tidal/rotational disturbing potential. The tidal
and rotational forces that cause planetary bodies to deviate
from spherical symmetry can be described in terms of disturbances in the potential field. For instance, the centrifugal
acceleration at a point on the surface of a spherical bodywith radius R and spin rate w is w2x, where x = R sin q x̂is a vector pointed outward from, and perpendicular to, theaxis of rotation and reaching the surface at colatitude q. This
acceleration can be expressed as the negative gradient of acentrifugal disturbing potential, –∇Vcf, where

V cf– 12 ω 2 R 2 sin 2 θ(6)For convenient comparison with equation (5), we canrewrite the centrifugal disturbing potential in terms of thedegree-2 Legendre polynomial, P2(cos q) = 12 (3 cos2 q–1),
as Vcf(q,f) = 13w2R2(P2(cos q)–1). Since the last term is aconstant, its gradient is zero, meaning that it gives rise to
no deforming forces. As such, many authors do not include
it, and we will likewise disregard it from here on. Finally,
we can write the centrifugal disturbing potential (excluding
the constant term) in terms of the degree-2 zonal sphericalharmonic function, Y20(q,f), as
V cf ( θ, φ ) = 13 ω 2 R 2 Y20 ( θ, φ )(7)

In addition to the effects of rotation, satellites experience a tidal disturbing potential due to the spatially varying
difference between their gravitational acceleration toward
the parent body and the outward centrifugal acceleration
associated with their orbital motion. To second order, it canbe shown (e.g., Murray and Dermott, 1999, p. 133) that the
net tidal disturbing potential is
V tid ( ψ ) = –

GmR 2
a3

P2 ( cos ψ )(8)where m is the mass of the parent body, a is the distance
between the centers of the two bodies, and y is the angle
from the axis connecting the centers of the two bodies to anarbitrary point on the satellite’s surface. If the parent bodyis much more massive than the satellite (m ? M), and thesatellite’s orbit is synchronous with its mean rotation rate,then w2 ≈ Gm/a3. And if the satellite’s spin axis is nearlynormal to its orbital plane (i.e., its obliquity is close to zero),then we can rewrite the angle y in terms of colatitude andlongitude as cos y ≈ cos f sin q. Both of these conditions
are satisfied in the present case because Saturn is more thana million times as massive as Enceladus, whose obliquity
is <0.001 (Chen and Nimmo, 2011; Baland et al., 2016).
Incorporating these approximations, it can be shown that
V tid ( θ, φ ) = ω 2 R 2( 12 Y20 (θ, φ) – 14 Y22 (θ, φ))(9)Adding equations (7) and (9) yields the combined effectof the tidal and rotational disturbing potentials
V ( θ, φ ) = ω 2 R 2( 56 Y20 (θ, φ) – 14 Y22 (θ φ))(10)

2.1.4. Equilibrium figure theory. To the extent that the
body can be treated as a hydrostatic fluid, its equilibrium
figure will conform to an equipotential. The effect of the
disturbing potential, V(q,f), is that the surface of the initiallyspherically symmetric body will no longer be equipotential.
To first order, and before any deformation has taken place,the resulting change in elevation of the equipotential surface
is –V(q,f)/g (since g is the gradient of the gravitational potential). In response to the disturbing potential, the body’sHemingway et al.: The Interior of Enceladus   61figure will relax toward the new equipotential surface.
However, in so doing, the body’s mass distribution is altered,effectively causing a small additional disturbance in thegravitational potential. Accounting for this self-gravitation
effect, the hydrostatic equilibrium figure can be written
H (θ φ) = R –V ( θ, φ )
hf
g

(11)

where hf is the fluid Love number for the figure, a scalar quantity describing the magnitude of the body’s longtimescale (zero frequency) deformation in response to the
disturbing potential. An infinitely rigid body would have hf =0, whereas a perfectly homogeneous fluid body has hf = 52
(the value is larger than unity because of the aforementioned
self-gravitation); fluid bodies with some internal density
stratification have 0 < hf < 52 . The resulting asymmetries in
the gravitational potential are given byU ( θ, φ ) = V ( θ, φ ) k f(12)which defines the potential fluid Love number kf. For
perfectly fluid bodies, hf = 1 + kf. The fluid Love numbers
are of great importance for interior modeling because they
can be related to the body’s internal density distribution, as
discussed in the next section. Since the disturbing potential
in equation (10) is a degree-2 spherical harmonic function,
we are here concerned with the degree-2 fluid Love numbers, h2f and k2f.
Note that the fluid Love numbers should not be confusedwith the tidal Love numbers, which characterize the body’s
viscoelastic response in both shape (h2t) and gravitational
potential (k2t) due to eccentricity tides — variations in tidal
potential associated with the body’s changing proximity toits parent body. Ignoring these short-timescale variationsin shape and gravity, a synchronous satellite in hydrostatic
equilibrium thus has a mean figure described by
H ( θ, φ ) = R – Rqh 2f

(

5
Y
6 20

( θ, φ )– 14 Y22( θ, φ ))(13)

and a gravitational field whose asymmetries are described by
U ( θ, φ ) =

GM
qk 2f
R

( 56 Y20 (θ, φ) – 14 Y22 (θ, φ))(14)

where we have introduced the commonly used parameter
q=

ω2R
mR 3
=
g
Ma 3

(15)

representing the relationship between the inward gravitational acceleration and the outward centrifugal acceleration
at the body’s equator.
It is clear from equation (14) that, for a synchronous
satellite in hydrostatic equilibrium, the non-central part of
the gravitational field described by equation (5), and resolved
at the reference radius, has just two non-zero coefficientsJ 2 = – C20 = 56 qk 2f
C22 = 14 qk 2f(16)

where we have followed the convention of expressing thezonal gravity coefficient according to Jl = – Cl0. For hydrostatic synchronous satellites, the ratio of the gravity coefficients is thus J2 /C22 = 10/3.
Likewise, if the body’s shape is expressed in sphericalharmonics as
H ( θ, φ ) =∞l∑ ∑ H Y (θ, φ)
lm

lm

(17)

l = 0 m = –l

where Hlm are the degree-l and order-m spherical harmonic
expansion coefficients representing the shape (with H00 =R), then the hydrostatic equilibrium figure is asymmetric inprecisely the same way, with the coefficients
H 20 = – 56 qRh 2f
H 22 = 14 qRh 2f(18)

which, again, have the characteristic ratio H20/H22 = –10/3.
It is also common to describe the figure using the semiaxes(a > b > c) of an (approximately) equivalent triaxial ellipsoid,by evaluating equation (17) at the pole and at two points on
the equator (at longitudes 0 and p/2), yielding()
b = R (1 – 13 qh 2f )c = R (1 – 56 qh 2f )a = R 1 + 76 qh 2f(19)This leads to another commonly referenced characteristic
ratio for synchronous satellites in hydrostatic equilibrium
a–c
=4
b– c

(20)

The forgoing describes the condition of hydrostatic equilibrium in terms of both the expected equilibrium figure andthe resulting asymmetries in the gravitational field. However, whereas the first-order approximation used to obtainequation (11) is adequate in most cases, and sufficient fora basic interpretation of the gravity data (e.g., Iess et al.,2014), higher-order approximations become important for
more precise modeling of fast rotating bodies like Enceladus(Tricarico, 2014; McKinnon, 2015), as we will discuss insection 2.3.
2.1.5. Moments of inertia. In the previous section, we
showed that, in the case of a body in perfect hydrostatic
equilibrium, measuring either of the two degree-2 gravity
coefficients, J2 or C22, is sufficient to recover the fluid Love
number k2f. This is very useful because k2f can be related
to the body’s moment of inertia using the Radau-Darwin
equation (e.g., Darwin, 1899; Murray and Dermott, 1999)

62   Enceladus and the Icy Moons of Saturn
C
MR

=

2

2  2 4 – k 2f 
1–
3  5 1 + k 2f (21)

where C is the polar moment of inertia (i.e., about the spinpole, or c-axis).
The asymmetric nature of the tidally and rotationally deformed body means that its principal moments of inertia are
not equal (see equation (3)). For a hydrostatic synchronous
satellite, A < B < C; in the case of a non-synchronously
rotating body, with no static tidal bulge, A = B < C.
It can be shown (e.g., Hubbard, 1984, p. 79) that the
moment of inertia asymmetries are related to the degree-2
gravity coefficients according to
J2 =
C22 =

C – 12 ( A + B)
MR 2
B–A

(22)

4MR 2

from which it follows that the mean moment of inertia is
I

C

2
=
– J2
2
2
3
MR
MR

(23)

Having obtained the mean moment of inertia factor, weare finally in a position to constrain the internal density
structure described by equations (1) and (2).
2.1.6. Nonhydrostatic effects. An additional complication arises from the fact that solid planetary bodies, evenlarge ones like Earth, have exteriors that are cold and rigid
enough to support some non-hydrostatic topography, evenover long timescales. In general, the measured shape andgravity are therefore a reflection of a mostly hydrostaticbody, superimposed with some relatively smaller nonhydrostatic topography (and corresponding non-hydrostatic
gravity). The challenge is that it is not necessarily obvious
how to separate the hydrostatic and non-hydrostatic partsof these signals.
In previous work [e.g., Europa (Anderson et al., 1998)],
it was not possible to obtain independent constraints on J2
and C22 or to measure H20 and H22 with sufficient precision
to confirm whether or not they conform to the hydrostatic
expectation. Instead, the condition of hydrostatic equilibrium
was assumed, allowing direct calculation of k2f from C22 via
equation (16), and yielding an estimate of the moment ofinertia via equation (21). The Cassini mission, however, hasenabled more precise and nearly independent measurementsof each of these four quantities for the saturnian satellites
Titan, Enceladus, Dione, and Rhea, and has revealed that thecondition of hydrostatic equilibrium is generally not satisfied. For example, Titan’s figure exhibits considerable excess
flattening, with –H20/H22 = 4.9 ± 0.1 (Zebker et al., 2012;
Hemingway et al., 2013a) in spite of its nearly hydrostatic
gravity, for which J2/C22 = 3.32 ± 0.02 (Iess et al., 2010,2012). Similarly, Enceladus has considerable excess flatten-

ing, with –H20/H22 = 4.2 ± 0.2 (Nimmo et al., 2011; Thomas
et al., 2016; Tajeddine et al., 2017), while its gravity field
is only modestly non-hydrostatic, with J2/C22 = 3.51 ± 0.05
(Iess et al., 2014). For Dione, both the figure and gravity
are considerably non-hydrostatic, with –H20/H22 = 5.2 ± 0.6
and J2/C22 = 4.00 ± 0.06 (Thomas et al., 2007; Nimmo et al.,2011; Hemingway et al., 2016). While Rhea’s shape is not
as well determined, with –H20/H22 = 3.4 ± 1.0, its gravity
field is substantially non-hydrostatic, with J2/C22 = 3.91 ±0.10 (Tortora et al., 2016).
This deviation from hydrostatic equilibrium complicatesthe interpretation, as we discuss below, but the combinationof both shape and gravity observations allows some of theambiguities to be resolved, and additionally provides information about how such unrelaxed topography is supported
in the relatively stiff exterior. For example, the large nonhydrostatic topography and small non-hydrostatic gravity ofTitan and Enceladus are suggestive of isostatic compensation(see section 2.3).
2.2. Observations
2.2.1. Gravity determination. Cassini carried out gravity
measurements of Enceladus during three flybys on April 28,
2010; November 30, 2010; and May 2, 2012 (labeled as E9,
E12, and E19). Flyby altitude and latitude were selected toenhance the signature of a hemispherical asymmetry of thegravity field, characterized mainly by the harmonic coefficient J3. E9 and E19 occurred over the south polar region(latitude 89°S and 72°S), at closest approach altitudes of
100 and 70 km, respectively. E12 was over the north polar
region (latitude 62°N), at an altitude of just 48 km. The
small closest approach distances were required by the needto maximize the spacecraft accelerations due to Enceladus’
gravity field.
The estimation of the mass and gravity field was obtained solely from measurements of the spacecraft range
rate (Doppler). Range rate is approximately equal to twice
the line-of-sight projection of the spacecraft velocity withrespect to a ground antenna. This observable quantity is astandard product of the radio tracking system and is crucial
for the accurate navigation of the spacecraft. Although otherradiometric data are used for orbit reconstruction (such asrange), range rate is by far the most valuable for geodesy
applications. The range rate of Cassini is measured at a
ground antenna of NASA’s Deep Space Network (DSN)
after establishing a coherent, two-way, microwave link at
X-band (7.2–8.4 GHz). In this configuration, the groundantenna transmits a radio signal (a monochromatic carrier,possibly modulated for telecommands and ranging) generated by a highly stable frequency reference. The signal isthen received onboard by means of a transponder, thencoherently retransmitted back to Earth. The frequency ofthe beat tone between the incoming and the outgoing signal
provides the Doppler shift and the (two-way) range rate.
The Cassini X-band radio system delivers measurement
accuracies up to 0.01 mm s–1 at 60-s integration times un-

Hemingway et al.: The Interior of Enceladus   63der favorable conditions (although 2–4 times higher noise
is more frequently encountered). On timescales relevant
to Cassini gravity investigations, the noise is dominatedby path delay variations due to interplanetary plasma andtropospheric water vapor. During the three gravity flybys ofEnceladus, the (two-way) range rate noise was 0.017 mm s–1(E9), 0.027 mm s–1 (E12), and 0.036 mm s–1 (E19) at 60 s
integration time.
The estimation of Enceladus’s gravity field with Cassini
poses several challenges, all traceable to three factors. First,
the small number of flybys strongly limits the samplingof the gravity field and does not permit breaking the correlations between some estimated parameters. Second, thespacecraft acceleration due to the moon’s gravity is quite
small (surface gravity is only ~0.11 m s–2), and so is thechange in the spacecraft range rate measured by the groundantenna. Third, the gravitational interaction time of thespacecraft with the small moon is short.
The velocity variation induced by the monopole gravity, and the zonal harmonics of degree 2 and 3 can becomputed from elementary arguments. Neglecting factors
on the order of unity, we simply multiply the acceleration
due to each harmonic by the interaction time r/V (r being
the Cassini-Enceladus distance and V the nearly constantrelative velocity) and obtaindV( 0) ≈GMrV2

GM  R dV( 2) ≈  J2rV  r dV( 3) ≈(24)

3

GM  R   J3rV  r 

Here, GM is the gravitational parameter of Enceladus and
R is the reference radius (252.1 km). Note that, for close flybys (r–R = R), the velocity variation scales with the productrR2, where r is the mean density. Using the values reported
by Iess et al. (2014), the velocity changes due to GM, J2,
and J3 are, respectively, about 3–4 m s–1, 1.1–1.4 cm s–1, and
0.2–0.3 mm s–1. Although all these values are well above
the noise threshold of the Cassini radio system, the various
contributions are not fully separable with only three flybys.
Gravity and orbits of the spacecraft and Enceladus wererecovered by fitting the range rate data using orbit determination software developed at JPL for deep space navigation.
These tools (ODP and the recently developed MONTE) generate computed observables using a fully relativistic modelof the spacecraft dynamics and radio signal propagation.
The residuals (observed range rate minus computed rangerate) are minimized through an iterative process where a setof free parameters is adjusted by linearizing the observation
equations. The output of the process is a set of estimatedparameters and their covariances. The estimated parameters
included the full satellite degree-2 harmonic coefficients, J3,
the position and velocity of Cassini at each flyby, and corrections to the mass and the orbital elements of Enceladus. In

the nominal solution, the tesseral and sectorial componentsof the degree-3 field were assumed to be zero. This set ofparameters was the minimum set able to fit the data. As a testof the stability of the model, an alternate solution, which also
estimated the remaining degree-3 terms, was obtained. The
result is statistically identical to the nominal result, although
the uncertainties increase. All the additional degree-3 terms
are compatible with zero within 2s or less, and thereby
add no new information (Iess et al., 2014, section S2.1).
The non-gravitational accelerations due to the anisotropic
thermal emission from Cassini’s radioisotope thermoelectric
generators and the solar radiation pressure were accounted
for using previous estimates obtained during the saturnian
tour. The associated uncertainties were considered in theconstruction of the covariance matrix.
A good orbital fit would not be possible without also
including the small but non-negligible aerodynamic dragexperienced by Cassini when it flew through the plumes
in the southern polar flybys E9 and E19. The dynamical
effect was modeled as a nearly impulsive, vectorial acceleration at closest approach, to be estimated together
with the other parameters. An equivalent method relies onmodels of the neutral particle density and estimates the
aerodynamic coefficient of the spacecraft. Both methods
lead to statistically identical solutions for the gravity field.
The estimated, aerodynamic DV is almost parallel to thespacecraft velocity V (as expected for a drag force), with amagnitude of 0.25 mm s–1 for E9 and 0.26 mm s–1 for E19.
The DV experienced by Cassini is comparable to the effectof the J3 harmonic.
As expected, the gravity field is dominated by the J2
and C22 harmonics, associated with the rotational and tidal
deformation of Enceladus (Table 1). The estimate of J2 andC22 is quite precise (about 1%, 1s). C2,1, C2,–1, and C2,–2 areconsistent with a null value at the 2s level. C2,2 is about oneorder of magnitude larger than J3, which is estimated to a
relative accuracy of about 20%. The positive sign of J3 implies a negative mass anomaly at the south pole (consistent
with the observed topography), but its small value indicatesthe presence of a compensating positive mass anomaly at
depth (see section 2.3).
2.2.2. Shape determination. The shape of Enceladus has
been determined based on analysis of limb profiles collected
over the course of the Cassini mission (Thomas et al., 2007;
TABLE 1. Gravity model (Iess et al., 2014) (unnormalized,
dimensionless potential coefficients, for Rref = 252.1 km).
Parameter

Value ±1s

GMJ2 = – C2,0
C2,1
C2,–1
C2,2
C2,–2
J3 = – C3,0
J2/C2,27.210443 ± 0.00003 km3 s–25526.1 ± 35.5 × 10–6
9.4 ± 11.8 × 10–6
40.5 ± 22.8 × 10–6
1575.7 ± 15.9 × 10–6
23.0 ± 7.5 × 10–6
–118.2 ± 23.5 × 10–6
3.51 ± 0.05

64   Enceladus and the Icy Moons of SaturnThomas, 2010; Nimmo et al., 2011; Tajeddine et al., 2017).
The best fitting triaxial ellipsoidal figure exhibits substantial polar flattening, with a–c ≈ 8.5 km and a significant
departure from hydrostatic equilibrium, with (a–c)/(b–c) ≈2.8 (recall that the hydrostatic ratio is 4). Similarly, in termsof spherical harmonic coefficients, the ratio of the degree-2
zonal to sectorial terms is ~4.2, much higher than the ~3.3
expected if Enceladus were in hydrostatic equilibrium(Table 2). The limb profiles provide sufficient spatial coverage to constrain global spherical harmonic models up todegree and order 8 (Nimmo et al., 2011); the recent additionof stereogrammetric measurements has extended the latest
topography models to degree and order 16 (Tajeddine et al.,2017). One of the most significant features of the shape isthe topographic depression associated with the SPT, having
an elevation that is ~1.1 km lower than at the north pole,an effect that is expressed mainly in the H30 term (Table 2).
TABLE 2. Shape model (Nimmo et al., 2011) (unnormalized).
Parameter

Value ±1s

R
H2,0
H2,1
H2,–1
H2,2
H2,–2
H3,0
–H2,0/H2,2

252.1 km
–3846 ± 179 m
0 ± 52 m– 65 ± 52 m
917 ± 19 m
–39 ± 19 m
384 ± 5 m
4.20 ± 0.22

2.3. Modeling and Interpretation
2.3.1. Basic interpretation. Since the gravity field ofEnceladus was first measured (Iess et al., 2014), its interpretation has been an active area of research (Iess et al., 2014;
McKinnon, 2015; Čadek et al., 2016, 2017; van Hoolst et al.,2016; Beuthe et al., 2016; Hemingway and Mittal, 2017). In
part because of Enceladus’ small size and rapid rotation rate,interpretation of the data is unusually sensitive to modelingdetails, as we discuss below. Nevertheless, there are a few
basic observations that are uncontroversial.
The first is that while the shape is substantially nonhydrostatic, with H20/H22 = – 4.20 ± 0.22 (1s), the gravity
field is only modestly so, with J2/C22 = 3.51 ± 0.05 (1s).
This situation is immediately suggestive of compensation.
Without compensation, the large non-hydrostatic topography
should give rise to a correspondingly large non-hydrostatic
gravity signal. Likewise, the observed J3 gravity anomaly,(–118 ± 23) × 10–6, is substantially smaller than would be
expected if the ~1-km north-south polar elevation difference
were completely uncompensated (≈–375 × 10–6).
A compensation mechanism is therefore required. Although there are several possibilities, including systematiclateral variations in density (i.e., Pratt isostasy) or nonhydrostatic topography on the surface of the rocky core (seesection 2.3.5), the most straightforward mechanism, which

naturally yields a small non-hydrostatic gravity signal inspite of the large non-hydrostatic topography, is Airy-type
isostatic compensation. That is, Enceladus’ long wavelength
topography appears to be supported, at least in part, by thedisplacement of a relatively low viscosity, higher density
material beneath the crust — e.g., a subsurface liquid ocean.
2.3.2. Admittance analysis. The presence of nonhydrostatic degree-2 gravity and topography prevents direct
calculation of the moment of inertia via equations (16) and(21). Instead, it is necessary to first separate the observed
degree-2 gravity and topography signals into their hydrostaticand non-hydrostatic parts (Iess et al., 2014). Assuming thesecomponents can be separated linearly (a reasonable assumption as long as the non-hydrostatic parts are small comparedto the hydrostatic parts), we can write
hyd
nh
J obs
2 = J2 + J2
hyd
nh
Cobs
22 = C 22 + C 22
hyd
nh
H obs
20 = H 20 + H 20

(25)

hyd
nh
H obs
22 = H 22 + H 22

where the superscripts refer to the observations (obs), thehydrostatic expectation (hyd) given by equations (16) and(18), and the non-hydrostatic components (nh). For a given
assumed moment of inertia, the only unknowns in equation (25) are the non-hydrostatic components. The relationship between these non-hydrostatic gravity and topography
signals can be quantified by the ratios
Z20 = – J nh
H nh
2
20
nh
Z22 = Cnh
22 H 22(26)

This is effectively a component-wise version of the spectral admittance (e.g., Wieczorek, 2015), a quantity related tothe degree and depth of isostatic compensation: Greater and/or shallower compensation results in a more muted gravity
signal and thus a smaller admittance. The advantage of thisapproach is that, provided the ice shell’s mechanical properties are not significantly variable laterally, we should expectthat Z20 = Z22. Hence, we obtain self-consistent results when
we find a moment of inertia that, via equations (21), (16),(18), (25), and (26), yields Z20 = Z22.
It is worth emphasizing, however, that, for a body like
Enceladus, the expectation that Z20 = Z22 is an assumptionthat has not been explicitly tested in the literature to date. Future work to better justify this assumption, or to investigate
the effect of relaxing it to some degree, may be valuable.
Considering a wide range of possible mean moments
of inertia, Iess et al. (2014) produced the equivalent ofFig. 2a. This result follows from the first-order equilibrium
figure theory described by equation (11) — an approximation that becomes increasingly poor for fast-rotating bodies
like Enceladus (Tricarico, 2014). To account for this effect,
McKinnon (2015) repeated the analysis of Iess et al. (2014)
using a fourth-order theory of figures approach (Tricarico,2014), yielding instead the equivalent of Fig. 2b. Whereas

Hemingway et al.: The Interior of Enceladus

the first-order results indicate a preferred admittance of Z2 ≈11.8 ± 1.6 mGal km–1 and a mean moment of inertia factorof ~0.333 ± 0.002, the fourth-order results suggest insteadZ2 ≈ 15.7 ± 2.0 mGal km–1 and ~0.331 ± 0.002. While thedifference may appear subtle, the implications for compensation depth (and hence ice shell thickness) are significant.

(a)Admittance (mGal km–1)

50
40

Z22

30
20
10

Z20

0
–10
0.31

0.32

0.33

0.34

0.35

0.36

Moment of Inertia Factor

Admittance (mGal km–1)

50

(b)

40

Z22

30
20
10

Z20

0
–10
0.31

0.32

0.33

0.34

0.35

0.36

Moment of Inertia Factor
Fig. 2. Comparison of degree-2 zonal (Z20, solid line) and
sectorial (Z22, dashed line) admittances across a range ofpossible mean moments of inertia. Shaded bands indicate
1s uncertainties propagated from both the shape and gravity
models. The point of intersection represents the admittanceand moment of inertia combination that yields self-consistent
results. (a) Replication of the Iess et al. (2014) results,
following from ﬁrst-order equilibrium ﬁgure theory (except
using gravity coefﬁcients from Table 1). (b) Replication of
McKinnon (2015) results, employing fourth-order methods ofTricarico (2014). Unlike in Iess et al. (2014) and McKinnon(2015), admittances are here given in the more commonlyused units of mGal km–1 (where 1 mGal = 10–5 m s–2), obtained by converting the dimensionless gravitational potential
coefﬁcients in equation (26) to acceleration by multiplying
them by (l+1)GM/R2.

65

In addition to the degree-2 signals, the zonal terms ofthe degree-3 gravity and topography are available and canbe used to obtain the zonal part of the degree-3 admittance,
Z30 ≈ 14.0 ± 2.8 mGal km–1. The degree-3 observations
provide an independent admittance estimate that also has theadvantage of not being complicated by tidal and rotational
deformation (the theoretical hydrostatic equilibrium figure
has no degree-3 components), and is therefore independentof moment of inertia.
To put these admittance values in context, assuming adensity of rc = 925 kg m–3 for the icy crust, the degree-2 and-3 admittances expected in the case of zero compensation are
Z2 ≈ 47 mGal km–1 and Z3 ≈ 44 mGal km–1. The fact that theobserved admittances are much smaller is another indicationthat the topography is significantly compensated. Although
the degree-2 admittance values depend on the assumed mean
moment of inertia, even without requiring that Z20 = Z22,
the small values of both Z20 and Z22 suggest compensation
independently. Although Z22 becomes compatible with
uncompensated topography for moments of inertia below~0.32, this seems unlikely, as it would require Enceladus
to be one of the most strongly differentiated solid bodies inthe solar system. This observation also provides a clue about
the spatial extent of the subsurface liquid layer (assuming
Airy type compensation): Whereas a south polar regionalsea could account for the smallness of the zonal admittanceterms (Z20 and Z30), the small value of Z22 requires a more
extensive, perhaps global, subsurface liquid layer.
2.3.3. Compensation model. Further interpretation of
the admittance values requires a compensation model. Whilemost authors to date have favored Airy isostasy to explain
the compensation of Enceladus’ long wavelength topography (Iess et al., 2014; McKinnon, 2015; Čadek et al., 2016;
Beuthe et al., 2016; Hemingway and Mittal, 2017), there are
multiple ways to conceive of Airy compensation, leading
to subtle but consequential differences between the models
and their implications. Moreover, the introduction of elastic
stresses to the models widens the parameter space considerably, and can have a significant effect on the results — anissue we discuss further below.
The non-hydrostatic topography gives rise to non-hydrostatic gravitational acceleration according to (Jeffreys, 1976;
Burša and Peč, 1993)
g lm ( R ) =

l +1
R 
4πG ∆ρi Hilm  i 
 R
2l + 1
i

∑

l+2

(27)

where R is the body’s mean radius, G is the gravitational
constant, Hilm describes the degree-l and order-m non-hydrostatic topography at the top of the ith layer, whose meanouter radius is Ri, and where Dri is the density contrast
between the ith layer and the layer above it.
Assuming Airy compensation, and ignoring the role ofelastic support for the moment, the relief at the base of theice shell mirrors the non-hydrostatic surface relief with theamplitude scaled by rc/Dr, where rc is the density of thecrust (or ice shell) and Dr is the density contrast between

66   Enceladus and the Icy Moons of Saturnthe ice shell and the underlying ocean. In this case, thenon-hydrostatic topography at the top (t) and bottom (b) of
the crust are related by
H blm = –H tlm

ρc
γ
∆ρ

(28)

Here, the dimensionless factor g is a placeholder allowing
for various conceptions of Airy isostasy and different typesof elastic support, as discussed below.
In general, the model admittance (gravity/topography
ratio) can thus be writtenl+2 Rb  
l +1
4π G ρc  1 – γ Z ( l) =
2l + 1
 R t  (29)

where Rt and Rb are the radii corresponding to the top andbottom of the ice shell, respectively (the mean shell thickness being d = Rt–Rb).
Equations (27) and (29) implicitly treat the non-hydrostatic topography as a surface density anomaly — an approximation that is good when the topography is sufficiently small. In
the present case, however, the small density contrast between
the ice shell and the ocean leads to substantial relief at theice/ocean interface. Taking this finite-amplitude into account(e.g., Martinec, 1994; Wieczorek and Phillips, 1998), one
obtains a slightly smaller admittance. Although we neglect
this effect in our discussion here, within the context of any
particular compensation model, the effect can be significant,and so some authors do take it into account (e.g., Čadek etal., 2016; Hemingway and Mittal, 2017).
In the limit of zero elastic support, and when isostasy
is defined as requiring equal masses in columns of equal
solid angle, then g = (Rt/Rb)2 (e.g., Wieczorek, 2015); therequirement of equal pressures at equal depths instead leads
to g = (gt/gb), where gt and gb are the mean gravitational
acceleration at the top and bottom of the ice shell, respectively (Hemingway and Matsuyama, 2017); when the total
weight of the surface topography is required to equal thebuoyancy force of the basal relief, then g = (gt/gb)(Rt/Rb)2(e.g., Čadek et al., 2016). Unfortunately, these different
definitions of isostasy lead to significantly different estimatesof the compensation depth.
Assuming a crustal density of rc = 925 kg m–3, the observed degree-3 admittance suggests a compensation depthof between 18 km and 30 km, depending on how one defines isostasy [the “equal masses” model leads to the largest
compensation depth estimate; the “equal pressures” modelleads to the smallest (Hemingway and Matsuyama, 2017)],in agreement with the result following from the degree-2
admittance obtained using first-order equilibrium figure
theory (Iess et al., 2014). The degree-2 admittance obtainedusing fourth-order equilibrium figure theory (Tricarico,
2014; McKinnon, 2015) instead suggests a compensation
depth of between 24 km and 47 km, again, depending on
how isostasy is defined. Accounting for uncertainties in theobserved shape and gravity, the range of possible compensa-

tion depths expands even further: from ~14 km at the low
end of the estimated degree-3 admittance to ~53 km at thehigh end of the estimated degree-2 admittance.
When elastic support is included, g decreases in the caseof surface loading (e.g., due to impacts or sedimentation/
erosion at the surface) and increases in the case of basalloading (e.g., freezing/melting/relaxation at the base of thecrust), and can have a significant effect on the resulting compensation depth estimate (Hemingway and Mittal, 2017). For
instance, adding even an effective 200-m elastic layer candecrease the estimated compensation depth by more than afactor of 2 (Čadek et al., 2016), assuming the shell thickness
variations are generated at the surface; the compensation
depth estimate increases if the shell thickness variations
are instead generated at the base of the shell (Beuthe et al.,2016; Hemingway and Mittal, 2017).
Further details of elastic/viscoelastic support models arebeyond the scope of this chapter, but it should be notedthat this remains an active area of research (Čadek et al.,2016; Soucek et al., 2016; Hemingway and Mittal, 2017).
In particular, it is not yet clear how elastic support of thelong-wavelength topography is affected by factors such aslateral variations in the shell’s elastic properties (Beuthe,
2008; Čadek et al., 2016, 2017), or the way bending andmembrane stresses are transmitted across the SPT given
the presence of the major fracture systems (i.e., the Tiger
Stripes) in that region (e.g., Soucek et al., 2016).
Although the range of possible mean ice shell thicknesses(i.e., mean compensation depths) can be narrowed by arguing for one or another compensation model, researchers have
yet to converge on the best approach, making it difficult tomake definitive statements about the ice shell thickness.
However, assuming some version of Airy compensation,a few things can be stated with confidence. First, the thinnest part of the ice shell must be located beneath the large
topographic depression at the south pole. Since the shell
thickness is necessarily greater than zero there, this provides an effective lower bound on the mean shell thickness.
Depending on the definition of isostasy, and the assumed
ice shell and ocean densities, one finds that the mean shell
thickness must be at least ~18 km in order to ensure nonzero shell thickness at the south pole. An upper bound on theshell thickness is harder to establish. However, as long as thedegree-3 observations are taken into account, the mean shell
thickness is unlikely to exceed ~44 km. Shell thicknesses
greater than 44 km would lead to degree-3 admittances morethan 2s larger than the observed value, even when using
the isostasy model that produces the smallest admittances[i.e., the “equal masses” model, with g = (Rt/Rb)2]. Again,
however, the possibility of substantial elastic support ofbasal topography can allow for larger mean shell thicknesses.
2.3.4. Internal structure. The simplest approach to modeling the internal density structure is to start by assuming a
two-layer body. If the density of the outer H2O layer (oceanplus ice shell) is prescribed, then with the known bulk density (1609 kg m–3) and the moment of inertia obtained fromthe admittance analysis, we can use equations (1) and (2) to67

Hemingway et al.: The Interior of Enceladus

solve for the core radius and density (Fig. 3). A mean coreradius of ~192 km leaves ~60 km for the mean thickness ofthe H2O layer. If the mean compensation depth (i.e., ice shell
thickness) is ~30 km, this leaves a mean ocean thickness of~30 km, neglecting the effect of the small difference in densities between the ice shell and the ocean (Iess et al., 2014;
McKinnon, 2015). More sophisticated three-layer models
have also been constructed (Čadek et al., 2016; Hemingway
and Mittal, 2017), yielding broadly similar results (Table 3).
A few basic conclusions emerge immediately. First, thelow density of the (presumably silicate) core suggests hydrothermal alteration (e.g., serpentinization) and/or substantial
porosity, which may be easy to maintain given the modest
internal pressures (Fig. 4). In particular, an unconsolidated
rubble pile core, with water-filled pores, is possible (Roberts,2015). The combined thickness of the H2O layers is reasonably well constrained, but the models permit some tradeoff
between the thicknesses of the ice shell and the ocean (e.g.,
Hemingway and Mittal, 2017). Although these tradeoffs
have little effect on the density and therefore the internalpressure and gravity profiles (Fig. 4), they do make it difficult to determine the shell thickness with precision. This isregrettable since the thickness of the ice shell has a numberof important implications. We discuss the ice shell’s possible
structure, dynamics, and related implications in section 4.
2.3.5. Additional considerations. In the models discussed so far, the core shape has been assumed to conformto an equipotential, as expected if the core shape was established when the interior was warm. The small internal
gravity and pressure, however, means that the surface of thecore could, in principle, support considerable non-hydrostatic
topography over long timescales. Indeed, an irregular coreshape had previously been proposed to account for the3600

non-hydrostatic shape of Enceladus’ surface (Thomas et al.,2007; McKinnon, 2013) and to explain the rotation state of
the similarly sized Mimas (Tajeddine et al., 2014). However, while the presence of non-hydrostatic core topography
would contribute to the observed gravity, there is no reason
to expect it to do so in precisely such a way as to offset thegravity anomalies associated with the surface topography.
That is, there is no reason to expect the non-hydrostatic core
topography to effectively mirror the non-hydrostatic surface
topography, as would be required for it to be substantially
responsible for the compensating effect observed in thegravity data. In fact, it is just as likely, if not more so, that
non-hydrostatic core topography would have the opposite
effect, as would be expected in the case of an overly oblateor tidally elongated core. While the core shape may not be
precisely hydrostatic, the hydrostatic assumption may bejustified on the grounds that it is the most parsimonious.
Although it has not been done in the literature to date, some
unknown non-hydrostatic core topography could be added
as a free parameter, and would increase uncertainties in theother estimated parameters.
Whereas all the models discussed so far assume someversion of Airy-type isostasy, in which the surface topography is supported in part by lateral variations in the ice shell’s
thickness, it should be noted that the topography could also
be supported in part by lateral variations in density [i.e.,
Pratt-type isostasy (e.g., Besserer et al., 2013; Tajeddine
et al., 2017)]. This mechanism is unlikely to be primarily
responsible for the observed compensation, however, because it would require parts of the heavily cratered plains
to exhibit a higher crustal density, which is the opposite ofwhat one would expect (e.g., McKinnon, 2015; Hemingway
and Mittal, 2017) and because it would require the density
variations to extend through an unrealistically large depthof ~40 km (Hemingway and Mittal, 2017). Nevertheless, thepotential compensating effect of lateral variations in density
cannot be ruled out entirely.

3000
2800
2600
2400
2200
2000
1800
0

50

100

150

200

250

Core Radius (km)Fig. 3. Core radius and density over a range of possible
mean moments of inertia (I/MR 2 = 0.30–0.40), assuming a two-layer body with outer layer density ﬁxed at r =950 kg m–3. The estimated mean moment of inertia (0.331 ±0.002) corresponds to a core radius of ~192 km and a coredensity of ~2450 kg m–3.

40

0.16

35

0.14

30

0.12

25

0.1

20

0.08

15

0.06

10

0.04

5

0.02

0

0

50

100

150

Radius (km)

200

0
250

Gravitational Acceleration (m s–2)

3200

Pressure (MPa)Core Density (kg m–3)

3400

Fig. 4. Interior pressure (solid line) and gravity (dashedline) proﬁles for a three-layer model of Enceladus with 22km mean thickness ice shell (900 kg m–3), 38-km mean
thickness ocean (1030 kg m–3), and 192-km radius core(~2450 kg m–3).

68   Enceladus and the Icy Moons of Saturn
TABLE 3. Comparison of published interior models.
Model

Core

Ocean

Ice Shell

Notes

Iess et al.
(2014)~190 km(~2400 kg m–3)

10–30 km(1000 kg m–3*)

30–40 km(920 kg m–3*)

Based on gravity and topography:
Hydrostatic terms from first-order equilibrium
figure theory; elastic support deemed
insignificant; both regional and global oceans
discussed, but left as an open question.

McKinnon(2015)

190–195 km(~2450 kg m–3)

~10 km(1007 kg m–3*)~50 km(925 kg m–3*)

Based on gravity and topography:
Hydrostatic terms from fourth-orderequilibrium figure theory (Tricarico, 2014);
elastic support deemed insignificant; ice shell
thickness determination is based on degree-2
data only (degree-3 data deemed to have only
regional significance).

Thomaset al. (2016)~200 km(~2300 kg m–3)

26–31 km(1000 kg m–3*)

21–26 km(~850 kg m–3)

Based on physical libration amplitude: Threelayer model with layer interfaces assumed
to be hydrostatic and defined by first-orderequilibrium figure theory; ice shell thickness
constrained by libration amplitude. Core andice shell densities adjusted to fit the moment ofinertia reported by Iess et al. [2014].

Čadek et al.
(2016)

180 –185 km(~2450 kg m–3)

~50 km(~1030 kg m–3)

18–22 km(925 kg m–3*)

Based on gravity, topography, and librations:
Elastic support (assuming top loading)introduced in order to bring gravity-based
shell thickness estimate into agreement with
Thomas et al. (2016); Airy compensation
model modified to account for radial variationin gravity; finite-amplitude correction
included (Martinec, 1994).

Van Hoolst
et al. (2016)

170 –205 km(2158–2829 kg m–3)

21–67 km(950–1100 kg m–3)

14–26 km(900–1000 kg m–3)

Based on physical libration amplitude:
Similar to Thomas et al. (2016) but
additionally accounting for uncertainties
on libration amplitude arising from rigidity
and viscoelastic behavior of the ice shell, awider range of ice and ocean densities, andthe possibility of core topography; relaxes
assumption of Airy compensation.

Beuthe et al.
(2016)

186–196 km(2350–2480 kg m–3)

34–42 km(1020 kg m–3*)

19–27 km(925 kg m–3*)

Based on gravity and topography: Threelayer model with hydrostatic terms fromsecond-order equilibrium figure theory(Zharkov, 2004); isostatic compensationachieved by minimizing deviatoric stresses(e.g., Dahlen, 1982).

Hemingway
and Mittal(2017)

188–205 km(2200–2450 kg m–3)

12–36 km(1000–1100 kg m–3)

22–41 km(850–950 kg m–3)

Based on gravity and topography: Threelayer model with equilibrium figures
determined numerically (Tricarico, 2014);
Airy compensation model modified to ensure
equal pressures at depth (Hemingway andMatsuyama, 2017); using newer shapemodel (Tajeddine et al., 2017); finiteamplitude correction included (Wieczorek andPhillips, 1998).

*Values that were prescribed rather than derived.
Thicknesses refer to mean layer thickness; lateral thickness variations are required by all models.

Hemingway et al.: The Interior of Enceladus

3. IMPLICATIONS FROM LIBRATIONSIn the previous section, we saw how tidal and rotational
forces deform synchronous satellites like Enceladus into
approximately triaxial ellipsoidal figures. We discussed how
the resulting asymmetries in the mass distribution lead todifferences among the principal moments of inertia, howthese moments of inertia are related to observable asymmetries in the gravitational field, and how the moments ofinertia are related to the satellite’s internal structure. Here
we look at how a satellite’s rotational dynamics depend on
the principal moments of inertia, focusing specifically on
how the amplitude of forced physical librations can help toconstrain the interior structure.
Most satellites in the solar system are locked in a 1:1
spin:orbit resonance, where the satellite’s mean rotation rateis equal to its orbital mean motion. In this configuration,
and assuming a circular orbit, the long axis of the satellite(i.e., its static tidal bulge) always points toward the parent
body (the central planet). The eccentricity of the satellite’s
orbit, however, causes it to move faster at pericenter andslower at apocenter, even as the rotation rate remains (nearly)
constant, leading to misalignments between the satellite’slong axis and the line connecting the two bodies. To a first
approximation, the long axis of the satellite points at theempty focus of the satellite’s elliptical orbit (Murray and
Dermott, 1999, p. 44). In a frame fixed on the planet, thesatellite thus appears to rock back and forth about its spin
axis as it orbits the planet — these apparent oscillations,called “optical librations,” explain in part why more than
50% of the nearside of the Moon is visible from Earth.
Because of the periodic misalignments between the satellite’s long axis and the satellite-planet line, the planet exerts
gravitational torques on the satellite’s static tidal bulge,
creating additional oscillations called “physical longitudinal
librations.” Whereas the amplitude of the optical libration
depends only on the orbital eccentricity, the physical libration amplitude additionally depends on the satellite’s moments of inertia, thus making its measurement a valuable
tool for constraining the internal structure [e.g., for the Moon
(Dickey et al., 1994), Mercury (Margot et al., 2007), Phobos(Willner et al., 2010; Nadezhdina and Zubarev, 2014), Janusand Epimetheus (Tiscareno et al., 2009), Mimas (Tajeddine
et al., 2014), and Enceladus (Thomas et al., 2016)].
3.1. Theory
The rotation of a satellite orbiting a central planet withmass m is described by the three Euler equations of motion(Danby, 1988; Murray and Dermott, 1999)Aω x – ( B – C ) ω y ω z = 3Gm ( C – B) YZ r 5
Bω y – ( C – A ) ω z ω x = 3Gm ( A – C ) ZX r 5
Cω z – ( A – B) ω x ω y = 3Gm ( B – A ) XY r 5

(30)

69

where A, B, and C are the satellite’s principal moments ofinertia; wx, wy, and wz are the projections of the spin vector onto the satellite’s principal axes, x, y, and z; X, Y, andZ represent the coordinates of the planet in the satellite’s
principal frame; and r is the distance between the satellite
and the planet.
When the satellite’s obliquity is negligible (i.e., when
its spin pole is normal to its orbital plane), then wx = wy =0, and the problem is reduced to two dimensions (Fig. 5).
The small obliquity expected for Enceladus (Chen andNimmo, 2011; Baland et al., 2016) makes this a reasonableapproximation for purposes of our discussion.
If we define y as the angle between the satellite’s longaxis and the line connecting the two bodies, then we have
X/r = cos y and Y/r = sin y, which reduces the third line
in equation (30) to
3 ( B – A ) Gmθ=
sin 2ψ
(31)
2 C
r3
..
where we now use q to represent the variation in the satellite’s
angular velocity. Here, q is the orientation of the satellite’slong axis with respect to a fixed inertial reference frame(a)Orbital direction
Satellite
x

ψ
y

ƒ

θ
Planet(b)

ω0 < nPericenter
Planetψ φ = 2e sin(nt)
γ

Empty
Focus

Satellite

Fig. 5. The geometry of a satellite subject to librations
illustrated (a) in a frame ﬁxed on the planet; and (b) in aframe ﬁxed on the satellite and oriented with the x-axis
pointing toward the empty focus of the satellite’s orbit. In(a), f represents the true anomaly for the satellite, while q
represents the orientation of its long axis with respect to aﬁxed inertial frame. In (b), f represents the angle between
the satellite-planet line and the line connecting the satellitewith the empty focus of its orbit [this is the optical libration
angle, not shown in (a)]; g is the physical libration. After
Tiscareno et al. (2009).

70   Enceladus and the Icy Moons of Saturn(where q =0 corresponds to pericenter). This variation in
angular velocity goes as (B–A)/C because the longitudinaltorques the planet exerts on the satellite are proportional to
the satellite’s moment of inertia difference, B–A, and because
the polar moment of inertia, C, represents the resistance tothose torques.
We also define f to be the angle between the line connecting the two bodies and the line connecting the satellite withthe empty focus of its orbit (Fig. 5b) — this is the opticallibration angle. In the absence of longitudinal torques, thesatellite’s long axis would remain pointed at the empty focus,
so that y would be equal to f. In this case, the satellite’sorientation would be given simply by the mean anomaly
q(t) = nt. The presence of longitudinal torques, however,causes small additional variations in the satellite’s orientation so that q(t) = nt + g(t), where g is the forced physicallibration. Assuming for now that n is constant (although a
discussion on perturbations
from other satellites is provided.. ..
below), we have q = g .
The satellite’s true anomaly (Fig. 5a), expanded as a Fourier series and limited to first order in eccentricity (Murray
and Dermott, 1999), is given by( )f = nt + 2e sin nt + 0 e 2(32)

It is clear from Fig. 5a that f–q = y, so that y = 2e sin nt–
g. Using the small angle approximation for y, equation (31)can be rewritten
3 ( B – A ) Gmθ=
( 2e sin nt – γ )
C
r3

(33)

.. With.. the approximation that Gm ≈ r n , and given thatq = g , this yields
3 2
γ = 2eω 02 si nt – γω 02(34)

where w0 is the natural frequency, defined as
ω0 = n3 (B – A)C(35)

Finally, substituting g = g0 sin nt into equation (34), thesolution becomesγ =2e

1 – ( n ω0 )2

sin nt(36)

The amplitude of the forced physical libration thus
depends on the natural frequency w0, which is a functionof the moments of inertia via equation (35), and hence thesatellite’s interior structure. Note that the phase of the physical libration is a function of n and w0: If n < w0, then thelibration would be in phase with the torque, while if n >w0, the libration and the torque would be 180° out of phase(Murray and Dermott, 1999, p. 216), which is the case forEnceladus and the case illustrated in Fig. 5b.

So far, we have considered only the two-body problem,
where the satellite’s motion is not perturbed by other satellites. The more general expression of the physical librationis (Rambaux et al., 2010)
γ = A d sin (ω d t + α d ) e – λt+Hi

∑ 1 – (ω ω ) sin (ω t + α )
i

i

0

2

i

i

(37)

The first term in equation (37) represents the free librationof the satellite, usually induced by episodic events such as
impacts. The free libration oscillates at ωd = ω 02 – λ 2 , also
called the free libration frequency, very close to the satellite’snatural frequency w0 [l is a damping constant (Rambaux etal., 2010)]. The amplitude Ad and the phase ad are integration
constants that are functions of the initial conditions related
to the perturbing event. Since dissipation in the satellite
will damp the free librations on a timescale of tens of years(Rambaux et al., 2010; Noyelles et al., 2011), we can assume
that the observed librations of Enceladus are forced.
The second term in equation (37) represents the forced
libration, including gravitational perturbations from othersatellites. Such perturbations introduce additional librations,where Hi, wi, and ai are the magnitude, frequency, and phaseof the ith perturbation, respectively. Equation (36) is a special
case of equation (37), corresponding to the orbital frequency,n, and an amplitude of Hi = 2e. Using the JPL Horizons
Ephemeris (http://ssd.jpl.nasa.gov/horizons.cgi), Rambaux
et al. (2010) identified three major libration frequencies
for Enceladus, with periods of 1.37 days, 3.89 years, and11.05 years, and respective amplitudes of –0.028°, 0.189°,
and 0.259° (based on a solid-body interior model). The
first signal, at the orbital period of Enceladus, is due to themean anomaly perturbation. The comparison to the semianalytical model of Enceladus’ orbit (Vienne and Duriez,1995) suggests that the remaining two signals are due to theEnceladus-Dione resonance: The first is related to Dione’s
precession of pericenter, and the second is related to theorbital libration argument of the Dione-Enceladus resonance.
Equation (37) indicates that the forced libration amplitudes, at all frequencies, depend on the natural frequency
w0, and thus on the satellite’s internal structure. However,
because the amplitude depends on 1–(wi/w0)2, the effect isonly significant for the libration at the orbital frequency, forwhich (wi/w0)2 is non-negligible (Table 4). For this reason,
in attempting to constrain the interior structure, we focus
only on the libration at the satellite’s orbital period.
3.2. Measuring Physical Librations
With the exceptions of Mercury (Margot et al., 2007) andthe Moon (Dickey et al., 1994), all measurements of physicallibration amplitudes have been carried out by photogrammetry using spacecraft imagery (Tiscareno et al., 2009; Willner
et al., 2010; Nadezhdina and Zubarev, 2014; Tajeddine et al.,2014; Thomas et al., 2016). This technique uses images toHemingway et al.: The Interior of Enceladus   71TABLE 4. Dominant signals involved in the forced librationsof Enceladus (Rambaux et al., 2010).
Period (days)

Amplitude (deg)(wi/w0)2

1.37
1418.934035.64–0.028
0.188
0.259

16.3
1.5 × 10–5
1.9 × 10–6

Because (wi/w0)2 is negligibly small for the longer-period librations, only
the libration at the orbital period is significantly determined by the internal
structure, via equation (35).

build three-dimensional coordinates of recognizable surface
features, called control points, by observing them from various angles. With sufficient coverage, a global control point
network can be built and used for mapping and for generating
digital terrain models (DTMs). Like any observation-based
technique, control points are subject to reconstruction errors.
While some are random, others are systematic, such as errors
in spacecraft position, camera pointing, and the body’s rotation model (including the effect of any physical libration).
During its time in the Saturn system, the Cassini spacecraft collected tens of thousands of images of Enceladus.
Observations were made from various phase angles anddistances, in order to study the plumes, surface geology, andorbital and rotational dynamics. The amplitude of the physicallibration of Enceladus is on the order of hundreds of metersof surface displacement at the equator. The spacecraft must
be close enough to obtain images of sufficiently high resolution to permit libration detection. However, the spacecraftmust not be too close, because when only a small part of thesurface is in view, all the predicted positions of control points
projected onto that image have nearly the same systematic
errors, preventing the fitting software from distinguishing
between libration-related offsets and camera-pointing errors. Provided that a sufficient latitudinal range is in view,
however, the ambiguity is removed because, unlike camera
pointing errors, the libration-related offsets vary with latitude.
To measure the physical librations of Enceladus, Thomas
et al. (2016) used a network of 488 control points (mostly
craters) across the surface, tracking them through a seriesof 340 Cassini Imaging Science Subsystem (ISS) images,
totaling 5873 measurements. Starting from the approximately
known positions of each control point in the satellite’s reference frame, the control points were first converted to J2000,
a fixed frame based on the initial conditions of Earth’s rotational state on January 1, 2000, at 12:00 Terrestrial Time. Thisrequires knowledge of the satellite’s rotational parameters(i.e., orientation of the spin axis, rotation rate, librations,etc.), accounting for the three major libration frequencies(Table 4). As discussed above, the amplitudes of the longperiod librations are well known since their dependence on
Enceladus’ interior is negligible. The amplitude of the signal
at the orbital period of Enceladus (including the physicallibration), however, is a parameter that needs to be fitted.
Next, the coordinates are converted to the reference frameof the Cassini ISS Narrow Angle Camera (NAC), requiring
instantaneous information of the spacecraft’s position andthe camera’s pointing angle. This information is available
at NASA’s Navigation and Ancillary Information Facility(NAIF) SPICE website (http://naif.jpl.nasa.gov/naif). The
twist angle (orientation about the axis along the line of sight)of the Cassini NAC is well determined, and the position of
the spacecraft is known to within 1–10 km. The camera’s
right ascension and declination, however, are not well determined and must be fitted.
Next, the three-dimensional coordinates of the control
points were adjusted in order to minimize the c2 residuals
between the observed and predicted positions of the control
points projected onto the camera image plane. To avoid
degeneracies associated with going from three dimensions
to two, each control point must have been observed fromat least two sufficiently different viewing geometries. The
more a control point is observed in different images, thelower the uncertainties in its position. Similarly, increasing
the number of control points in an image reduces the uncertainties on camera-pointing angle. Finally, the reconstructed
control point network still depends on the assumed rotationalmodel for the satellite. When the model does not describe
the satellite’s rotation accurately, additional errors will resultin the reconstructed positions of control points. Therefore,
the libration amplitude at the orbital period must be varied
in order to minimize the reconstruction errors.
Following this procedure, Thomas et al. (2016) obtained
a best-fit physical libration amplitude of –0.120° (the minus
sign indicates that the libration is out of phase by 180°) witha 2s uncertainty of 0.014°. This translates to 528 ± 60 m insurface displacement at the equator. The smallness of theuncertainty is a function of the fact that the control points
are based on ellipses fitted to craters, permitting subpixel
accuracy in their positions and, more importantly, the large
number of data points (since the uncertainty goes as N–1/2).
3.3. Interior Modeling and Interpretation
Having measured the libration amplitude, we can now
compute the satellite’s natural frequency via equation (34)and hence determine its dynamical triaxiality, (B–A)/C, via
equation (35). The value corresponding to the observed libration amplitude of 0.120° is (B–A)/C ≈ 0.0607. Althoughfixing this quantity does not yield a unique interior model,
various simple trial models can be tested to determine theapproximate internal structure.
Assuming a multi-layered body consisting of nested
triaxial ellipsoids, the principal moments of inertia, definedin equation (3), become
A=
B=
C=

4π
15

∑∆ρ ( b + c ) a b c4π
15∑∆ρ ( a + c ) a b c4π
15∑∆ρ ( a + b ) a b c

i

2
i

2
i

i i i

i

i

2
i

2
i

i i i

i

i

i

2
i

2
i

i i i

(38)

72   Enceladus and the Icy Moons of Saturnwhere ai, bi, and ci are the semiaxes of the ith ellipsoid, andDri is the density contrast between layer i and the layer
above it. In the simplest case of a homogeneous interior,
the dynamical triaxiality is simply
B–A
a 2 – b2
= 2
C
a + b2(39)

Based on their shape model, in which a = 256.2 km,
and b = 251.4 km, Thomas et al. (2016) found that the libration amplitude expected for a homogenous interior was~0.032° — more than 10s smaller than the observed value.
In the case of a two-layer body, the triaxiality is given by

(
(

)
)

(
(

Vρs a s2 – bs2 Vc (ρc – ρs ) a c2 – bc2
B –A=CVρs a s2 + bs2 + Vc (ρc – ρs ) a c2 + bc2

)
)

(40)

where the subscripts c and s indicate the rocky core and theicy shell, respectively, V is the volume, and r is the density.
Unlike the case of a homogeneous interior, models with two
or more layers require additional assumptions about theshape and densities of each layer. The known bulk density(1609 kg m–3 for Enceladus) and the inferred moment ofinertia factor constrain, via equations (1) and (2), the layer
thicknesses such that the core radius can be computed once
the layer densities are specified. Computing the shapes ofinternal layers is less straightforward but can be done in a
number of ways. Thomas et al. (2016) and van Hoolst et al.
(2016), for instance, integrate the Clairaut equation (Clairaut, 1743; Danby, 1988) to first order in terms of the polarand equatorial flattening to obtain expressions that depend
only on the densities and mean radii of the layers, which
can then be related to the semiaxes of the triaxial ellipsoidal
shape (for details, see Tajeddine et al., 2014). Similarly,
Tricarico (2014) describes methods for computing higherorder nested ellipsoidal figures recursively or numerically interms of polar and equatorial eccentricities, again, startingfrom the specified densities and mean radii of each layer.
Making the assumption that the core shape conforms to
the hydrostatic expectation (i.e., its surface is equipotential), and considering a range of ice shell densities between
700 kg m–3 and 930 kg m–3, Thomas et al. (2016) determinedthat the libration amplitude for the two-layer model would
be between 0.032° and 0.034°, still much smaller than theobserved value. Thomas et al. further showed that a regional
subsurface sea centered on the south pole has little effect on
the libration amplitude because of its symmetry about theaxis of rotation. Similar results were obtained by van Hoolst
et al. (2016), confirming the incompatibility of the libration
observations with interior models involving a core that is
physically coupled to the icy mantle.
As another check, one can use the above dynamical triaxiality in combination with the observed sectorial gravity
harmonic, C22 (Table 1), and equation (22) to compute thepolar moment of inertia. The result is C ≈ 0.104M R2. This
unrealistically small polar moment of inertia is not compatible with the internal density structure inferred from thegravity observations [roughly 0.22M R2 (Iess et al., 2014;
McKinnon, 2015)] and is another indication that Enceladus
is not behaving like an entirely solid body.
The large observed libration amplitude evidently requiresa more radically different interior model. Thomas et al.
(2016) thus argued for a model with a global subsurface
ocean that completely decouples the rocky core from the icy
shell. In this case, the icy shell and the rocky core experienceand respond to gravitational torques nearly independently.
In the limit of a spherical core, for example, only the shellitself experiences external gravitational torques and, because
the polar moment of inertia for the shell alone, Cs, is smaller
than that of the body as a whole, C, the libration amplitude
should be larger.
Thomas et al. (2016) and van Hoolst et al. (2016) applieda more generalized analysis, which includes the effects of
torques between the shell and the core, as well as the effectsof the ocean pressure on both the core and the shell. For abody with a subsurface ocean decoupling the shell from thecore, the libration amplitude of the shell is given by (e.g.,van Hoolst et al., 2009; Rambaux et al., 2011; Richard etal., 2014)()

2e  K s K c 2K int – n 2 Cc + 2K int K c γs =Cc Cs n 2 – ω12 n 2 – ω 22

(

)(

)

(41)

where Ks, Kc, and Kint are the planet-shell, planet-core, andcore-shell torques, respectively. The first two torques arefunctions of the moments of inertia of the shell and the core,
as well as the ocean pressure on each layer; Kint depends
on the densities and the dimensions of each layer; and w1
and w2 are the system’s natural frequencies (for details, seeRichard et al., 2014).
This expression shows that the libration amplitude once
again depends on the densities assumed for the internal layers as well as their shapes. Assuming hydrostatic figures foreach layer and assuming the ice shell and ocean densities tobe 850 kg m–3 and 1000 kg m–3, respectively, Thomas et al.
(2016) found the best agreement with the observed librations
for interior models with mean ice shell thicknesses in therange 21–26 km and mean ocean thicknesses in the range
26–31 km. The densities of the ice shell and the core were
adjusted (within the constraints noted above) to 850 kg m–3
and 2300 kg m–3, respectively, in order to be consistent withthe moment of inertia inferred from the gravity observations(Iess et al., 2014). Using similar methods, but considering
different ranges of core, ice shell, and ocean densities (2158–
2829 kg m–3, 900–1000 kg m–3, and 950–1100 kg m–3,
respectively), and allowing for non-hydrostatic figures forthe core and the ocean/ice interface, van Hoolst et al. (2016)
obtained mean ice shell and ocean thicknesses in the ranges
14–21 km and 24–67 km, respectively.

Hemingway et al.: The Interior of Enceladus   73Although results vary slightly with the assumptions that
go into such models, it is clear that the observed physical libration amplitude places powerful constraints on theinternal structure. Most importantly, these libration studies
demonstrate that the ice shell must be fully decoupled fromthe deeper interior, requiring the subsurface ocean to be
global. Additionally, the libration studies provide a meansof estimating the ice shell’s thickness in a manner that is
independent from the gravity-based analysis discussed in
section 2, increasing confidence in both approaches.
4. ICE SHELL STRUCTURE AND DYNAMICS
4.1. Structure
The structure of the ice shell is of particular interest because it bears on the thermal state of Enceladus, the means
for supporting its topography, the mechanics of the major fissures in the south polar region, and the nature of the ongoing
eruptions. Determining the precise thickness of the ice shell,
however, is not necessarily straightforward. For example,
analysis of Enceladus’ gravity field (section 2) has yieldedmean shell thickness estimates ranging from 20–40 km (Iesset al., 2014) at the low end, to 40–60 km (McKinnon, 2015)
at the high end (although the latter estimate is based only
on the degree-2 gravity). The large amplitude of the diurnal
forced physical librations (section 3) (Thomas et al., 2016;van Hoolst et al., 2016), on the other hand, suggests a meanshell thickness of 15–25 km. The apparent discrepancy is
not as problematic as it may seem, however, as it can be
resolved in a number of ways involving details of the elastic/
isostatic compensation model (Čadek et al., 2016; Beuthe
et al., 2016; Hemingway and Mittal, 2017), as discussed in
section 2.3.3.
What is common among the various models is thatthey take the long-wavelength topography to be related tolateral variations in the thickness of the icy shell, which is,
to some degree, supported isostatically (i.e., it is “floating”
on a subsurface ocean). Hence, the large topographic basin
at the south pole implies a regional thinning of the ice
shell, as predicted by Collins and Goodman (2007) basedon the observed shape (Thomas et al., 2007) and the localized geologic activity. The gravity observations (Iess et al.,2014) confirmed this prediction by showing that the small
magnitude of the corresponding south polar gravity anomaly
implied substantial compensation (i.e., thinning at the base of
the ice shell). Given the ~2-km topographic depression at thesouth pole (measured with respect to an equilibrium figure)(Nimmo et al., 2011), and making reasonable assumptions
about the ice shell and ocean densities, the total crustal thinning (relative to the mean shell thickness) at the south polemust be roughly 16–18 km, assuming complete Airy-type
isostatic compensation (see equation (28) in section 2.3.3)(Čadek et al., 2016; Hemingway and Mittal, 2017). Since
the shell thickness is evidently greater than zero at thesouth pole, this may be regarded as an approximate lower

bound on the mean shell thickness. However, the precise
amplitude of lateral shell thickness variations also depends
on the assumed compensation model (i.e., pressure vs. force
balance, inclusion of partial elastic support, lateral variations
in crustal density, etc.) (Čadek et al., 2016; Hemingway andMatsuyama, 2017; Hemingway and Mittal, 2017; Tajeddine
et al., 2017).
The various models also agree that the thickest part of theshell is around the equator, and especially near the sub- andanti-saturnian points, where the thickness may be 30–40 km,
or even more, depending on the preferred compensation
model (Čadek et al., 2016; Hemingway and Mittal, 2017).
The overall pattern (Fig. 6) may be the result, at least inpart, of tidal heating. Tidal dissipation is strongest at thepoles and weaker in the equatorial regions (Ojakangas andStevenson, 1989). The equilibrium shell thickness basedon tidal dissipation and insolation results in a pattern withpower at degrees 2 and 4 (Ojakangas and Stevenson, 1989;
Hemingway and Mittal, 2017). Higher-order structure in
the shell thickness variations may be a result of heterogeneities in the shell and/or mode coupling associated with
non-Newtonian flow of the ice as it relaxes (Nimmo, 2004).
4.2. Thermal and Mechanical Stability
The structure of the ice shell raises questions about thestability of the current configuration. Lateral shell thickness
variations induce stresses that will tend to relax away those
thickness variations over time, with the ice flowing fromthe thicker regions (i.e., the equator) to the thinner regions(i.e., the poles). The observed structure is thus stable only
if the relaxation rate is counterbalanced by spatial variationsin ice melting (associated with tidal dissipation) and ocean
crystallization rates (e.g., Collins and Goodman, 2007;
Kamata and Nimmo, 2017). The rate at which the deflected
ice/ocean interface relaxes toward an equipotential surface
is controlled by the viscosity structure of the ice shell. Following the approach of Lefevre et al. (2014), developed forTitan, Čadek et al. (2016) determined that the relaxation
velocity remains smaller than a few centimeters per year
only if the ice shell is conductive with a bottom viscosity ofat least 1014−1015 Pa s. Counterbalancing such a relaxation
rate requires heat flux variations between the equatorial andpolar regions of at least 100–200 mWm–2, which is comparable with the estimated heat flux in the SPT.
It is also possible that the current configuration istransient, with the observed surface topography being thepartially relaxed remnants of a previous state with larger
thickness variations (Čadek et al., 2016). Čadek et al. (2017)
have argued, using a viscoelastic relaxation model, that theobserved surface topography may be presently relaxing
toward a new equilibrium if Enceladus experiences changes
in heat production and ice melting/ocean crystallization ontimescales comparable to the relaxation timescale (which
may range between 1 and 100 m.y. for lithospheric viscosities between 1022 and 1024 Pa s). Indeed, the intense

74

Enceladus and the Icy Moons of Saturn
90

Latitude (ºN)

60
30
0
–30
–60
–90

0

45

90

135

180

225

270

315

360

Longitude (ºE)
Fig. 6. Lateral variations in ice shell thickness assuming complete Airy compensationof all known topography [up to spherical harmonic degree 8 (Nimmo et al., 2011)].
Contours indicate shell thickness in kilometers. In this example, the ice shell thickness ranges from ~6 km at the south pole to ~36 km at the sub- and anti-saturnian
points along the equator, with a mean thickness of 22 km. While the mean shell
thickness is somewhat model-dependent, the amplitude of lateral variation does notvary signiﬁcantly among the models (e.g., Čadek et al., 2016; Hemingway and Mittal,
2017). The addition of partial elastic support, however, would alter the amplitude ofshell thickness variations, especially at the shortest wavelengths.

surface activity in the south polar region is an indicationof a dynamical world where the ice shell structure may be
evolving on even shorter timescales.
4.3. Implications
The shell thickness implies a small Rayleigh number,making convection unlikely and suggesting a thermallyconductive ice shell. Heat flux, F, can be related to thetemperature structure and thickness of a conductive ice
shell according toF=c  Ts lnd  Tb (42)where d is the shell thickness, Ts and Tb are the surface andbasal temperatures, respectively, and where c is an empirically derived constant, taken to be 567 W/m (e.g., Klinger,
1980; Nimmo, 2004). For an ice shell thickness of about
30–40 km in the equatorial region, the diffusive heat flux
typically ranges between 15 and 20 mWm–2. As the diffusive
heat loss is inversely proportional to the ice shell thickness,cooling is even faster away from the equator, where the shellis thinner. For an average ice shell thickness between 20 and
25 km, the total diffusive heat loss outside the SPT is about
20 GW. In addition to the passive heat loss, one must also
consider the heat loss associated with the intense activity
at the south pole, which is estimated to be up to 15 GW(Howett et al., 2011). For ice shell thicknesses ranging be-

tween 3 and 5 km in the SPT, about 5 to 10 GW is lost justby thermal diffusion through the ice shell, with the rest of
the power being associated with the eruption activity itself.
In total, then, some ~35 GW must be generated to maintainthe present-day state of Enceladus.
Such a large power is consistent with the latest estimates ofthe dissipation function in Saturn (Lainey et al., 2012; Lainey,2016), which indicates that a strong dissipation (20 GW oreven more) may be sustained during relatively long periodsof time before the orbital eccentricity would be damped.
However, the mechanism for generating such power within
Enceladus is less clear (see the chapter in this volume by
Nimmo et al.). Such heat production in the ice shell would
require a very low bottom viscosity [<5 × 1013 Pa s (Roberts
and Nimmo, 2008)], which would result in rapid relaxation of
the basal topography, incompatible with the observed thickness variations. On the other hand, for a bottom viscosity
ranging between 1014 and 1015 Pa s, heat flow due to tidal
dissipation in the conductive ice shell would barely exceed
a few milliwatts per square meter (Roberts and Nimmo,
2008; Čadek et al., 2016), too weak to counterbalance thediffusive heat loss and to prevent the ocean from freezing(the timescale to crystallize 10 km of ice in the thickest partof the ice shell is on the order of 10 m.y.). The present-day
state may result from enhanced tidal dissipation due to higher
eccentricities in the recent past, as suggested by Běhounková
et al. (2012). This would imply that the ocean is presently
crystallizing with a maximum crystallization rate in the equatorial region where tidal heating is at a minimum.

Hemingway et al.: The Interior of Enceladus   75An alternative solution would be strong dissipation in
the deeper interior that counterbalances the diffusive heatloss through the ice shell. Two candidate processes would
be dissipation of resonant tidal waves in the ocean (Tyler,
2011; Matsuyama, 2014; Hay and Matsuyama, 2017) or
tidal friction in the unconsolidated core (Roberts, 2015). The
former process requires a thin ocean, contrary to the ocean
thickness inferred by recent studies (Čadek et al., 2016;van Hoolst et al., 2016; Beuthe et al., 2016; Hemingway
and Mittal, 2017). For the latter process, an unconsolidated
rocky core filled with water ice, the power produced is toolow (Roberts, 2015). Dissipation in a water-filled porous
core might be more efficient, but more work is needed todemonstrate whether this would be sufficient to counterbalance the diffusive heat loss. Interestingly, strong dissipationin the core may also provide the energy source to power
the hydrothermal activity, as suggested by the detection ofnanosilica emitted from Enceladus (Hsu et al., 2015), andto sustain long-term circulation of hot water in the core, as
modeled by Travis and Schubert (2015).
An interesting consequence of the reduced ice shell thickness in the SPT is a strong increase of tidal deformation andassociated tidal friction. Compared to models with constant
thickness of 25 km, the reduction in the SPT to less than
5 km in thickness results in an amplification of the tidal
stresses by a factor of 4 (Soucek et al., 2017). This is aconsequence of the small size of Enceladus, which makes
the amplitude of tidal deformation much more sensitive tothe ice shell thickness than in larger moons like Europa andTitan, where tidal deformation depends only slightly on theice shell thickness (Tobie et al., 2005). Soucek et al. (2016)
have also demonstrated that the presence of faults further
enhances the tidal deflections by at least a factor of 2. The
resulting stress patterns are much more complex than thosepredicted from standard tidal deformation models based ona thin-shell approximation and neglecting the presence offaults (e.g., Hurford et al., 2007; Nimmo et al., 2007). Enhanced tidal stresses in the SPT, together with lithosphericstresses resulting from ice shell melting and subsequentrelaxation, may help explain the tectonic patterns observedin the SPT (e.g., Patthoff and Kattenhorn, 2011; Yin et al.,2016). However, future modeling efforts are required to better understand the complex interplay between the ice shellevolution and the tectonic patterns.
5. SUMMARY AND OPEN QUESTIONSThe Cassini mission has helped to answer many questions about the interior of Enceladus. In particular, as we
discussed in this chapter, Cassini-derived measurements of
the low-order gravity field and rotational dynamics have
led to a good, albeit basic, understanding of the internal
structure (Fig. 1). We now know, for instance, that Enceladus
has a large, low-density core that is in contact with a global
subsurface liquid water ocean containing some 107 km3 ofwater (section 2.3.4) — comparable to the volume of Earth’s
Arctic Ocean. The icy shell that covers the ocean may bevanishingly thin at the south pole, but is at least a few tensof kilometers thick in the equatorial regions (Fig. 6).
Nevertheless, the thickness, structure, and dynamics of
the ice shell are not yet sufficiently well constrained to give
us a clear understanding of the thermal history and evolution of the interior. Is the current configuration stable, or
are we seeing Enceladus in the midst of a transition? If thepresent configuration is stable, how is the long-wavelength
topography maintained in spite of the relatively high temperatures at the base of the ice shell? How and where is
the internal heat generated, and can we match the predicted
input power with the estimated heat loss? In the near term,
additional modeling efforts may help to better address thesequestions. In the longer term, when spacecraft next visit
Enceladus, and especially with dedicated missions that canorbit or even land on Enceladus, new observations — which
may include higher-order gravity measurements, improved
heat flow measurements, ice penetrating radar, and evenseismology — are sure to be the most powerful drivers inadvancing our understanding of the interior of Enceladus.
Acknowledgments. D.H. was supported by the Miller Institute for Basic Research in Science at the University of California, Berkeley. R.T. was supported by the Cassini mission. L.I.
acknowledges support from the Italian Space Agency. G.T. was
supported by the European Research Council under the European
Community’s Seventh Framework Programme, FP7/20077-2013,
grant agreement 259285. We thank F. Nimmo and A. Ermakov forconstructive reviews and M. Manga for helpful comments.

REFERENCES
Anderson J. D., Schubert G., Jacobson R. A., Lau E. L., Moore W. B.,
and Sjogren W. L. (1998) Europa’s differentiated internal structure:
Inferences from four Galileo encounters. Science, 281(5385), 2019–
2022, DOI: 10.1126/science.281.5385.2019.
Baland R. M., Yseboodt M., and Van Hoolst T. (2016) The obliquity ofEnceladus. Icarus, 268, 12–31, DOI: 10.1016/j.icarus.2015.11.039.
Baum W. A., Kreidl T., Westphal J. A., Danielson G. E., Seidelmann
P. K., Pascu D., and Currie D. G. (1981) Saturn’s E ring:
I. CCD observations of March 1980. Icarus, 47(1), 84–96, DOI:
10.1016/0019-1035(81)90093-2.
Běhounková M., Tobie G., Choblet G., and Čadek O. (2012)
Tidally-induced melting events as the origin of south-pole
activity on Enceladus. Icarus, 219(2), 655–664, DOI: 10.1016/j.
icarus.2012.03.024.
Besserer J., Nimmo F., Roberts J. H., and Pappalardo R. T. (2013)
Convection-driven compaction as a possible origin of Enceladus’slong wavelength topography. J. Geophys. Res.–Planets, 118(5),908–915, DOI: 10.1002/jgre.20079.
Beuthe M. (2008) Thin elastic shells with variable thickness forlithospheric flexure of one-plate planets. Geophys. J. Intl., 172(2),817–841, DOI: 10.1111/j.1365-246X.2007.03671.x.
Beuthe M., Rivoldini A., and Trinh A. (2016) Enceladus’ and Dione’s
floating ice shells supported by minimum stress isostasy. Geophys.
Res. Lett., 43, DOI: 10.1002/2016GL070650.
Brown R. H. et al. (2006) Composition and physical properties ofEnceladus’ surface. Science, 311(5766), 1425–1428, DOI: 10.1126/
science.1121031.
Buratti B. and Veverka J. (1984) Voyager photometry of Rhea,
Dione, Tethys, Enceladus and Mimas. Icarus, 58, 254–264, DOI:
10.1016/0019-1035(84)90042-3.
Burša M. and Peč K. (1993) Gravity Field and Dynamics of the Earth.
Springer, New York. 333 pp.

76   Enceladus and the Icy Moons of SaturnČadek O. et al. (2016) Enceladus’s internal ocean and ice shell
constrained from Cassini gravity, shape and libration data. Geophys.
Res. Lett., 43, DOI: 10.1002/2016GL068634.
Čadek O., Běhounková M., Tobie G., and Choblet G. (2017)
Viscoelastic relaxation of Enceladus’s ice shell. Icarus, 291, 31–35,
DOI: 10.1016/j.icarus.2017.03.011.
Chen E. M. A. and F. Nimmo (2011) Obliquity tides do not significantly
heat Enceladus. Icarus, 214(2), 779–781, DOI: 10.1016/j.
icarus.2011.06.007.
Clairaut A. C. (1743) Theorie de la figure de la Terre, Tirée desprincipes de l’Hydrostatique, Paris.
Clark R. N., Brown R. H., Nelson M. L., and Hayashi J. N. (1983)
Surface composition of Enceladus. Bull. Am. Astron. Soc.,
15, 853.
Collins G. C. and Goodman J. C. (2007) Enceladus’ south polar sea.
Icarus, 189(1), 72–82, DOI: 10.1016/j.icarus.2007.01.010.
Dahlen F. A. (1982) Isostatic geoid anomalies on a sphere. J. Geophys.
Res., 87(B5), 3943–3947, DOI: 10.1029/JB087iB05p03943.
Danby J. M. A. (1988) Fundamentals of Celestial Mechanics. WillmanBell, Richmond. 348 pp.
Darwin G. H. (1899) The theory of the figure of the Earth carried tothe second order of small quantities. Mon. Not. R. Astron. Soc., 60,
82–124.
Dickey J. O. et al. (1994) Lunar Laser Ranging: A continuing legacy of
the Apollo program. Science, 265(5171), 482–490.
Dougherty M. K., Khurana K. K., Neubauer F. M., Russell C. T., SaurJ., Leisner J. S., and Burton M. E. (2006) Identification of a dynamic
atmosphere at Enceladus with the Cassini magnetometer. Science,
311(5766), 1406–1409, DOI: 10.1126/science.1120985.
Gladman B., Quinn D. D., Nicholson P., and Rand R. (1996)
Synchronous locking of tidally evolving satellites. Icarus, 122,
166–192, DOI: 10.1006/icar.1996.0117.
Haff P. K., Eviatar A., and Siscoe G. L. (1983) Ring and plasma: The
enigmae of Enceladus. Icarus, 56(3), 426–438, DOI: 10.1016/00191035(83)90164-1.
Hansen C. J., Esposito L., Stewart A. I. F., Colwell J., Hendrix A., Pryor
W., Shemansky D., and West R. (2006) Enceladus’ water vapor
plume. Science, 311, 1422–1425, DOI: 10.1126/science.1121254.
Hay H. C. F. C. and Matsuyama I. (2017) Numerically modelling tidal
dissipation with bottom drag in the oceans of Titan and Enceladus.
Icarus, 281, 342–356, DOI: 10.1016/j.icarus.2016.09.022.
Hemingway D. J. and Matsuyama I. (2017) Isostatic equilibrium inspherical coordinates and implications for crustal thickness on theMoon, Mars, Enceladus, and elsewhere. Geophys. Res. Lett., 44,7695–7705, DOI: 10.1002/2017GL073334.
Hemingway D. J. and Mittal T. (2017) What explains the structure of
Enceladus’s ice shell and can it be in equilibrium? Abstract P43B2877 presented at 2017 Fall Meeting, AGU, New Orleans, Louisiana,
11–15 December.
Hemingway D., Nimmo F., Zebker H., and Iess L. (2013a) A rigid andweathered ice shell on Titan. Nature, 500(7464), 550–552, DOI:
10.1038/nature12400.
Hemingway D., Nimmo F., and Iess L. (2013b) Enceladus’ internal
structure inferred from analysis of Cassini-derived gravity andtopography. Abstract P53E-03 presented at 2013 Fall Meeting, AGU,
San Francisco, California, 9–13 December.
Hemingway D. J., Zannoni M., Tortora P., Nimmo F., and Asmar S. W.
(2016) Dione’s internal structure inferred from Cassini gravity andtopography. In Lunar and Planetary Science XLVII, Abstract #1314.
Lunar and Planetary Institute, Houston.
Howett C. J. A., Spencer J. R., Pearl J., and Segura M. (2011) Highheat flow from Enceladus’ south polar region measured using 10–
600 cm–1 Cassini/CIRS data. J. Geophys. Res.–Planets, 116(3), 1–15,
DOI: 10.1029/2010JE003718.
Hsu H.-W. et al. (2015) Ongoing hydrothermal activities within
Enceladus. Nature, 519(7542), 207–210, DOI: 10.1038/nature14262.
Hubbard W. B. (1984) Planetary Interiors. Van Nostrand Reinhold, New
York. 334 pp.
Hurford T. A., Helfenstein P., Hoppa G. V., Greenberg R., and Bills
B. G. (2007) Eruptions arising from tidally controlled periodic
openings of rifts on Enceladus. Nature, 447(7142), 292–294, DOI:
10.1038/nature05821.
Iess L., Rappaport N. J., Jacobson R. A., Racioppa P., Stevenson D. J.,
Tortora P., Armstrong J. W., and Asmar S. W. (2010) Gravity field,
shape, and moment of inertia of Titan. Science, 327(5971), 1367–
1369, DOI: 10.1126/science.1182583.

Iess L., Jacobson, Ducci M., Stevenson D. J., Lunine J. I., Armstrong
J. W., Asmar S. W., Racioppa P., Rappaport N. J., and Tortora
R. A. P. (2012) The tides of Titan. Science, 337, 457–459, DOI:
10.1126/science.1219631.
Iess L. et al. (2014) The gravity field and interior structure of Enceladus.
Science, 344(6179), 78–80, DOI: 10.1126/science.1250551.
Jacobson R. A. (2004) The orbits of the major saturnian satellites
and the gravity field of Saturn from spacecraft and Earth-based
observations. Astron. J., 128, 492–501, DOI: 10.1086/421738.
Jeffreys H. (1976) The Earth: Its Origin, History, and Physical
Constitution, 6th edition. Cambridge Univ., New York. 612 pp.
Kamata S. and Nimmo F. (2017) Interior thermal state of Enceladus
inferred from the viscoelastic state of the ice shell. Icarus, 284,
387–393, DOI: 10.1029/2011JE003835.
Klinger J. (1980) Influence of a phase transition of ice on the heat andmass balance of comets. Science, 209(4453), 271–272.
Kozai Y. (1976) Masses of satellites and oblateness parameters ofSaturn. Publ. Astron. Soc. Japan, 28, 675–691.
Lainey V. (2016) Quantification of tidal parameters from solar systemdata. Cel. Mech. Dyn. Astron., 126(1), 1–12, DOI: 10.1007/s10569016-9695-y.
Lainey V. et al. (2012) Strong tidal dissipation in Saturn and constraints
on Enceladus’ thermal state from astrometry. Astrophys. J., 752(14),
DOI: 10.1088/0004-637X/752/1/14.
Lefevre A., Tobie G., Choblet G., and Cadek O. (2014) Structure anddynamics of Titan’s outer icy shell constrained from Cassini data.
Icarus, 237, 16–28, DOI: 10.1016/j.icarus.2014.04.006.
Margot J. L., Peale S. J., Jurgens R. F., Slade M. A, and Holin I. V.
(2007) Large longitude libration of Mercury reveals a molten core.
Science, 316(5825), 710–714, DOI: 10.1126/science.1140514.
Martinec Z. (1994) The density contrast at the Mohorovičic̀
discontinuity. Geophys. J. Intl., 117(2), 539–544, DOI:
10.1111/j.1365-246X.1994.tb03950.x.
Matsuyama I. (2014) Tidal dissipation in the oceans of icy satellites.
Icarus, 242, 11–18, DOI: 10.1016/j.icarus.2014.07.005.
McKinnon W. B. (2013) The shape of Enceladus as explained by anirregular core: Implications for gravity, libration, and survival of its
subsurface ocean. J. Geophys. Res.–Planets, 118, 1775–1788, DOI:
10.1002/jgre.20122.
McKinnon W. B. (2015) Effect of Enceladus’s rapid synchronous spinon interpretation of Cassini gravity. Geophys. Res. Lett., 42, DOI:
10.1002/2015GL063384.
Meyer J. and Wisdom J. (2007) Tidal heating in Enceladus. Icarus,
188(2), 535–539, DOI: 10.1016/j.icarus.2007.03.001.
Murray C. D. and Dermott S. F. (1999) Solar System Dynamics.
Cambridge Univ., Cambridge.
Nadezhdina I. E. and Zubarev A. E. (2014) Formation of a reference
coordinate network as a basis for studying the physical parametersof phobos. Solar System Res., 48(4), 269–278, DOI: 10.1134/
S003809461404008X.
Nimmo F. (2004) Non-Newtonian topographic relaxation on Europa.
Icarus, 168(1), 205–208, DOI: 10.1016/j.icarus.2003.11.022.
Nimmo F., Spencer J. R., Pappalardo R. T., and Mullen M. E.
(2007) Shear heating as the origin of the plumes and heat flux on
Enceladus. Nature, 447(7142), 289–291, DOI: 10.1038/nature05783.
Nimmo F., Bills B. G., and Thomas P. C. (2011) Geophysicalimplications of the long-wavelength topography of theSaturnian satellites. J. Geophys. Res., 116, E11001, DOI:
10.1029/2011JE003835.
Noyelles B., Karatekin O., and Rambaux N. (2011) The rotationof Mimas. Astron. Astrophys., 536, 1–13, DOI: 10.1051/00046361/201117558.
Ojakangas G. W. and Stevenson D. J. (1989) Thermal state of an ice
shell on Europa. Icarus, 81(2), 220–241, DOI: 10.1016/00191035(89)90053-5.
Pang K. D., Voge C. C., Rhoads J. W., and Ajello J. M. (1984) The
E ring of Saturn and satellite Enceladus. J. Geophys. Res., 89(B11),
9459, DOI: 10.1029/JB089iB11p09459.
Patthoff D. A. and Kattenhorn S. A. (2011) A fracture history on
Enceladus provides evidence for a global ocean. Geophys. Res. Lett.,
38(18), 1–6, DOI: 10.1029/2011GL048387.
Porco C. C. et al. (2006) Cassini observes the active south poleof Enceladus. Science, 311(5766), 1393–1401, DOI: 10.1126/
science.1123013.
Postberg F., Kempf S., Schmidt J., Brilliantov N., Beinsen A., Abel B.,
Buck U., and Srama R. (2009) Sodium salts in E-ring ice grains

Hemingway et al.: The Interior of Enceladus   77from an ocean below the surface of Enceladus. Nature, 459, 1098–
1101, DOI: 10.1038/nature08046.
Postberg F., Schmidt J., Hillier J., Kempf S., and Srama R. (2011) A
salt-water reservoir as the source of a compositionally stratified
plume on Enceladus. Nature, 474(7353), 620–622, DOI: 10.1038/
nature10175.
Rambaux N., Castillo-Rogez J. C., Williams J. G., and Karatekin Ö.
(2010) Librational response of Enceladus. Geophys. Res. Lett., 37(4),
1–6, DOI: 10.1029/2009GL041465.
Rambaux N., Van Hoolst T., and Karatekin Ö. (2011) Librational
response of Europa, Ganymede, and Callisto with an ocean
for a non-Keplerian orbit. Astron. Astrophys., 527, A118, DOI:
10.1051/0004-6361/201015304.
Rappaport N. J., Iess L., Tortora P., Anabtawi A., Asmar S. W., Somenzi
L., and Zingoni F. (2007) Mass and interior of Enceladus from
Cassini data analysis. Icarus, 190(1), 175–178, DOI: 10.1016/j.
icarus.2007.03.025.
Richard A., Rambaux N., and Charnay B. (2014) Librational responseof a deformed 3-layer Titan perturbed by non-Keplerian orbit andatmospheric couplings. Planet. Space Sci., 93–94, 22–34, DOI:
10.1016/j.pss.2014.02.006.
Roberts J. H. (2015) The fluffy core of Enceladus. Icarus, 258, 54–66,
DOI: 10.1016/j.icarus.2015.05.033.
Roberts J. H. and Nimmo F. (2008) Tidal heating and the long-term
stability of a subsurface ocean on Enceladus. Icarus, 194(2), 675–
689, DOI: 10.1016/j.icarus.2007.11.010.
Smith B. A. et al. (1982) A new look at the Saturn system: The
Voyager 2 images. Science, 215(4532), 504–537.
Soucek O., Hron J., Behounkova M., and Cadek O. (2016) Effect of
the tiger stripes on the deformation of Saturn’s moon Enceladus.
Geophys. Res. Lett., 43, 7417–7423, DOI: 10.1002/2016GL069415.
Soucek O., Behounkova M., Cadek O., Tobie G., and Choblet G. (2017)
Tidal deformation of Enceladus’ ice shell with variable thickness
and Maxwell rheology. In 19th EGU General Assembly, EGU201716357.
Spahn F. et al. (2006) Cassini dust measurements at Enceladus andimplications for the origin of the E ring. Science, 311(5766), 1416–
1418, DOI: 10.1126/science.1121375.
Spencer J. R., Pearl J. C., Segura M., Flasar F. M., Mamoutkine A.,
Romani P., Buratti B. J., Hendrix A. R., Spilker L. J., and Lopes
R. M. C. (2006) Cassini encounters Enceladus: Background and thediscovery of a south polar hot spot. Science, 311(5766), 1401–1405,
DOI: 10.1126/science.1121661.
Spitale J. N. and Porco C. C. (2007) Association of the jets of Enceladus
with the warmest regions on its south-polar fractures. Nature,449(7163), 695–697, DOI: 10.1038/nature06217.
Squyres S. W., Reynolds R. T., Cassen P. M., and Peale S. J. (1983) The
evolution of Enceladus. Icarus, 53(2), 319–331, DOI: 10.1016/00191035(83)90152-5.
Tajeddine R., Rambaux N., Lainey V., Charnoz S., Richard A., RivoldiniA., and Noyelles B. (2014) Constraints on Mimas’ interior from
Cassini ISS libration measurements. Science, 346(6207), 322–324.
Tajeddine R., Soderlund K. M., Thomas P. C., Helfenstein P., Hedman
M. M., Burns J. A., and Schenk P. M. (2017) True polar wander
of Enceladus from topographic data. Icarus, 295, 46–60, DOI:
10.1016/j.icarus.2017.04.019.
Thomas P. C. (2010) Sizes, shapes, and derived properties of thesaturnian satellites after the Cassini nominal mission. Icarus, 208(1),395–401, DOI: 10.1016/j.icarus.2010.01.025.
Thomas P., Burns J., Helfenstein P., Squyres S., Veverka J., Porco C.,
Turtle E., McEwen A., Denk T., and Giese B. (2007) Shapes of thesaturnian icy satellites and their significance. Icarus, 190(2), 573–
584, DOI: 10.1016/j.icarus.2007.03.012.
Thomas P. C., Tajeddine R., Tiscareno M. S., Burns J. A., JosephJ., Loredo T. J., Helfenstein P., and Porco C. (2016) Enceladus’s

measured physical libration requires a global subsurface ocean.
Icarus, 264, 37–47, DOI: 10.1016/j.icarus.2015.08.037.
Tiscareno M. S., Thomas P. C., and Burns J. A. (2009) The rotation ofJanus and Epimetheus. Icarus, 204(1), 254–261, DOI: 10.1016/j.
icarus.2009.06.023.
Tobie G., Mocquet A., and Sotin C. (2005) Tidal dissipation within large
icy satellites: Applications to Europa and Titan. Icarus, 177(2),534–549, DOI: 10.1016/j.icarus.2005.04.006.
Tortora P., Zannoni M., Hemingway D., Nimmo F., Jacobson R. A., Iess
L., and Parisi M. (2016) Rhea gravity field and interior modelingfrom Cassini data analysis. Icarus, 264, 264–273, DOI: 10.1016/j.
icarus.2015.09.022.
Travis B. J. and Schubert G. (2015) Keeping Enceladus warm. Icarus,
250, 32–42, DOI:10.1016/j.icarus.2014.11.017.
Tricarico P. (2014) Multi-layer hydrostatic equilibrium of planets andsynchronous moons: Theory and application to Ceres and to solar
system moons. Astrophys. J., 782(2), 99, DOI: 10.1088/0004637X/782/2/99.
Tyler R. (2011) Tidal dynamical considerations constrain the state ofan ocean on Enceladus. Icarus, 211, 770–779, DOI: 10.1016/j.
icarus.2010.10.007.
van Hoolst T. (2015) Rotation of the terrestrial planets. In Treatise on
Geophysics, 2nd edition, Vol. 10: Physics of Terrestrial Planets andMoons (G. Schubert, ed.), pp. 121–151. Elsevier, Amsterdam.
van Hoolst T., Rambaux N., Karatekin Ö., and Baland R. M. (2009)
The effect of gravitational and pressure torques on Titan’s lengthof-day variations. Icarus, 200(1), 256–264, DOI: 10.1016/j.
icarus.2008.11.009.
van Hoolst T., Baland R.-M., and Trinh A. (2016) The diurnal librationand interior structure of Enceladus. Icarus, 277, 311–318, DOI:
10.1016/j.icarus.2016.05.025.
Vienne A. and Duriez L. (1995) Ephemerides of the major saturnian
satellites. Astron. Astrophys., 297, 588–605.
Waite J. H. et al. (2006) Cassini Ion and Neutral Mass Spectrometer:
Enceladus plume composition and structure. Science, 311(5766),
1419–1422, DOI: 10.1126/science.1121290.
Wieczorek M. A. (2015) Gravity and topography of the terrestrial
planets. In Treatise on Geophysics, 2nd edition, Vol. 10: Physics
of Terrestrial Planets and Moons (G. Schubert, ed.), pp. 153–193.
Elsevier, Amsterdam.
Wieczorek M. A. and Phillips R. J. (1998) Potential anomalies on asphere: Applications to the thickness of the lunar crust. J. Geophys.
Res., 103(E1), 1715–1724.
Willner K., Oberst J., Hussmann H., Giese B., Hoffmann H., Matz
K. D., Roatsch T., and Duxbury T. (2010) Phobos control point
network, rotation, and shape. Earth Planet. Sci. Lett., 294(3–4),
541–546, DOI: 10.1016/j.epsl.2009.07.033.
Yin A., Zuza A. V., and Pappalardo R. T. (2016) Mechanics of evenly
spaced strike-slip faults and its implications for the formation oftiger-stripe fractures on Saturn’s moon Enceladus. Icarus, 266, 204–
216, DOI: 10.1016/j.icarus.2015.10.027.
Zebker H. A., Iess L., Wall S. D., Lorenz R. D., Lunine J. I., andStiles B. W. (2012) Titan’s figure fatter, flatter than its gravity
field. Abstract P23F-01 presented at 2012 Fall Meeting, AGU, San
Francisco, California, 3–7 December.
Zharkov V. N. (2004) A theory of the equilibrium figureand gravitational field of the Galilean satellite Io: The
second approximation. Astron. Lett., 30(7), 496–507, DOI:
10.1134/1.1774402.
Zharkov V. N., Leontjev V. V., and Kozenko A. V. (1985) Models,
figures, and gravitational moments of the Galilean satellites ofJupiter and icy satellites of Saturn. Icarus, 61, 92–100, DOI:
10.1016/0019-1035(85)90157-5.