The Planetary Science Journal, 2:157 (14pp), 2021 August

https://doi.org/10.3847/PSJ/ac06d2© 2021. The Author(s). Published by the American Astronomical Society.

A Recipe for the Geophysical Exploration of Enceladus
Anton I. Ermakov1 , Ryan S. Park2 , Javier Roa2 , Julie C. Castillo-Rogez2 , James T. Keane2 , Francis Nimmo3 ,Edwin S. Kite4 , Christophe Sotin5 , T. Joseph W. Lazio2, Gregor Steinbrügge6 , Samuel M. Howell2 , Bruce G. Bills2 ,Douglas J. Hemingway7,8 , Vishnu Viswanathan9,10 , Gabriel Tobie11 , and Valery Lainey12
1
University of California, Berkeley, McCone Hall 307, CA 94720, USA; eai@berkeley.edu
2
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA
3
Department of Earth and Planetary Sciences, University of California, Santa Cruz, Santa Cruz, CA 95064, USA
4
Department of Geophysical Sciences, University of Chicago, Chicago, IL 60637, USA
5
Laboratoire de Planétologie et Géodynamique, Université de Nantes, F-44322 Nantes, France
6

Department of Geophysics, Stanford University, Stanford, CA 94305, USA
Earth & Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA
8
Maxar Technologies, 3875 Fabian Way, Palo Alto, CA 94303, USA
9NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA
10
University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250, USA
11
Laboratoire de Planétologie et Géodynamique, CNRS/Université de Nantes, France
12
IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, Univ. Lille, 77 Avenue Denfert-Rochereau, F-75014 Paris, France
Received 2021 February 10; revised 2021 May 25; accepted 2021 May 28; published 2021 August 10
7

Abstract
Orbital geophysical investigations of Enceladus are critical to understanding its energy budget. In this paper, we
identiﬁed key science questions for the geophysical exploration of Enceladus, answering which would support
future assessment of Enceladus’ astrobiological potential. Using a Bayesian framework, we explored how science
requirements map to measurement requirements. We performed mission simulations to study the sensitivity ofsingle-spacecraft and dual-spacecraft conﬁgurations to static gravity and tidal Love numbers of Enceladus. We ﬁndthat mapping Enceladus’ gravity ﬁeld, improving the accuracy of the physical libration amplitude, and measuring
Enceladus’ tidal response would provide critical constraints on the internal structure and establish a framework forassessing Enceladus’ long-term habitability. This kind of investigation could be carried out as part of a life search
mission at little additional resource requirements.
Uniﬁed Astronomy Thesaurus concepts: Planetary interior (1248); Gravitational ﬁelds (667); Saturnian satellites(1427); Space observatories (1543); Astrodynamics (76)
2. To offer recommendations for optimizing the futurecollection of geophysical data;
3. To suggest implementation options that would address
priority science questions as part of an Enceladus mission
concept within NASA’s medium New Frontiers–class ora large Flagship-class mission (National ResearchCouncil 2011).

1. Introduction
Enceladus—a cryovolcanically active and apparently habitable satellite in the Saturn system—challenges our understanding of geodynamical processes governing the evolution ofocean worlds (Waite et al. 2017; McKay et al. 2018). Startingin 2005, the Cassini spacecraft revealed active eruptions inEnceladus’ southern hemisphere (Porco et al. 2006), conﬁrmed
the presence of a deep, global ocean (Iess et al. 2014; Thomas
et al. 2016) with complex organic molecules (Postberg et al.
2018), and provided direct evidence for recent hydrothermal
activity (Hsu et al. 2015) that can produce redox disequilibria
necessary for supporting life (Waite et al. 2017; McKay et al.
2018). These characteristics make Enceladus a high-priority
target for astrobiology-driven exploration.
A focused geophysical investigation would underpin Enceladus’ astrobiological potential. Key to understanding how
Enceladus works is the spatial and temporal distribution of
the energy dissipated through tidal ﬂexing. This energy isnecessary for Enceladus’ ocean to persist over geologic
timescales (Roberts & Nimmo 2008).
The goals of this paper are:

The Cassini geophysical data, despite its limited resolution
and lack of global coverage, have yielded valuable constraints
on the interior state of Enceladus (e.g., McKinnon 2013, 2015;
Beuthe et al. 2016; Hemingway & Mittal 2019). The gravity
data along with shape models were used to constrain
Enceladus’ state of differentiation (Iess et al. 2014; McKinnon2015) and the long-wavelength variations in ice shell thickness(Beuthe et al. 2016; Čadek et al. 2019; Hemingway & Mittal
2019). The physical libration data revealed that the icy shell is
decoupled from the core by a liquid layer (Thomas et al. 2016).
The measured heat ﬂux in the South Polar Terrain provided
a lower bound on the amount of energy currently emitted
from Enceladus (Howett et al. 2011; Spencer et al. 2013).
The coherent spatiotemporal pattern of cryovolcanic activity(Hedman et al. 2013; Nimmo et al. 2014) provided constraints
on the rheology of and heat production within the cryovolcanically active region (Běhounková et al. 2015; Spitale et al.
2015; Kite & Rubin 2016) as well as provided hints of alonger-term variability (Ingersoll et al. 2020). The current stateof Enceladus’ geophysical data is summarized in Table 1 andillustrated in Figure 1.

1. To provide an overview of the science questions thatorbit-based geophysical data at Enceladus could answer;
Original content from this work may be used under the termsof the Creative Commons Attribution 4.0 licence. Any further
distribution of this work must maintain attribution to the author(s) and the titleof the work, journal citation and DOI.

1

The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

Table 1
Summary of the Orbit-based Enceladus Geophysical Data
Quantity

Current Knowledge

Information Provided

Gravity ﬁeld

Up to degree 3 ﬁeld measured by Cassini (Iess et al.
2014).

Constraint on the internal structure and compensation
mechanism in the ice shell.

Shape

Estimated up to degree 16 (Tajeddine et al. 2017).

Constraint on the internal structure and ice shell thickness
variations.

Gravity-topography admittancea

Degree 3 admittance derived (Iess et al. 2014).

Constraint on compensation state of topography, total shell
thickness, elastic shell thickness, and shell density and shellocean density contrast.

Obliquity

Cassini-derived upper bound of 0°. 05 (Giese 2014).
Upper theoretical bound 4 × 10−4° (Baland et al.
2016).

Enables determination of moment of inertia and whether
Enceladus is tidally damped into a Cassini state.

Physical libration amplitude

Measured amplitude at the orbital frequency of
0°. 120 ± 0°. 014 (Thomas et al. 2016) or
0°. 155 ± 0°. 014 (Nadezhdina et al. 2016).

Libration amplitude is strongly sensitive to the shell thickness. Observed large amplitude requires a decoupled ice shelland hence implies a global subsurface ocean.

Precession and nutation

Not currently measured. Precession rate estimated
2.6 rad yr−1 (Baland et al. 2016).

In combination with degree 2 gravity, enables moment ofinertia determination without the need for hydrostatic equilibrium assumption.

Tidally driven orbital migration

Measured from historical astrometric data (Lainey
et al. 2012).

Provides a constraint on the tidal dissipation factor Q of
Saturn.

Potential Love number k2; radial and lateral displacement Love numbers h2 and l2,
respectively.

Not currently measured.

Real parts are mostly sensitive to the thickness of the icy shelland its rigidity. Imaginary parts are sensitive to the viscosity
proﬁle and, thus, provide a constraint on total tidal dissipation
within Enceladus.

Radar subsurface mapping

Not currently achieved.

Enables independent shell thickness determination.

Magnetic induction

Not currently measured. The magnitude of the forcing
ﬁeld is ≈10 nT.

Sensitive to the subsurface ocean conductivity and thickness.

Thermal IR mapping

15.8 ± 3.1 GW from SPT (Howett et al. 2011); 4.2
GW from Tiger Stripes (Spencer et al. 2013).

Constraint on the total heat ﬂux.

Note.
aGravity-topography admittance is not an independent geophysical observable. It is derived as the ratio of the gravity amplitude to topography amplitude.

Tidal dissipation within Enceladus’ interior is a function of
its orbital characteristics (e.g., eccentricity and proximity toSaturn) and its internal structure (e.g., thickness of the ice
shell). The thermal and orbital evolutions of Enceladus are
therefore coupled (Meyer & Wisdom 2007, 2008; Běhounková
et al. 2012; Neveu & Rhoden 2019). That coupling has a strong
effect on the long-term evolution of the satellite. The presentday internal structure determines the instantaneous spatiallyvariable tidal dissipation rate. Depending on the efﬁciency of
the transport of heat from the interior to the surface, Enceladus
might or might not be in a thermal steady state. In that state, theheat produced within Enceladus (mostly from tidal dissipation)
would equal the heat Enceladus outputs to space. In addition,
tidal migration occurs, driven by dissipation within Saturn. Fast
tidal migration of the Saturnian moons determined fromastrometric observations indicates strong dissipation within
Saturn (Lainey et al. 2012). This observed fast migration
implies either that the major satellites are young or that the
migration rates have not been steady; the latter is predicted by
the resonance-locking mechanism with Saturn’s normal modes(Fuller et al. 2016; Lainey et al. 2020). Thus, by studying thecurrent internal structure and dissipation within Enceladus, we
can place constraints on its tidal migration history, which isinterlocked with the history of tidal heating.

We have identiﬁed the following interrelated priority science
questions that should be addressed by future geophysicalobservations of Enceladus:
1. What is the internal structure of Enceladus?
2. Where is the heat generated and how is it transported?
3. Is Enceladus currently in a steady state?
4. What are the feedbacks between volcanism and tectonics
that regulate Enceladus’ cryovolcanism?

2. Priority Science Questions
2.1. What is the Internal Structure of Enceladus?
Thanks to Cassini, we have a general understanding of
Enceladus’ internal structure (Hemingway et al. 2018). Gravity(Iess et al. 2014) and shape data (Nimmo et al. 2011; Tajeddine
et al. 2017) collected by Cassini, coupled with libration data
(Nadezhdina et al. 2016; Thomas et al. 2016), revealed that
Enceladus has a low-density, rocky core with an estimatedradius of 191–198 km, overlain by a global ocean with anestimated thickness of 30–39 km and an icy shell with anestimated average thickness of 19–24 km, but thinner at thesouth pole (4–12 km; Hemingway & Mittal 2019).
2

The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

Figure 1. Schematic illustration of Enceladus’ internal structure along with geophysical measurements and science priority questions identiﬁed in this paper.

topographic depression (Nimmo et al. 2011; Tajeddine et al.
2017) and hosts a set of prominent fractures, called Tiger
Stripes, that emit plumes of H2O, salt, organics, and othervolatiles indicating a connection to a reservoir of ocean-derived
material (Glein et al. 2018). Regional gravity and topography
mapping would allow the crustal structure in the South Polar
Terrain to be characterized. If gravity and topography data arecomplemented by radar sounding, an independent constraint on
spatial variations in shell thickness could be derived.
Mapping gravity and topography would also be useful inmore heavily cratered terrains in Enceladus’ northern hemisphere. The relaxed state of craters in that region indicateshigher heat ﬂuxes in the past as such heating episodes would
facilitate viscous ﬂow within the shell by reducing its viscosity(Bland et al. 2012). Higher-resolution gravity data would
enable the shell viscosity structure (Akiba et al. 2021) to beassessed and would offer an independent probe of the intense
past heating episode proposed by Bland et al. (2012).
If the precession of Enceladus’ pole could be measured, thatwould provide an independent constraint on the moment ofinertia of the shell. However, the very small amplitude of thedeﬂection (order of 1 m) would be challenging to measure,owing to Enceladus’ small obliquity (see Table 1).
Partitioning nonhydrostaticity between the shell and thecore—Separating the nonhydrostatic signal of the shell from thatof the core could be achieved by mapping the shell thickness
variations using a combination of topography, radar sounding,and gravity data. Radar sounding would enable tectonics andconvection zones to be identiﬁed, allowing the signals originating
within the ice shell to be separated. A direct detection of the ice–

However, these inferences rely on plausible, yet untested,assumptions. In particular, the core is assumed to be inhydrostatic equilibrium, which was challenged by McKinnon(2013) and Monteux et al. (2016). The gravity-topography
admittance, deﬁned as the wavelength-dependent ratio ofgravity amplitude to topography amplitude, was assumed
isotropic, which would be violated if there are signiﬁcantdeviations from spherical symmetry. Indeed, the inferred shell
thickness variations are on the order of the mean shell thickness
itself (Hemingway & Mittal 2019). If these assumptions arerelaxed, it becomes impossible to separate geophysical signalsthat originate in the icy shell (e.g., freezing/melting of theshell, convection, impact cratering, tectonics) from thoseoriginating in the core (e.g., core’s nonhydrostaticity and coretidal dissipation; Roberts 2015). Below, we summarize how afuture mission could sharpen our picture of Enceladus.
Icy shell structure—Currently, the mean ice shell thickness
is best constrained by libration data (Nadezhdina et al. 2016;
Thomas et al. 2016), while shell thickness variations are
deduced from gravity and topography (Hemingway et al. 2018;
Hemingway & Mittal 2019). Improving the estimate of thelibration amplitude is an easy objective for an Enceladus orbiter
mission. Continuous surface observations with a stereo imager
can bring down the libration amplitude uncertainty to 2 m (Parket al. 2020b)—a factor of 30 improvement over the current
uncertainty. The gravity and shape data determined to a higher
degree would allow to extend the shell thickness inversion tosmaller spatial scales, which would be especially valuable in
the South Polar Terrain. The South Polar Terrain is a broad
region to the south of 60° S. It is characterized by a 500 m deep3

The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

et al. (2014) and Běhounková et al. (2017) are exceptions.
Běhounková et al. (2017) found that degree 2 Love numbers ofdifferent order could vary by a factor of 2 due to the nonuniform
shell thickness and faults in the South Polar Terrain.

ocean interface would allow subtraction of the respective gravity
signal of the ice shell and inference of the properties of the core.
Improving the libration amplitude measurement at multiple
frequencies ranging from hours to tens of days could be used tobetter constrain the shell structure, speciﬁcally its moment ofinertia differences (such as the difference between the shell’sequatorial moments of inertia). Finally, measuring the gravitational signal of the core–shell misalignment (Buffett 1996) could
help eliminate the core–shell degeneracy.
Ocean thickness and composition—Ocean thickness anddensity are indirectly constrained by gravity and shape data bysatisfying the mass balance. The current uncertainty on theocean thickness as constrained by gravity, shape, and libration
data is a factor of ≈2 larger than that of the shell thickness(Hemingway & Mittal 2019). The gravity-topography admittance is sensitive to the density contrast between the ocean andthe shell. Thus, if the shell density is constrained from the highdegree admittance, the lower-degree admittance can help
constrain ocean density, which is mostly controlled by salinity.
A measurement of the magnetic induction response convolves ocean thickness and salinity. In the Jovian system,
magnetic induction has been used as a constraint on the ocean
thickness (Kivelson et al. 2000). However, unlike Jupiter’s,
Saturnʼs magnetic ﬁeld is axisymmetric. Thus, the changing
magnetic ﬁeld that Enceladus experiences is primarily due to its
eccentric orbit, resulting in a much smaller forcing ﬁeld
amplitude (≈10 nT). This signal is of roughly the same
magnitude as time-variable ﬁelds arising from the plumes(Dougherty et al. 2006). Thus, detecting an induction responseand disentangling it from the plume-induced variability in themagnetic ﬁeld would require long-term observations andcharacterization of the plume activity.
Rheology of the material—The tidal response of Enceladus isdescribed by its Love numbers knm, hnm, and lnm (e.g., Wahr et al.
2006), characterizing potential change, radial, and lateraldisplacements, respectively, where n is spherical harmonic degreeand m is order. For a spherically symmetric body, the Love
numbers are degenerate with respect to m. Thus, for such a body:
k20 = k21 = k22. Thus, we can drop the second subscript for
shorter notation if spherical symmetry is assumed. Separating thetidal signal corresponding to different orders within the same
degree would indicate deviations of the body’s mechanicalproperties from spherical symmetry. The Love numbers arecomplex quantities and quantify how Enceladus’ gravity ﬁeld andshape respond to time-varying tidal forces. The real parts of theLove numbers depend primarily on its shell thickness and shell
rigidity. The imaginary parts depend on the rheology of thematerial and quantify the lag in the tidal response with respect totidal forcing, which is related to the current total dissipation within
Enceladus. Independent measurement of both h2 (the degree 2
radial displacement Love number, omitting the order index for
simplicity) and k2 (the degree 2 gravitational potential Love
number) would reduce the correlation between the shell thickness
and rigidity (e.g., Wahr et al. 2006). In addition, libration
amplitude is sensitive to the shell rigidity in a way different fromtidal deformation (Van Hoolst et al. 2013). Thus, joint
measurement of the tidal response and libration would provide
independent constraints on the shell thickness and its rigidity.
One complicating factor for Enceladus is that the shell
thickness varies laterally, while almost all models (including
those presented in Section 3) have assumed a sphericallysymmetric shell when modeling the expected tidal response. A

2.2. Where Is the Heat Generated and How Is It Transported?
Tidal heating is thought to be the main heat source, preventing
quick freezing of the subsurface ocean (Roberts & Nimmo 2008).
The heat ﬂux of Enceladus has been measured by mapping out thesurface thermal emission using the Composite Infrared Spectrometer (CIRS) on board Cassini. The estimated value of the heat
ﬂux has varied over the years of Cassini data analysis by a factor
of several. All studies to date have focused only on the heat ﬂux
localized in the South Polar Terrain. Spencer et al. (2006)
estimated a heat ﬂux of 5.8 ± 1.9 GW. A subsequent study by
Howett et al. (2011) yielded a larger estimate of 15.8 ± 3.1 GW
using observations from a higher-wavelength CIRS detector.
Spencer et al. (2013) used higher spatial resolution data andderived a total emitted power of 4.2 GW localized at the Tiger
Stripes. The variability of these estimates derived from different
data sets indicates the challenges of deriving the global heat ﬂux
from ﬂyby observations of heterogeneous quality.
The conductive and advective heat ﬂuxes are approximately
equal in the South Polar Terrain (Kite & Rubin 2016). Advection
likely dominates within the Tiger Stripes, causing prominent
temperature anomalies. The heat ﬂux away from the South Polar
Terrain remains unconstrained. Current estimates on global heatﬂux assuming a conductive shell range from 25 to 40 GW(Hemingway & Mittal 2019). Radiogenic heating within the rocky
core can account for <0.3 GW (Hemingway & Mittal 2019),
hence it is a minor contributor to the observed heat ﬂow.
The local instantaneous tidal heat generation is determinedby the product of tidal stress and strain rate. The mutual
relation between these two quantities is, in turn, set by therheology of the material that is highly sensitive to temperature.
Heat is likely generated in the warmer ice at the base of theshell or in the core, especially if the core is unconsolidated(e.g., sandy/muddy; Roberts 2015). Souček et al. (2019) used a
3D model with a variable ice shell thickness including theSouth Polar Terrain faults and concluded that dissipation in thesolid ice cannot exceed 2.1 GW, implying that an additional
heat source is needed to explain the observed heat ﬂow.
Additional dissipation could occur within the liquid slots in theice (Kite & Rubin 2016). Over 10 GW can be generated bytidal friction inside the unconsolidated rocky core (Choblet
et al. 2017) and comparable amounts by tidal ﬂushing of water
through the porous core (Liao et al. 2020). Theoretical models
predict that the ocean heat production is likely negligible (Chenet al. 2014; Hay & Matsuyama 2017; Rekier et al. 2019;
Rovira-Navarro et al. 2019) because of Enceladus’ low
obliquity and thick ocean, although dissipation may beenhanced in narrow liquid-water conduits within the ice shell(Kite & Rubin 2016). However, more recently, Tyler (2020)
indicated that there could be resonant conﬁgurations, where
ocean tidal dissipation could create signiﬁcant heating.
If Enceladus is in thermal equilibrium, there exists a
relationship between the shell thickness and the conductive heat
ﬂux. Thus, the scientiﬁc requirement driving the uncertainty inthe mean shell thickness measurement can be tied to the desired
uncertainty in the heat ﬂux. It follows from Hemingway & Mittal(2019) that in order to determine the conductive heat ﬂux towithin 10%, the average shell thickness needs to be determined
4The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

to within 2 km. In the South Polar Terrain, where the ice shell
could be <5 km thick, radar observations may more easily reachthe ice–ocean interface and return the shell thickness with anaccuracy better than 100 m. This would constrain the regional
dynamics with implications for heat transport.
For a fully conductive and homogeneous shell, the thickness
derived from various techniques (i.e., radar returns, tidal Love
numbers, libration amplitude, or gravity-topography admittance) should be comparable. These complementary techniques
offer joint advantages as the deviations from homogeneity and/or from a conductive geothermal gradient can be obtained bytheir cross-analysis.
Distinguishing between convection and conduction is ofgreat interest. Conduction is likely dominant in the outer part ofthe shell. Convection, if it occurs within Enceladus’ icy shell,
would dominate heat transport in its deepest part, making the
conductive region thinner and with a steeper temperature
gradient and therefore greater conductive heat loss. Efﬁcientheat transport by convection might lead to quick ocean
freezing. Large-amplitude shell thickness variations inferredfrom the gravity and topography data (Nimmo et al. 2011;
Tajeddine et al. 2017; Hemingway & Mittal 2019) are at odds
with a global convective layer as it would lead to fast viscous
relaxation of shell thickness variations. However, a conclusive
determination of whether or not convection occurs within theicy shell would require joint analysis of multiple data sets.
A convective layer would affect the gravity-topography
admittance, leading to higher admittance values expected for
uncompensated topography (Watts 2001). Thus, higher-resolution gravity and topography data can help identify a convective
layer through their sensitivity to the viscosity proﬁle. For aconvective shell, warmer ice at shallower depth could lead tohigher radar signal attenuation, especially if the salinity of the ice
is high. A loss of signal resulting in no detected radar reﬂectionsat depth or weaker radar reﬂections from the warm and salty
bottom ice would indicate a saline, convective layer. In addition,a convecting shell would lead to an apparent mismatch between
the total ice shell thickness (derived from induction or gravitytopography admittance) and the thickness derived from the Love
numbers, because the latter is sensitive only to the elastic part ofthe shell and not the weak convecting part. If the convecting partof the ice shell is less saline, assuming that parts with higher
salinity would have undergone melting at some stage during
convection, the ice–ocean interface could be visible as a strong
reﬂection in the radar return. However, this direct signal could befurther covered up by an accretion zone at the ice–oceaninterface, which could produce ice of marine composition or
result in a mushy layer with high attenuation.
Estimating the elastic thickness would provide a criticalconstraint on the heat transport within the shell. Elastic thickness
can be derived either from the gravity-topography admittance(McGovern et al. 2002) or by mapping ﬂexural proﬁles oftectonic features that would require accurate topography knowledge (e.g., Giese et al. 2008). A ﬂexural proﬁle amplitude of
120 m at ≈10 km scale is predicted in the vicinity of Tiger
Stripes (Hemingway et al. 2020). Deriving spatial variations of
the elastic thickness and cross-correlating them with the heat
ﬂow mapped in the thermal IR emission could help validate thedissipation pattern within the shell. Measuring ﬂexural signals
would require regional topography knowledge to at least 10 m
vertical accuracy, which would require a dedicated stereo
mapping campaign. Stereo-derived digital terrain model can befurther improved and geodetically referenced to the center-ofmass frame by laser or radar altimeter data.
2.3. Is Enceladus Currently in a Steady State?
The long-term orbital evolution of Enceladus depends on
both the dissipation within Saturn and Enceladus. The
dissipation is described by the so-called quality factor Q,
which is proportional to the ratio of energy stored in tidal
motion to the energy dissipated over one tidal cycle. Q depends
on the frequency of tidal forcing (e.g., Wu 2005) and thus, canchange over the tidal migration history. In addition, Q may
evolve as the temperature of the body evolves due to secularcooling. Recent astrometric efforts indicate that Saturn is moredissipative than previously thought. Lainey et al. (2012) foundSaturn’s Q to be as small as ≈2000 at Enceladus’ tidal
frequency, permitting equilibrium tidal dissipation within
Enceladus to be as much as ≈25 GW (Meyer & Wisdom 2007).
The feedback between orbital and thermal evolution leads totwo distinct kinds of steady state. First, in the thermal steadystate, the present-day heat production is equal to the present-day
heat loss. Second, in an orbital steady state, damping of
eccentricity, e, due to dissipation within Enceladus is balancedby the pumping of its eccentricity by the 2:1 resonance between
Dione and Enceladus (i.e., de/dt = 0; Meyer & Wisdom 2007).
If not in a steady state, Enceladus could exhibit a periodic
behavior (Ojakangas & Stevenson 1986), in which energy isstored at one time and released at a later time (perhaps resulting
in cyclical variations of ocean thickness), or it might be in a net
freezing or melting state. Indeed, topographic evidence forviscous relaxation of impact craters in ancient terrains (Bland
et al. 2012) and mapping of fault patterns away from the South
Polar Terrain (Patterson et al. 2018) both show that the location
and/or intensity of crustal heat ﬂow and associated tectonic
activity has changed over geologic time. Geophysical investigation of these tectonic anomalies may provide clues as to nonsteady-state behavior in Enceladus’ past: for example, theprocesses that drove the onset of hyperactive resurfacing near theSouth Pole, but not the North Pole (e.g., Kang & Flierl 2020).
Both kinds of steady states critically depend on the tidal phaselag, which is characterized by the imaginary part of the potential
Love number, Im (k 2 ). Larger values of Im (k 2 ) correspond tolarger internal tidal dissipation (see Equation (5)) and faster tidal
damping of eccentricity. Im (k 2 ) can be derived in two ways.
First, a tidal phase lag would affect the orbit of the spacecraft:the measured response would lag that expected from theperturbing potential. Thus, it can be derived by radio-trackingof the spacecraft in the same way the gravity ﬁeld is derived.
Second, the tidal phase lag affects the orbit of Enceladus,
causing the damping of its eccentricity. Thus, it could be derived
by determining an accurate ephemeris model using ground- orspacecraft-based data (e.g., radio ranging or astrometry).
Pursuing both ways of deriving the tidal phase lag would
provide a robustness check. In addition, the change of
eccentricity should be accompanied by the corresponding
evolution of the Enceladus–Dione resonance libration angle,
the determination of which would require precise ephemerides ofboth Enceladus and Dione.
In conclusion, long-baseline continuous ground-based astrometry at the kilometer level and future radio ranging and highaccuracy astrometry in the Saturnian system (20 yr after
Cassini) would be needed to reveal whether Enceladus is in
the orbital steady state. Thermal IR mapping of Enceladus,
5

The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

including the conductive heat ﬂux away from the Tiger Stripes,
is needed to assess if Enceladus is in the thermal steady state.

Table 2
Prior Probability Distribution of Enceladus’ Internal Structure Parameters
Minimum
Value

Parameter

2.4. What Are the Feedbacks between Volcanism and Tectonics
That Regulate Enceladus’ Cryovolcanism?
Enceladus cryovolcanism leads to signiﬁcant mass loss fromthe moon (150–350 kg s−1 or 4%–10% of Enceladus’ mass per
Gyr; Hansen et al. 2006) and also offers the opportunity tosample material from Enceladus’ ocean. However, the causesof the tectonic features of the volcanically active South Polar
Terrain are not well understood (e.g., Yin & Pappalardo 2015;
Hemingway et al. 2020). The resurfacing mechanisms as wellas rates of resurfacing and surface–interior exchange, which areimportant for sustaining habitability (e.g., Soderlund et al.
2020), are not well constrained (Bland et al. 2015; Spencer
et al. 2018). As a result, we do not know whether the presentday rate of cryovolcanic activity is representative of the historyof the South Polar Terrain (O’Neill & Nimmo 2010).
Enceladus operates in a regime intriguingly different fromthe other worlds known to have active volcanism such as Earthand Io. Earth’s oceanic lithosphere is dominantly cooled byconduction (ratio of advected to conducted heat ≈0.1). On Io,
volcanism drives the cooling and, thus, the tectonics (ratio of
advected to conducted heat ≈10). By contrast, on Enceladus,
the ratio of volcanic heat to conducted heat appears to be ≈1,
implying that the South Polar Terrain evolution cannot be
understood without considering both processes. The strong
potential coupling between volcanism and tectonics highlights
the importance of a high-resolution gravity and topography
mapping within the South Polar Terrain as well as higher
spatial resolution mapping of heat ﬂow. Estimates of the elastic
thickness from the ﬂexural proﬁles collected in the vicinity of
the cracks might be signiﬁcantly different from the estimatesusing localized gravity-topography admittance as a signiﬁcantadvection of heat occurs in the fractures. Better constraints on
the spatiotemporal pattern of eruptions (Nimmo et al. 2014;
Spitale et al. 2015) might constrain the heat production within
the Tiger Stripes and the plumbing system, as well as tectonic
mechanisms operating within the South Polar Terrain (Kite &Rubin 2016). High-resolution mapping of how surface thermal
emission falls off with distance from the Tiger Stripes would
provide additional constraints on the thermal structure of theplumbing system (e.g., Abramov & Spencer 2009).
Radar sounding would provide a direct way to investigate
tectonic activity. Faulting processes can generate dielectric
discontinuities in the ice shell, which would then act as reﬂectors
for electromagnetic waves and manifest themselves in the radar
return. Vertical fractures could be identiﬁed by point scattering,usually leaving characteristic hyperbolas in the (nonfocused)radar return. Furthermore, compositional variations associatedwith tectonic activity, injection of water, or even variations in
crystalline fabric can be detected in radar sounding data. Due tothe complex surface topography of Enceladus, the radarsounding data would need to be complemented by highresolution topography data to mitigate surface clutter.
Higher-resolution gravity, shape, and heat ﬂux, as well asradar sounding and data within the South Polar Terrain, would
open the prospect of studying a new kind of tectonics as richand unusual by terrestrial standards as plate tectonics when ﬁrst
understood in the 1960s and 1970s. High-resolution gravity(50 km resolution, or degree ≈30, to fully resolve the South
Polar Terrain) in combination with stereo digital elevation

Maximum
Value

Distribution
Type

Shell shear modulus [GPa]110log-uniform

Shell Poissonʼs ratio

0.25

0.45

uniform

Shell viscosity [Pa s]

14

19

10

10

Shell density [kg m ]

600

1100

uniform

Shell thickness [km]240

uniform

Ocean compressibility [GPa]23log-uniform

Ocean density[kg m−3]

920

1300

uniform

Ocean thickness [km]280

uniform

Core shear modulus [GPa]

10

100

log-uniform

Core Poissonʼs ratio

0.25

0.45

uniform

Core viscosity [Pa s]

20

22

−3

−3

Core density [kg m ]Core radius [km]log-uniform

1010log-uniform

2000

3000

uniform

2

252

uniform

Note. In addition to these bounds, the densities were constrained to increase
monotonically with depth. The thickness of the three layers were constrained to
add up to the outer radius of Enceladus.

models could be particularly valuable. Crucial constraints on
the rates of both magmatic and tectonic activity, as well as links
between surface and subsurface activity, might be provided
by Interferometric Synthetic Aperture Radar (InSAR) measurements—analogous to what has been done successfully for
Earth such as mapping the ice ﬂow velocities over Antarctica(Rignot et al. 2011) or mapping active volcanism-induced
deformations (see Segall 2010 and references therein).
3. Development of Measurement Requirements
Geophysical measurement requirements can be derived in
multiple ways as several combinations of measurements canyield identical accuracy for a recovered parameter. For
example, gravity and radar measurements can both yield shell
thickness estimates. In addition, the measurement requirements
might be dependent on the (yet unknown) value of therecovered parameter. For example, the libration amplitude is aninverse function of shell thickness.
In order to develop traceable measurement requirements, weused the Markov Chain Monte Carlo (MCMC) approach todevelop a framework for connecting the science requirements
to the measurement requirements (Matsuyama et al. 2016). We
use central values for the internal structure parameters (such as
layer thicknesses and densities) from Hemingway & Mittal(2019) as the truth values. In addition, we assumed a priorprobability distribution for all model parameters, including thecurrently unconstrained (due to a lack of Love number
measurements) viscoelastic moduli. The form of the priorprobability distribution is summarized in Table 2.
6

The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

Figure 2. Corner plot visualizing the posterior distribution of internal structure model parameters. Histograms for layer thicknesses (hi) and shear modulus of the icy
shell (μshell) are shown. Note that μshell is shown on a logarithmic scale. μshell is constrained only when Love numbers are measured. The contours show the 95th
percentile. The black contour corresponds to the current state of knowledge, while the colored contours show how the posterior distribution changes depending on the
additional measurements included in the inversion. Zl stands for the gravity topography admittance at degree l. The table at the top shows the measurement
uncertainties for the corresponding contours. The inﬁnity sign in the table indicates that the quantity was not included in the inversion.

To generate synthetic observations, we compute a multilayerhydrostatic equilibrium model using Tricarico (2014). Viscoelastic moduli and densities are assumed constant within eachlayer. This gives us the hydrostatic shape and gravity sphericalharmonic coefﬁcients. Degree 2 complex Love numbers arecomputed in a way similar to Kamata et al. (2015). Enceladus
is assumed to have a small nonhydrostatic topography, which iscompensated with an Airy isostasy mechanism (Watts 2001).
This yields modeled gravity-topography admittance. The
correlation between gravity and topography is assumed to beunity. The libration amplitude at the orbital frequency iscomputed for rigid shells (Van Hoolst et al. 2008), which is a
simpliﬁcation of our model. As the libration measurement
accuracy is improved, it becomes sensitive to the rigidity of theshell. Čadek et al. (2016) estimated that varying the shell
rigidity between 1 and 5 GPa leads to a libration amplitudedifference of 50 m, which is just below current observational
uncertainty. The variation of libration amplitude due to theshell viscosity proﬁle is yet smaller: <20 m for bottom shellviscosity ranging from 1012 to 1015 Pa s (Van Hoolst et al.

2016). More details on the MCMC geophysical inversion are
given in the Appendix.
In assuming the measurement uncertainties, we impose that
future measurements of physical libration and gravity
coefﬁcients will not be worse than the current ones (Iesset al. 2014; Thomas et al. 2016). We explore the combined
sensitivity of physical libration amplitude at the orbital
period; tidal Love numbers k2, h2, and l2; and the gravitytopography admittance spectrum Zn, to the internal structure
parameters. For each set of measurement errors, we derive theposterior distribution of internal structure model parameters.
In order to visualize the multidimensional posterior distribution, we marginalize it over the parameters of interest (e.g.,shell thickness or shell density).
We have explored several combinations of measurement
uncertainties and mapped them with MCMC into the internal
structure parameter posterior distributions. Figure 2 shows theposterior distributions for selected internal structure parametersin the form of a corner plot (Foreman-Mackey 2016). The top
boxes show the 1D histograms of the parameters of interest. In
7

The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

Figure 3. Corner plot visualizing the posterior distribution of Enceladus’ real parts of Love numbers as constrained by the gravity (Iess et al. 2014), shape (Tajeddine
et al. 2017), and libration data (Thomas et al. 2016). The darker colors in the 2D histograms indicate higher probability values. The 2.5th–97.5th percentile ranges are
given at the top.

addition, we show 95% conﬁdence regions for selectedparameter pairs. The black contours in Figure 2 represent thecurrent state of knowledge. The other (colored) contours show
how the conﬁdence regions shrink as more data are included in
the MCMC inversion.
We ﬁnd that improving the accuracy of admittance (mostly)
and libration (less so) can decrease the shell thickness
uncertainty. Improving the accuracy of the gravity-topography
admittance to the level of 1 mGal km−1 up to degree 10 andlibration amplitude to an accuracy of 6 m can decrease that
uncertainty down to 2 km, which corresponds to a 10%
uncertainty in the conductive heat ﬂux. The Love numbers
provide sensitivity to viscoelastic moduli. If k2 is measured to10−2 or better, it provides a constraint on shell rigidity andimproves the accuracy of the shell thickness determination. If
k2 is measured to 10−3, it becomes sensitive to the shellviscosity.

3.1. Estimating the Amplitude of Tidal Deformations
Having estimated Enceladus’ Love numbers using an
MCMC inversion, we can statistically assess the magnitudeof surface tidal deformation. The posterior distribution for theLove numbers is shown in Figure 3. The derived conﬁdence
intervals for the Love numbers allow an estimation of the tidal
displacement ranges and surface gravity changes to bemeasured by future spacecraft missions. The tidal surface
displacement is proportional to the displacement Love numbers
h2 and l2. The maximum range of tidal displacement fromdegree 2 tides is given by Park et al. (2020b) asmax (DR) = 12

2 h2 R 2w 2e,7gl R 2w 2emax (DN ) = 9 2
,
g

8

(1 )
(2 )

The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

Figure 4. Radial, northerly, and easterly tidal displacement ranges as well as the surface gravity range of Enceladus given the median values of the Love numbers asconstrained by the gravity (Iess et al. 2014), shape (Tajeddine et al. 2017), and libration data (Thomas et al. 2016).

l R 2w 2emax (DE ) = 12 2
,
g

(3 )

for the radial, northerly, and easterly displacements, respectively.
Here, R is the satellite’s radius, e is its orbital eccentricity, ɷ isthe orbital frequency, and g is the surface gravity. In addition,
the maximum surface gravity variation range is given bymax (Dg) = 242⎛3
1 - k 2 + h2 ⎞ Rw 2e.
7⎝
2
⎠

(4 )

Taking the median values for the Love numbers from Figure 3(k2 = 0.0167, h2 = 0.0429, l2 = 0.0096), we derive the maximum
range of tidal displacement of 2.0, 0.6, and 0.9 m in the radial,
northerly, and easterly directions. The map of tidal displacement
ranges for these median values of Love numbers is shown in
Figure 4. In addition, Figure 4 shows the range of surface gravity
changes, which has the same pattern as the vertical tidal
displacement range. It can be seen that in our spherically symmetricmodel, radial tidal displacements are maximized at the equator.
The maximum radial tidal deformation range is 0.5–6.5 m(2.5th–95th conﬁdence interval) in the equatorial region but
can be ampliﬁed in the vicinity of the Tiger Stripes(Běhounková et al. 2017; Marusiak et al. 2021). We note thattidal displacements are expected to be much smaller comparedto the libration amplitude at the orbital frequency (≈530 m at
the equator, Thomas et al. 2016). Thus, an accurate libration
model would be required to tease out tidal deformations. The
combined measurement of the gravity-topography admittance,
libration amplitude, and tidal deformation can effectively
reduce the shell thickness uncertainty. Finally, measuring
obliquity and precession would require a submeter accuracy on
the Enceladus orientation.

Figure 5. Tidal dissipation of Enceladus as a function of the imaginary part of
the Love number.

Given the observed heat ﬂux (Howett et al. 2011), we canderive an estimate of the imaginary part of the potential Love
number using the tidal dissipation formula (Peale et al. 1979):
21 GMSaturn REnceladus weE˙ = Im (k 2)
,
2
a6
2

5

2

(5 )

where G is the universal gravitational constant, and a is thesemimajor axis of Enceladus’ orbit. Thus, the measurement
requirement on Im (k 2 ) should be tied to the requirement on
the heat ﬂux measurement. Figure 5 shows the total tidal
9

The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

dissipation as a function of Im (k 2 ), from which follows thatmeasuring Im (k 2 ) to 10−3–10−2 would be required for
detecting the tidal lag and constraining the total dissipation,
assuming Enceladus is currently in thermal equilibrium.
Finally, the recovery of h2 and l2 with the same accuracy as fork2 can help mitigate the ambiguity of the ice rheology thatarises when measuring k2 only (Wahr et al. 2006).

spacecraft remaining in orbit without escaping or impacting
Enceladus. In practice, a certain orbit may be considered stableas long as its lifetime, although ﬁnite, is long enough tocomplete the scientiﬁc exploration phase according to themission requirements.
Tidal perturbations due to Enceladus’ tidal deformation aresmall compared to the perturbation from Saturn—evencompared to the third-order terms of the static gravity ﬁeld—
and do not pose a risk to orbital stability. For this reason, it isrelevant to ﬁnd the orbit conﬁgurations that maximize the tidal
signature when designing a mission to recover the Love
numbers.
Low-altitude, near-circular orbits are a versatile option
suitable for achieving different scientiﬁc goals, including
gravity recovery and surface mapping. However, such orbits
are only stable at relatively low inclinations and low altitudes
over Enceladus. The Lidov–Kozai (LK) mechanism (Kozai
1962; Lidov 1962) produces a long-period exchange between
inclination and eccentricity, effectively increasing the eccentricity of the orbit while reducing its inclination and vice versa.
Assuming that the radius of the spacecraft orbit is small
compared to the radius of the orbit of Enceladus, a ﬁrst-order
expansion of the three-body potential leads to the quadrupole
approximation, which predicts the maximum eccentricity
during an LK cycle to be

4. Mission Design Considerations
4.1. Gravity and Tides Sensitivity
We have studied the sensitivity of single spacecraft with
Earth-based Doppler tracking and dual spacecraft withinterspacecraft tracking (i.e., GRACE/GRAIL-like conﬁguration, Tapley et al. 2004; Zuber et al. 2013) to map the static andtemporally variable gravity ﬁeld of Enceladus. Polar orbits with
global coverage are preferred for global static gravity mapping,surface imaging, and shape reconstruction. However, such
orbits at Enceladus are unstable and other special types of highinclination orbits need to be considered (Russell & Lara 2009;
Massarweh & Cappuccio 2020). High-inclination orbits arealso of interest for gravity mapping of the South Polar Terrain
of Enceladus that lies below 60° S.
Tidal deformation causes the gravity ﬁeld of Enceladus tovary in time, which affects the motion of spacecraft around it.
When operating a formation of multiple spacecraft aroundEnceladus, the tidal investigation would rely both on the effectof tides on each spacecraft but also on the differential effect ofthe perturbation, which affects the spacecraft relative motion.
The effect of tides on the motion of a spacecraft can be split
into two components:

emax =

1-

5
cos2 i ,3(6 )where i is the orbital inclination. This approximation is only
valid for cos2 i < 3 5, or for inclinations between the Kozai
angles imin = 39 . 2 and imax = 140 . 8. The eccentricity of
initially circular orbits with inclinations between these criticalvalues would grow due to the LK mechanism, effectively
lowering its periapsis eventually leading to an impact. Given
the semimajor axis of the orbit, Equation (6) can be used toestimate the inclination at which the periapsis reaches thesurface of Enceladus. However, Lara et al. (2010) proved thatthe quadrupole expansion is not sufﬁcient to accurately modelthe dynamics at Enceladus and found that a higher-order
expansion predicts stable near-circular orbits at inclinations ashigh as 50°.
The maximum altitude that low-inclination near-circularorbits can reach while remaining stable is driven by the Coriolis
asymmetry (Innanen 1979). When modeling orbital motion in
the Saturn–Enceladus synodic frame, the Coriolis effectopposes the gravitational attraction from Enceladus for
prograde motion whereas it supplements Enceladus’ gravity
for retrograde motion. As a result, the maximum stable orbital
radius of prograde orbits is only half of that for retrograde
orbits. Stable prograde circular equatorial orbits are stablebelow altitudes of 200 km, whereas their retrograde counterparts are stable up to 700 km above the surface.
An alternative design approach is to seek periodic orbits in
the Saturn–Enceladus circular restricted three-body problem(CR3BP), which may present larger eccentricities and alsoreach higher inclinations (Russell & Lara 2009). The different
families of periodic orbits in the CR3BP provide mission
planners with a wide variety of options for pursuing speciﬁc
objectives. For example, the low-altitude near-circular equatorial orbits already discussed are members of the family of distantprograde orbits (DPOs) and distant retrograde orbits (DROs),

1. Direct (short-term): tides cause short-period variations in
the orbital elements of a single spacecraft or periodic
changes in the range rate between two spacecraft over
one orbital revolution.
2. Indirect (long-term): tides can cause long-term variationsin the orbital elements or long-term changes in the range
rate between two spacecraft. The indirect effect may be
ampliﬁed by resonances between the spacecraft’s meanmotion and tidal harmonics.
The amplitude of the tidal deformation and tidal gravity
perturbation is at the maximum at the equator (Figure 4). For
this reason, equatorial orbits are in general better for opticallyobserving tidal deformation and detecting tidally induced
gravity variations as they would maximize the observablesignal. In addition, displacement due to forced libration is also
maximized at the equator. Deviations from spherical symmetry
would cause Love numbers to be nondegenerate with respect toorder m, with degree 2 Love numbers differing by up to a factorof 2 (Běhounková et al. 2017). In order to separate thenondegenerate Love numbers, a nonequatorial orbit would beneeded. Thus, the choice of an optimal orbit should depend on
the expected internal structure and should be driven by ahypothesis to maximize the sensitivity to the phenomenon inquestion.
4.2. Orbital Stability
The stability of orbits around Enceladus is driven by the tidal
force due to Saturn. Orbital stability in this context refers to a
10

The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.
−9−1

and 5 × 10 km s over a 10 s integration time (Asmar et al.
2005; Bills & Ermakov 2019). The impact of the frequency
used for spacecraft Doppler tracking depends on the S/N of thetwo-way radio link, which would depend on the spacecraft
thermal noise, transmitter power, received signal noise, antenna
type, etc. For the single-spacecraft tracking scenario, the radio
wave will go through Earth media and Earth-to-Saturn distanceencountering spatially and temporally variable solar plasma.
On the other hand, for dual-spacecraft tracking, the radio wave
will only have to travel a few hundreds of kilometers. In thelatter case, the dual-spacecraft tracking data quality would
mainly be affected by the spacecraft thermal noise, which
would be substantially better than in the single-spacecraft case,i.e., increasing accuracy by at least a factor of 10, which would
yield at least an order of magnitude better result than the singlespacecraft case.
We assumed an X-band tracking accuracy for the singlespacecraft case, which is typical for a deep-space mission. We
only considered X-band for single spacecraft because the Kaband uplink capability is currently available only at theGoldstone DSN station. In addition, the Ka-band uplinkcapability is costly to maintain as it requires a water vapor
radiometer to get the full accuracy. If a Ka-band tracking of asingle spacecraft is assumed, it can lead to a factor of 4–10
measurement accuracy improvement.
We assumed a stable orbit around Enceladus with 60°
inclination and periapsis and apoapsis altitudes of 150 km and200 km, respectively. Both poles can be covered with ﬂybys
prior to being placed into this stable orbit, depending on thescience requirements, which was not modeled in our simulations. We have simulated how accurately we can recover
Enceladus’ geophysical parameters, such as Love numbers,
gravity ﬁeld, and libration amplitude. The resulting radial
gravity acceleration error spectra are shown in Figure 6. In
addition, Figure 6 shows the expected amplitude of the gravity
signal as well as the spatial scales of several geologic featuresof interest.
Assuming a single orbiter with two-way X-band tracking
capability and assuming the Doppler data accuracy of
10−7 km s−1 at a 60 s count time with a 28 day data collection
duration, a gravity ﬁeld up to degree ≈10 can be recovered andk2 can be recovered with an accuracy of ≈10−2. A highresolution context imager with accurate pointing knowledge
and image-motion-compensation capability, such as theAdvanced Pointing Imaging Camera (APIC; Park et al.
2020b), can measure h2 and l2 with accuracies ≈0.01
and ≈0.002, respectively. This corresponds to a relative
uncertainty of the displacement Love numbers of 20%–25%.
The libration amplitude can be recovered with an accuracy of<1 m at the equator (a relative uncertainty of <0.2%). The
accuracy of h2, l2, and libration amplitude can be improved
with the optimal distribution of crossover points (Park et al.
2015) or lower orbit altitude.
Recovering h2 can further be achieved by laser altimetry(Steinbrügge et al. 2015) or by radar in combination with stereo
imaging at crossover points (Steinbrügge et al. 2018).
Equivalent measurements are already planned with the
Ganymede Laser Altimeter (GALA; Hussmann et al. 2019)
and the Europa Clipper’s instrument Radar for Europa: Ocean
to Near-Surface (REASON; Blankenship et al. 2018). State-ofthe-art laser altimeters can achieve vertical resolutions of 10 cm(Araki et al. 2019). However, the error budget for therespectively. The continuation of the family of periodic DPOs
yields orbits with different geometries but these are also
unstable. Russell & Lara (2009) noted that northern halo orbits,
which are eccentric and highly inclined, provide unique
opportunities for exploring the South Polar Terrain andobserving the Tiger Stripes, a conﬁguration similar to Dawn’s
Second Extended Mission (Park et al. 2020a). A subset of these
halo orbits is stable with periapsis altitudes ranging from 30 km
all the way down to the surface. Although stable in thesimpliﬁed CR3BP model, the actual orbit lifetime of this set of
halo orbits is less than one week. Davis et al. (2018) presented
additional periodic orbits with interesting opportunities toobserve the South Polar Terrain. Although such highly inclined
orbits tend to be unstable in the full-ephemeris model, Daviset al. (2018) demonstrated that halo orbits ﬂying over the South
Polar Terrain at altitudes lower than 200 km can be controlledwith an estimated ΔV cost of 20 m s−1 per month. Nuclear
electric propulsion can signiﬁcantly increase the maneuvering
capabilities for orbit transfer and control (Casani et al. 2020),
potentially enabling a wider range of science orbits.
A mission consisting of more than one spacecraft requires
not only the orbit of each spacecraft to be stable, but also
ensures that the relative formation is preserved. The potential offormation ﬂying concepts for deep-space exploration has been
discussed in the past, including design and control considerations (Gurﬁl et al. 2003; Howell & Marchand 2005). Examples
of proposed applications include close-proximity exploration ofsmall bodies (Baresi et al. 2016; Lippe & D’Amico 2020) andformations along generic periodic orbits in the CR3BP (Gurﬁl& Kasdin 2004). Still, formation ﬂying in close proximity to
tidally locked satellites remains an understudied area.
For the case of Enceladus, the perturbation from Saturn’sand Enceladus’ nonuniform gravity can produce strong
differential accelerations on the spacecraft and destabilize theformation. Passively mitigating the differential effect of theperturbations during the design stage is crucial for formulating
an efﬁcient control plan that minimizes the overall cost of
station-keeping. Active control may still be required tomaintain the formation within operational requirements. Precise
autonomous guidance, navigation, and control are required toensure that the spacecraft can react to and counteract
perturbations efﬁciently. Constant maneuvering can interferewith the continuity of the radio-tracking data. Thus, stationkeeping maneuvers should be kept at a minimum whenpossible.
4.3. Mission Simulations
Based on the aforementioned considerations, we conducted aseries of mission simulations and performed detailed covariance analyses for two mission conﬁgurations: a single orbiterwith radiometric tracking to Earth and dual spacecraft (GRAILlike conﬁguration) with intersatellite range-rate measurements.
The covariance analysis is based on a least-squares principle
and is a powerful tool for assessing the expected uncertainties
in the estimated parameters (Park et al. 2011, 2012, 2015).
For simplicity, we assumed that both single and dual
spacecraft are tracked continuously. In reality, the tracking of
single spacecraft is limited by the Deep Space Network (DSN)
availability. On interplanetary ranges, the S band (2.3 GHz), X
band (8.4 GHz), and Ka band (32 GHz) have been used forspacecraft radio-tracking Doppler measurements. These bands
respectively yield ranging accuracies of ≈5 × 10−7, 5 × 10−8,11The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

Figure 6. Summary of the gravity anomaly signals due to various geologic landmarks. The vertical axis shows the rms magnitude of the gravity signal. The top
horizontal axis shows the spherical harmonic degree n. The bottom horizontal axis shows the corresponding spatial wavelength, λ ≈ 2πR/n. The light blue region
indicates the expected range of the gravity signal. The upper bound is given by a power law ﬁtted to an uncompensated gravity-from-shape using the shape model byTajeddine et al. (2017) with a shell density of 920 kg m−3. The lower bound corresponds to gravity from compensated topography using the same shape model with ashell density of 600 kg m−3. The large-scale depressions refer to the chain of basins identiﬁed by Tajeddine et al. (2017).
Table 3
A Summary of the Potential Love Numbers and Gravity Recovery
Range-rate Accuracy (km s−1)
Single spacecraft

−7

10

−8

10
Dual spacecraft

−9

10
10

−10

k20

k21
−3

9 × 10−4

9 × 10−4

2 × 10−5

2 × 10k22

5 × 10−1

2 × 10−2

3 × 10−3

4 × 10−4

Global Degree Strength

4 × 10−3

2 × 10−4

2 × 10−5

3 × 10

−6

9
17
23
30

Note. Standard deviations of the degree 2 Love numbers are shown. The degree strength refers to the spherical harmonic degree at which the gravity signal is equal to
the noise in the gravity data. Thus, degree strength deﬁnes the global resolution of the gravity model.

than the inherent range resolution can be achieved (Steinbrügge
et al. 2018).
A dual-spacecraft architecture with intersatellite Ka-bandtracking (i.e., GRAIL-like scenario) was also considered for thesame orbit conﬁguration described above. We assumed rangerate accuracies of 10−8, 10−9, and 10−10 km s−1. These or
better (down to 3 × 10−11 km s−1) ranging accuracies wereachieved by the GRAIL mission (Konopliv et al. 2013). Suchmeasurement conﬁguration can signiﬁcantly improve theaccuracy of the k2 determination yielding a k2 error down
to ≈2 × 10−5 (or relative accuracy of ≈0.1%; Table 3). The
gravity ﬁeld can be determined to degrees 20–30 (see Figure 6),
depending on the ranging accuracy. This level of accuracy
would allow a deﬁnitive determination of the tidal phase lagand total tidal dissipation within Enceladus. The coreviscoelastic moduli and ocean compressibility would remain
virtually unconstrained for the studied mission conﬁgurations,although a lower bound on the core viscosity could be retrieved
from the Love numbers if they are measured to an accuracyof 0.1%.

measurement of radial tides is dominated by interpolation
errors at crossover points and orbit determination uncertainties
of the spacecraft. Having an orbiter with an accurate radio
science experiment is therefore beneﬁcial for the h2 measurement as well. Typical accuracies for such a conﬁguration would
be on the order of 1 m. The accuracy of the h2 inversion thendepends on the number of orbits, norb, a spacecraft completes
within its mission lifetime. Because the number of crossover
points increases with n2orb, millions of crossover points can becollected over mission lifetimes of a few months.
Using a radar sounder instead of a dedicated altimeterrequires an additional stereo camera to mitigate the ambiguity
induced by surface clutter. The radar altimetric accuracy
depends on the radar bandwidth, which is limited by half thecenter frequency. Ice is particularly transparent to frequencies
between 1 and 300 MHz (Blankenship et al. 2009). To achievea 3 m vertical resolution, a 50 MHz bandwidth is required,
which is larger than the bandwidth used by typical radarsounders. However, by using stereo imaging in combination
with radar sounders, resolutions up to a factor of 4–5 better
12

The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

Advances in small radios based on the IRIS transponder on
the Mars CubeSat One mission and Artemis-1 CubeSats (e.g.,
software programmable universal Space Transponder-Lite, or
UST-Lite ≈1 kg, ≈15 W, two-way Allan deviation of 10−14 at
1000 s; Pugh et al. 2017) would enable a GRAIL-like scenario
at little resource cost to the primary mission but with lessaccurate satellite-to-satellite ranging accuracy. The CubeSat
could be carried to the destination and network with themothership via a deep-space deployer providing power/
thermal during cruise. This CubeSat would be enabled by
new technologies in thermal management and propulsion
subsystems.

Markov chain to generate the positions at the following step.
The walker positions are updated based on the likelihood
function. The likelihood function tells how well a model
reproduces the observations. We use the likelihood function in
the following form:
1
log L µ - (X - Y )T S-1(X - Y ) ,2(A1)

where X is the vector of observations and Y is the vector ofmodel predictions, which includes degree 2 gravity and shape
coefﬁcients, Love numbers, and libration amplitude at theorbital frequency. Σ is the covariance matrix that contains
contributions from the observational and model covariances:Σ = Σmodel + Σobs. The observation covariance is given by thegravity, shape, tidal Love numbers, and libration determination
accuracies. The observations were assumed independent, thus
yielding a diagonal Σobs. Each Markov chain was inspectedvisually and the initial burn-in steps were discarded. The
Markov chains were run until convergence as informed byvisual inspection and computed chain autocorrelation times.

5. Conclusions
A focused geophysical investigation of Enceladus within a
New Frontiers– or Flagship-class Enceladus Orbiter mission
concept can address the priority science questions outlined inthis paper, as well as other compelling science questions. Thus,
we conclude with the following recommendations:
Recommendation 1: Geophysical measurements should bean essential component of future Enceladus exploration.
Geophysical data can shed light on the mechanism of tidal
dissipation and heat transport. Distinguishing between various
locations of tidal dissipation can be achieved by mapping thevariations of the total ice shell thickness, which can be further
augmented by mapping the elastic shell thickness either by
localizing gravity-topography admittance or by measuring
ﬂexural proﬁles at various locations. Thermal IR mappingaway from the South Polar Terrain is needed to estimate theglobal heat ﬂux. Laser or radar altimetry could be used toderive a high-resolution shape model and measure tidal
deformations.
Recommendation 2: A dedicated gravity mapping investigation should be considered for Enceladus. A GRAIL-likemission would be able to measure tidal phase lag andsigniﬁcantly improve the accuracy of the gravity-topography
admittance. This would enable investigating the energetics ofEnceladus, thus establishing a foundation for understanding
Enceladus’ long-term habitability. That investigation could be
added to New Frontiers– or Flagship-class mission concepts to
Enceladus for little additional resource requirements byleveraging current and upcoming developments in CubeSat
technologies and small radios.

ORCID iDs
Anton I. Ermakov https://orcid.org/0000-0002-7020-7061
Ryan S. Park https://orcid.org/0000-0001-9896-4585
Javier Roa https://orcid.org/0000-0002-0810-1549
Julie C. Castillo-Rogez https://orcid.org/0000-00030400-1038
James T. Keane https://orcid.org/0000-0002-4803-5793
Francis Nimmo https://orcid.org/0000-0003-3573-5915
Edwin S. Kite https://orcid.org/0000-0002-1426-1186
Christophe Sotin https://orcid.org/0000-0003-3947-1072
Gregor Steinbrügge https://orcid.org/0000-0002-1050-7759
Samuel M. Howell https://orcid.org/0000-0002-5126-3228
Bruce G. Bills https://orcid.org/0000-0002-1978-8243
Douglas J. Hemingway https://orcid.org/0000-00015617-207X
Vishnu Viswanathan https://orcid.org/0000-00029027-8588
Gabriel Tobie https://orcid.org/0000-0002-4494-0294
Valery Lainey https://orcid.org/0000-0003-1618-4281
References
A, G., Wahr, J., & Zhong, S. 2014, JGRE, 119, 659
Abramov, O., & Spencer, J. R. 2009, Icar, 199, 189
Akiba, R., Ermakov, A., & Militzer, B. 2021, PSJ, submitted (arXiv:2105.
02790)
Araki, H., Ishibashi, K., Namiki, N., et al. 2019, JSAST, 17, 150
Asmar, S., Armstrong, J., Iess, L., & Tortora, P. 2005, RaSc, 40, RS2001
Baland, R.-M., Yseboodt, M., & Van Hoolst, T. 2016, Icar, 268, 12
Baresi, N., Scheeres, D., & Schaub, H. 2016, AcAau, 123, 364
Beuthe, M., Rivoldini, A., & Trinh, A. 2016, GeoRL, 43, 10,088
Bills, B., & Ermakov, A. 2019, P&SS, 179, 104744
Bland, M. T., McKinnon, W. B., & Schenk, P. M. 2015, Icar, 260, 232
Bland, M. T., Singer, K. N., McKinnon, W. B., & Schenk, P. M. 2012, GeoRL,
39, L17204
Blankenship, D., Ray, T., Plaut, J., et al. 2018, 42nd COSPAR Scientiﬁc
Assembly, 42, B5
Blankenship, D. D., Young, D. A., Moore, W. B., & Moore, J. C. 2009, in
Europa, ed. R. T. Pappalardo et al. (Tucson, AZ: Univ. Arizona Press), 631
Buffett, B. 1996, GeoRL, 23, 2279
Běhounková, M., Souček, O., Hron, J., & Čadek, O. 2017, AsBio, 17, 941
Běhounková, M., Tobie, G., Čadek, O., et al. 2015, NatGe, 8, 601
Běhounková, M., Tobie, G., Choblet, G., & Čadek, O. 2012, Icar, 219, 655
Čadek, O., Souček, O., Běhounková, M., et al. 2019, Icar, 319, 476

A portion of this research was carried out at the Jet
Propulsion Laboratory, California Institute of Technology,
under a contract with the National Aeronautics and Space
Administration (80NM0018D0004).
Appendix
Markov Chain Monte Carlo Internal Structure Inversion
We use the afﬁne-invariant ensemble sampler (Goodman &Weare 2010) implemented in the publicly available emcee
Python library (Foreman-Mackey et al. 2013). First, an
ensemble of internal structure models was generated bysampling the prior probability distribution of the model
parameters. An ensemble consists of individual model realizations or so-called walkers. Typically, from 500 to 1500 walkers
were used. More walkers were used for the MCMC runs withmore stringent measurement errors. The afﬁne-invariant
ensemble sampler uses walkers from the previous step in the13The Planetary Science Journal, 2:157 (14pp), 2021 August

Ermakov et al.

Čadek, O., Tobie, G., Van Hoolst, T., et al. 2016, GeoRL, 43, 5653
Casani, J. R., Strange, N. J., Gibson, M. A., et al. 2020, Enabling a New
Generation of Outer Solar System Missions: Engineering Design Studies for
Nuclear Electric Propulsion, Tech. Rep., 20-2360, NASA JPL
Chen, E. M. A., Nimmo, F., & Glatzmaier, G. A. 2014, Icar, 229, 11
Choblet, G., Tobie, G., Sotin, C., et al. 2017, NatAs, 1, 841
Davis, D. C., Phillips, S. M., & McCarthy, B. P. 2018, AcAau, 143, 16
Dougherty, M. K., Khurana, K. K., Neubauer, F. M., et al. 2006, Sci, 311, 1406
Foreman-Mackey, D. 2016, JOSS, 1, 24
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP,
125, 306
Fuller, J., Luan, J., & Quataert, E. 2016, MNRAS, 458, 3867
Giese, B. 2014, EPSC, 9, 419
Giese, B., Wagner, R., Hussmann, H., et al. 2008, GeoRL, 35, L24204
Glein, C., Postberg, F., & Vance, S. 2018, in Enceladus and the Icy Moons ofSaturn, ed. P. M. Schenk et al. (Tucson, AZ: Univ. Arizona Press), 39
Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics andComputational Science, 5, 65
Gurﬁl, P., Idan, M., & Kasdin, N. 2003, JGCD, 26, 491
Gurﬁl, P., & Kasdin, N. 2004, AcAau, 54, 433
Hansen, C., Esposito, L., Stewart, A., et al. 2006, Sci, 311, 1422
Hay, H. C. F. C., & Matsuyama, I. 2017, Icar, 281, 342
Hedman, M. M., Gosmeyer, C. M., Nicholson, P. D., et al. 2013, Natur,
500, 182
Hemingway, D., Iess, L., Tajeddine, R., & Tobie, G. 2018, in Enceladus andthe Icy Moons of Saturn, ed. P. M. Schenk et al. (Tucson, AZ: Univ.
Arizona Press), 57
Hemingway, D., & Mittal, T. 2019, Icar, 332, 111
Hemingway, D., Rudolph, M., & Manga, M. 2020, NatAs, 4, 234
Howell, K., & Marchand, B. 2005, Dynamical Systems, 20, 149
Howett, C. J. A., Spencer, J. R., Pearl, J., & Segura, M. 2011, JGRE, 116,
E03003Hsu, H.-W., Postberg, F., Sekine, Y., et al. 2015, Natur, 519, 207
Hussmann, H., Lingenauber, K., Kallenbach, R., et al. 2019, CEAS, 11, 381
Iess, L., Stevenson, D. J., Parisi, M., et al. 2014, Sci, 344, 78
Ingersoll, A. P., Ewald, S. P., & Trumbo, S. K. 2020, Icar, 344, 113345
Innanen, K. A. 1979, AJ, 84, 960
Kamata, S., Matsuyama, I., & Nimmo, F. 2015, JGRE, 120, 1528
Kang, W., & Flierl, G. 2020, PNAS, 117, 14764
Kite, E., & Rubin, A. 2016, PNAS, 113, 3972
Kivelson, M., Khurana, K., Russell, C., et al. 2000, Sci, 289, 1340
Konopliv, A., Park, R., Yuan, D.-N., et al. 2013, JGRE, 118, 1415
Kozai, Y. 1962, AJ, 67, 591
Lainey, V., Casajus, L. G., Fuller, J., et al. 2020, NatAs, 4, 1053
Lainey, V., Karatekin, O., Desmars, J., et al. 2012, ApJ, 752, 14
Lara, M., Palacián, J., & Russell, R. 2010, CeMDA, 108, 1
Liao, Y., Nimmo, F., & Neufeld, J. A. 2020, JGRE, 125, e06209
Lidov, M. L. 1962, P&SS, 9, 719
Lippe, C., & D’Amico, S. 2020, AdSpR, 67, 3426
Marusiak, A., Vance, S., Panning, M., et al. 2021, PSJ, in press
Massarweh, L., & Cappuccio, P. 2020, in AIAA Scitech 2020 Forum (Reston,
VA: AIAA), 0467, https://doi.org/10.2514/6.2020-0467
Matsuyama, I., Nimmo, F., Keane, J. T., et al. 2016, GeoRL, 43, 8365
McGovern, P. J., Solomon, S. C., Smith, D. E., et al. 2002, JGRE, 107, 5136
McKay, C., Davila, A., Glein, C., Hand, K., & Stockton, A. 2018, Enceladus
and the Icy Moons of Saturn, ed. P. M. Schenk et al., 437
McKinnon, W. 2013, JGRE, 118, 1775
McKinnon, W. B. 2015, GeoRL, 42, 2137
Meyer, J., & Wisdom, J. 2007, Icar, 188, 535
Meyer, J., & Wisdom, J. 2008, Icar, 193, 213Monteux, J., Collins, G., Tobie, G., & Choblet, G. 2016, Icar, 264, 300
Nadezhdina, I., Zubarev, A., Brusnikin, E., & Oberst, J. 2016, ISPAr, 41, 459
National Research Council 2011, Vision and Voyages for Planetary Science inthe Decade 2013–2022 (Washington, DC: The National Academies Press),
https://doi.org/10.17226/13117
Neveu, M., & Rhoden, A. R. 2019, NatAs, 3, 543
Nimmo, F., Bills, B. G., & Thomas, P. C. 2011, JGRE, 116, E11001
Nimmo, F., Porco, C., & Mitchell, C. 2014, AJ, 148, 46
Ojakangas, G. W., & Stevenson, D. J. 1986, Icar, 66, 341
O’Neill, C., & Nimmo, F. 2010, NatGe, 3, 88
Park, R., Asmar, S., Fahnestock, E., et al. 2012, JSpRo, 49, 390
Park, R. S., Asmar, S. W., Bufﬁngton, B. B., et al. 2011, GeoRL, 38, L22802
Park, R. S., Bills, B., Bufﬁngton, B. B., et al. 2015, P&SS, 112, 10
Park, R. S., Konopliv, A. S., Ermakov, A. I., et al. 2020a, NatAs, 4, 748
Park, R. S., Riedel, J. E., Ermakov, A. I., et al. 2020b, P&SS, 194, 105095
Patterson, G. W., Kattenhorn, S. A., Helfenstein, P., Collins, G. C., &Pappalardo, R. T. 2018, Enceladus and the Icy Moons of Saturn, ed.
P. M. Schenk et al., 95
Peale, S., Cassen, P., & Reynolds, R. 1979, Sci, 203, 892
Porco, C. C., Helfenstein, P., Thomas, P. C., et al. 2006, Sci, 311, 1393
Postberg, F., Khawaja, N., Abel, B., et al. 2018, Natur, 558, 564
Pugh, M., Kuperman, I., Aguirre, F., et al. 2017, in 2017 IEEE Aerospace
Conf. (Piscataway, NJ: IEEE), 1, doi:10.1109/AERO.2017.7943866
Rekier, J., Trinh, A., Triana, S. A., & Dehant, V. 2019, JGRE, 124, 2198
Rignot, E., Mouginot, J., & Scheuchl, B. 2011, Sci, 333, 1427
Roberts, J. 2015, Icar, 258, 54
Roberts, J. H., & Nimmo, F. 2008, Icar, 194, 675
Rovira-Navarro, M., Rieutord, M., Gerkema, T., et al. 2019, Icar, 321, 126
Russell, R. P., & Lara, M. 2009, AcAau, 65, 27
Segall, P. 2010, Volcano and Earthquake Deformation (Cambridge: Cambridge
Univ. Press), doi:10.1515/9781400833856
Soderlund, K. M., Kalousová, K., Buffo, J. J., et al. 2020, SSRv, 216, 80
Souček, O., Běhounková, M., Čadek, O., et al. 2019, Icar, 328, 218
Spencer, J., Pearl, J., Segura, M., et al. 2006, Sci, 311, 1401
Spencer, J. R., Howett, C. J. A., Verbiscer, A., et al. 2013, EPSC, 8, 840
Spencer, J. R., Nimmo, F., Ingersoll, A. P., et al. 2018, in Enceladus and theIcy Moons of Saturn, ed. P. M. Schenk et al. (Tucson, AZ: Univ. Arizona
Press), 163
Spitale, J. N., Hurford, T. A., Rhoden, A. R., Berkson, E. E., & Platts, S. S.
2015, Natur, 521, 57
Steinbrügge, G., Schroeder, D., Haynes, M., et al. 2018, E&PSL, 482, 334
Steinbrügge, G., Stark, A., Hussmann, H., Sohl, F., & Oberst, J. 2015, P&SS,
117, 184
Tajeddine, R., Soderlund, K. M., Thomas, P. C., et al. 2017, Icar, 295, 46
Tapley, B., Bettadpur, S., Watkins, M., & Reigber, C. 2004, GeoRL, 31, L09607
Thomas, P. C., Tajeddine, R., Tiscareno, M. S., et al. 2016, Icar, 264, 37
Tricarico, P. 2014, ApJ, 782, 99
Tyler, R. 2020, Icar, 348, 113821
Van Hoolst, T., Baland, R., & Trinh, A. 2016, Icar, 277, 311
Van Hoolst, T., Baland, R.-M., & Trinh, A. 2013, Icar, 226, 299
Van Hoolst, T., Rambaux, N., Karatekin, O., Dehant, V., & Rivoldini, A. 2008,
Icar, 195, 386
Wahr, J. M., Zuber, M. T., Smith, D. E., & Lunine, J. I. 2006, JGRE, 111,
E12005
Waite, J. H., Glein, C., Perryman, R., et al. 2017, Sci, 356, 155
Watts, A. B. 2001, Isostasy and Flexure of the Lithosphere (Cambridge:
Cambridge Univ. Press)
Wu, Y. 2005, ApJ, 635, 688
Yin, A., & Pappalardo, R. T. 2015, Icar, 260, 409
Zuber, M., Smith, D., Watkins, M., et al. 2013, Sci, 339, 668

14