<figure><img src="image_1.png" /></figure>

HAL
open science

# A Recipe for the Geophysical Exploration of Enceladus

Anton Ermakov, Ryan Park, Javier Roa, Julie Castillo-Rogez, James Keane, Francis Nimmo, Edwin Kite, Christophe Sotin, T. Joseph W. Lazio, Gregor Steinbrügge, et al.

► To cite this version:

Anton Ermakov, Ryan Park, Javier Roa, Julie Castillo-Rogez, James Keane, et al. A Recipe for the Geophysical Exploration of Enceladus. The Planetary Science Journal, 2021, 2 (4), pp.157. 10.3847/PSJ/ac06d2. hal-03468856

HAL Id: hal-03468856

https://hal.science/hal-03468856v1

Submitted on 30 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Recipe for the Geophysical Exploration of Enceladus

Anton I. Ermakov<fnref n="1" />, Ryan S. Park<fnref n="2" />, Javier Roa<fnref n="2" />, Julie C. Castillo-Rogez<fnref n="2" />, James T. Keane<fnref n="2" />, Francis Nimmo<fnref n="3" />,
Edwin S. Kite<fnref n="4" />, Christophe Sotin<fnref n="5" />, T. Joseph W. Lazio<fnref n="2" />, Gregor Steinbrügge<fnref n="6" />, Samuel M. Howell<fnref n="2" />, Bruce G. Bills<fnref n="2" />,
Douglas J. Hemingway<fnref n="7" />,<fnref n="8" />, Vishnu Viswanathan<fnref n="9" />,<fnref n="10" />, Gabriel Tobie<fnref n="11" />, and Valery Lainey<fnref n="12" />

## Abstract

Orbital geophysical investigations of Enceladus are critical to understanding its energy budget. In this paper, we identified 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 of single-spacecraft and dual-spacecraft configurations to static gravity and tidal Love numbers of Enceladus. We find that mapping Enceladus' gravity field, 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 for assessing Enceladus' long-term habitability. This kind of investigation could be carried out as part of a life search mission at little additional resource requirements.

*Unified Astronomy Thesaurus concepts: Planetary interior (1248); Gravitational fields (667); Saturnian satellites (1427); Space observatories (1543); Astrodynamics (76)*

## 1. Introduction

Enceladus—a cryovolcanically active and apparently habitable satellite in the Saturn system—challenges our understanding of geodynamical processes governing the evolution of ocean worlds (Waite et al. 2017; McKay et al. 2018). Starting in 2005, the Cassini spacecraft revealed active eruptions in Enceladus' southern hemisphere (Porco et al. 2006), confirmed 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 flexing. This energy is necessary for Enceladus' ocean to persist over geologic timescales (Roberts & Nimmo 2008).

The goals of this paper are:

1. To provide an overview of the science questions that orbit-based geophysical data at Enceladus could answer;

2. To offer recommendations for optimizing the future collection 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 or a large Flagship-class mission (National Research Council 2011).

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; McKinnon 2015) 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 flux 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 a longer-term variability (Ingersoll et al. 2020). The current state of Enceladus' geophysical data is summarized in Table 1 and illustrated in Figure 1.

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ University of California, Berkeley, McCone Hall 307, CA 94720, USA; eai@berkeley.edu</footnote>

<footnote n="2">² Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="3">³ Department of Earth and Planetary Sciences, University of California, Santa Cruz, Santa Cruz, CA 95064, USA</footnote>

<footnote n="4">⁴ Department of Geophysical Sciences, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="5">⁵ Laboratoire de Planétologie et Géodynamique, Université de Nantes, F-44322 Nantes, France</footnote>

<footnote n="6">⁶ Department of Geophysics, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="7">⁷ Earth & Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA</footnote>

<footnote n="8">⁸ Maxar Technologies, 3875 Fabian Way, Palo Alto, CA 94303, USA</footnote>

<footnote n="9">⁹ NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA</footnote>

<footnote n="10">¹⁰ University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250, USA</footnote>

<footnote n="11">¹¹ Laboratoire de Planétologie et Géodynamique, CNRS/Université de Nantes, France</footnote>

<footnote n="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</footnote>
<table><caption>**Table 1**
Summary of the Orbit-based Enceladus Geophysical Data</caption><thead><tr><td>Quantity</td><td>Current Knowledge</td><td>Information Provided</td></tr></thead><tbody><tr><td>Gravity field</td><td>Up to degree 3 field measured by Cassini (Iess et al. 2014).</td><td>Constraint on the internal structure and compensation mechanism in the ice shell.</td></tr><tr><td>Shape</td><td>Estimated up to degree 16 (Tajeddine et al. 2017).</td><td>Constraint on the internal structure and ice shell thickness variations.</td></tr><tr><td>Gravity-topography admittance<sup>a</sup></td><td>Degree 3 admittance derived (Iess et al. 2014).</td><td>Constraint on compensation state of topography, total shell thickness, elastic shell thickness, and shell density and shell-ocean density contrast.</td></tr><tr><td>Obliquity</td><td>Cassini-derived upper bound of 0°05 (Giese 2014). Upper theoretical bound 4 × 10<sup>-4</sup>° (Baland et al. 2016).</td><td>Enables determination of moment of inertia and whether Enceladus is tidally damped into a Cassini state.</td></tr><tr><td>Physical libration amplitude</td><td>Measured amplitude at the orbital frequency of 0°120 ± 0°014 (Thomas et al. 2016) or 0°155 ± 0°014 (Nadezhdina et al. 2016).</td><td>Libration amplitude is strongly sensitive to the shell thickness. Observed large amplitude requires a decoupled ice shell and hence implies a global subsurface ocean.</td></tr><tr><td>Precession and nutation</td><td>Not currently measured. Precession rate estimated 2.6 rad yr<sup>-1</sup> (Baland et al. 2016).</td><td>In combination with degree 2 gravity, enables moment of inertia determination without the need for hydrostatic equilibrium assumption.</td></tr><tr><td>Tidally driven orbital migration</td><td>Measured from historical astrometric data (Lainey et al. 2012).</td><td>Provides a constraint on the tidal dissipation factor Q of Saturn.</td></tr><tr><td>Potential Love number k<sub>2</sub>; radial and lateral displacement Love numbers h<sub>2</sub> and l<sub>2</sub>, respectively.</td><td>Not currently measured.</td><td>Real parts are mostly sensitive to the thickness of the icy shell and its rigidity. Imaginary parts are sensitive to the viscosity profile and, thus, provide a constraint on total tidal dissipation within Enceladus.</td></tr><tr><td>Radar subsurface mapping</td><td>Not currently achieved.</td><td>Enables independent shell thickness determination.</td></tr><tr><td>Magnetic induction</td><td>Not currently measured. The magnitude of the forcing field is ≈10 nT.</td><td>Sensitive to the subsurface ocean conductivity and thickness.</td></tr><tr><td>Thermal IR mapping</td><td>15.8 ± 3.1 GW from SPT (Howett et al. 2011); 4.2 GW from Tiger Stripes (Spencer et al. 2013).</td><td>Constraint on the total heat flux.</td></tr></tbody></table>

**Note.**

Tidal dissipation within Enceladus' interior is a function of its orbital characteristics (e.g., eccentricity and proximity to Saturn) 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 present-day internal structure determines the instantaneous spatially variable tidal dissipation rate. Depending on the efficiency 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, the heat 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 from astrometric 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 the current internal structure and dissipation within Enceladus, we can place constraints on its tidal migration history, which is interlocked with the history of tidal heating.

We have identified the following interrelated priority science questions that should be addressed by future geophysical observations 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 estimated radius of 191–198 km, overlain by a global ocean with an estimated thickness of 30–39 km and an icy shell with an estimated average thickness of 19–24 km, but thinner at the south pole (4–12 km; Hemingway & Mittal 2019).

<footnote><sup>a</sup> Gravity-topography admittance is not an independent geophysical observable. It is derived as the ratio of the gravity amplitude to topography amplitude.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic illustration of Enceladus' internal structure along with geophysical measurements and science priority questions identified in this paper.</figcaption></figure>

However, these inferences rely on plausible, yet untested, assumptions. In particular, the core is assumed to be in hydrostatic equilibrium, which was challenged by McKinnon (2013) and Monteux et al. (2016). The gravity-topography admittance, defined as the wavelength-dependent ratio of gravity amplitude to topography amplitude, was assumed isotropic, which would be violated if there are significant deviations 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 are relaxed, it becomes impossible to separate geophysical signals that originate in the icy shell (e.g., freezing/melting of the shell, convection, impact cratering, tectonics) from those originating in the core (e.g., core's nonhydrostaticity and core tidal dissipation; Roberts 2015). Below, we summarize how a future 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 the libration 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 (Park et 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 to smaller 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 deep

topographic depression (Nimmo et al. 2011; Tajeddine et al. 2017) and hosts a set of prominent fractures, called Tiger Stripes, that emit plumes of H₂O, salt, organics, and other volatiles 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 are complemented by radar sounding, an independent constraint on spatial variations in shell thickness could be derived.

Mapping gravity and topography would also be useful in more heavily cratered terrains in Enceladus' northern hemisphere. The relaxed state of craters in that region indicates higher heat fluxes in the past as such heating episodes would facilitate viscous flow 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 be assessed 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, that would provide an independent constraint on the moment of inertia of the shell. However, the very small amplitude of the deflection (order of 1 m) would be challenging to measure, owing to Enceladus' small obliquity (see Table 1).

**Partitioning nonhydrostaticity between the shell and the core**—Separating the nonhydrostatic signal of the shell from that of 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 and convection zones to be identified, allowing the signals originating within the ice shell to be separated. A direct detection of the ice-
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 to better constrain the shell structure, specifically its moment of inertia differences (such as the difference between the shell's equatorial 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 and density are indirectly constrained by gravity and shape data by satisfying the mass balance. The current uncertainty on the ocean 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 and the shell. Thus, if the shell density is constrained from the high-degree 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 field is axisymmetric. Thus, the changing magnetic field that Enceladus experiences is primarily due to its eccentric orbit, resulting in a much smaller forcing field amplitude (≈10 nT). This signal is of roughly the same magnitude as time-variable fields arising from the plumes (Dougherty et al. 2006). Thus, detecting an induction response and disentangling it from the plume-induced variability in the magnetic field would require long-term observations and characterization of the plume activity.

**Rheology of the material**—The tidal response of Enceladus is described by its Love numbers $k_{nm}$, $h_{nm}$, and $l_{nm}$ (e.g., Wahr et al. 2006), characterizing potential change, radial, and lateral displacements, respectively, where *n* is spherical harmonic degree and *m* is order. For a spherically symmetric body, the Love numbers are degenerate with respect to *m*. Thus, for such a body: $k_{20} = k_{21} = k_{22}$. Thus, we can drop the second subscript for shorter notation if spherical symmetry is assumed. Separating the tidal signal corresponding to different orders within the same degree would indicate deviations of the body's mechanical properties from spherical symmetry. The Love numbers are complex quantities and quantify how Enceladus' gravity field and shape respond to time-varying tidal forces. The real parts of the Love numbers depend primarily on its shell thickness and shell rigidity. The imaginary parts depend on the rheology of the material and quantify the lag in the tidal response with respect to tidal forcing, which is related to the current total dissipation within Enceladus. Independent measurement of both $h_2$ (the degree 2 radial displacement Love number, omitting the order index for simplicity) and $k_2$ (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 from tidal 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 spherically symmetric shell when modeling the expected tidal response. A

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

## 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 flux of Enceladus has been measured by mapping out the surface thermal emission using the Composite Infrared Spectrometer (CIRS) on board Cassini. The estimated value of the heat flux has varied over the years of Cassini data analysis by a factor of several. All studies to date have focused only on the heat flux localized in the South Polar Terrain. Spencer et al. (2006) estimated a heat flux of $5.8 \pm 1.9$ GW. A subsequent study by Howett et al. (2011) yielded a larger estimate of $15.8 \pm 3.1$ GW using observations from a higher-wavelength CIRS detector. Spencer et al. (2013) used higher spatial resolution data and derived 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 flux from flyby observations of heterogeneous quality.

The conductive and advective heat fluxes are approximately equal in the South Polar Terrain (Kite & Rubin 2016). Advection likely dominates within the Tiger Stripes, causing prominent temperature anomalies. The heat flux away from the South Polar Terrain remains unconstrained. Current estimates on global heat flux 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 flow.

The local instantaneous tidal heat generation is determined by the product of tidal stress and strain rate. The mutual relation between these two quantities is, in turn, set by the rheology of the material that is highly sensitive to temperature. Heat is likely generated in the warmer ice at the base of the shell 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 the South Polar Terrain faults and concluded that dissipation in the solid ice cannot exceed 2.1 GW, implying that an additional heat source is needed to explain the observed heat flow. Additional dissipation could occur within the liquid slots in the ice (Kite & Rubin 2016). Over 10 GW can be generated by tidal friction inside the unconsolidated rocky core (Choblet et al. 2017) and comparable amounts by tidal flushing of water through the porous core (Liao et al. 2020). Theoretical models predict that the ocean heat production is likely negligible (Chen et 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 be enhanced in narrow liquid-water conduits within the ice shell (Kite & Rubin 2016). However, more recently, Tyler (2020) indicated that there could be resonant configurations, where ocean tidal dissipation could create significant heating.

If Enceladus is in thermal equilibrium, there exists a relationship between the shell thickness and the conductive heat flux. Thus, the scientific requirement driving the uncertainty in the mean shell thickness measurement can be tied to the desired uncertainty in the heat flux. It follows from Hemingway & Mittal (2019) that in order to determine the conductive heat flux to within 10%, the average shell thickness needs to be determined
to within 2 km. In the South Polar Terrain, where the ice shell could be <5 km thick, radar observations may more easily reach the ice-ocean interface and return the shell thickness with an accuracy 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 by their cross-analysis.

Distinguishing between convection and conduction is of great interest. Conduction is likely dominant in the outer part of the 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. Efficient heat transport by convection might lead to quick ocean freezing. Large-amplitude shell thickness variations inferred from 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 the icy 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 profile. For a convective shell, warmer ice at shallower depth could lead to higher radar signal attenuation, especially if the salinity of the ice is high. A loss of signal resulting in no detected radar reflections at depth or weaker radar reflections 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 gravity-topography admittance) and the thickness derived from the Love numbers, because the latter is sensitive only to the elastic part of the shell and not the weak convecting part. If the convecting part of 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 reflection in the radar return. However, this direct signal could be further covered up by an accretion zone at the ice-ocean interface, which could produce ice of marine composition or result in a mushy layer with high attenuation.

Estimating the elastic thickness would provide a critical constraint 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 flexural profiles of tectonic features that would require accurate topography knowledge (e.g., Giese et al. 2008). A flexural profile 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 flow mapped in the thermal IR emission could help validate the dissipation pattern within the shell. Measuring flexural 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 be

further improved and geodetically referenced to the center-of-mass 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, can change over the tidal migration history. In addition, *Q* may evolve as the temperature of the body evolves due to secular cooling. Recent astrometric efforts indicate that Saturn is more dissipative than previously thought. Lainey et al. (2012) found Saturn'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 to two distinct kinds of steady state. First, in the thermal steady state, 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 balanced by 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 is stored 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 for viscous 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 flow and associated tectonic activity has changed over geologic time. Geophysical investigation of these tectonic anomalies may provide clues as to non-steady-state behavior in Enceladus' past: for example, the processes that drove the onset of hyperactive resurfacing near the South Pole, but not the North Pole (e.g., Kang & Flierl 2020).

Both kinds of steady states critically depend on the tidal phase lag, which is characterized by the imaginary part of the potential Love number, $Im(k_2)$. Larger values of $Im(k_2)$ correspond to larger 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 the perturbing potential. Thus, it can be derived by radio-tracking of the spacecraft in the same way the gravity field 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- or spacecraft-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 of both Enceladus and Dione.

In conclusion, long-baseline continuous ground-based astrometry at the kilometer level and future radio ranging and high-accuracy 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,
including the conductive heat flux away from the Tiger Stripes, is needed to assess if Enceladus is in the thermal steady state.

## 2.4. What Are the Feedbacks between Volcanism and Tectonics That Regulate Enceladus’ Cryovolcanism?

Enceladus cryovolcanism leads to significant mass loss from the moon (150–350 kg s⁻¹ or 4%–10% of Enceladus’ mass per Gyr; Hansen et al. 2006) and also offers the opportunity to sample material from Enceladus’ ocean. However, the causes of 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 well as rates of resurfacing and surface-interior exchange, which are important 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 present-day rate of cryovolcanic activity is representative of the history of the South Polar Terrain (O’Neill & Nimmo 2010).

Enceladus operates in a regime intriguingly different from the other worlds known to have active volcanism such as Earth and Io. Earth’s oceanic lithosphere is dominantly cooled by conduction (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 flow. Estimates of the elastic thickness from the flexural profiles collected in the vicinity of the cracks might be significantly different from the estimates using localized gravity-topography admittance as a significant advection 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 the plumbing 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 reflectors for electromagnetic waves and manifest themselves in the radar return. Vertical fractures could be identified by point scattering, usually leaving characteristic hyperbolas in the (nonfocused) radar return. Furthermore, compositional variations associated with tectonic activity, injection of water, or even variations in crystalline fabric can be detected in radar sounding data. Due to the complex surface topography of Enceladus, the radar sounding data would need to be complemented by high-resolution topography data to mitigate surface clutter.

Higher-resolution gravity, shape, and heat flux, as well as radar sounding and data within the South Polar Terrain, would open the prospect of studying a new kind of tectonics as rich and unusual by terrestrial standards as plate tectonics when first 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

<table><caption>**Table 2**

Prior Probability Distribution of Enceladus’ Internal Structure Parameters</caption><thead><tr><th>Parameter</th><th>Minimum Value</th><th>Maximum Value</th><th>Distribution Type</th></tr></thead><tbody><tr><td>Shell shear modulus [GPa]</td><td>1</td><td>10</td><td>log-uniform</td></tr><tr><td>Shell Poisson's ratio</td><td>0.25</td><td>0.45</td><td>uniform</td></tr><tr><td>Shell viscosity [Pa s]</td><td>10<sup>14</sup></td><td>10<sup>19</sup></td><td>log-uniform</td></tr><tr><td>Shell density [kg m<sup>-3</sup>]</td><td>600</td><td>1100</td><td>uniform</td></tr><tr><td>Shell thickness [km]</td><td>2</td><td>40</td><td>uniform</td></tr><tr><td>Ocean compressibility [GPa]</td><td>2</td><td>3</td><td>log-uniform</td></tr><tr><td>Ocean density [kg m<sup>-3</sup>]</td><td>920</td><td>1300</td><td>uniform</td></tr><tr><td>Ocean thickness [km]</td><td>2</td><td>80</td><td>uniform</td></tr><tr><td>Core shear modulus [GPa]</td><td>10</td><td>100</td><td>log-uniform</td></tr><tr><td>Core Poisson's ratio</td><td>0.25</td><td>0.45</td><td>uniform</td></tr><tr><td>Core viscosity [Pa s]</td><td>10<sup>20</sup></td><td>10<sup>22</sup></td><td>log-uniform</td></tr><tr><td>Core density [kg m<sup>-3</sup>]</td><td>2000</td><td>3000</td><td>uniform</td></tr><tr><td>Core radius [km]</td><td>2</td><td>252</td><td>uniform</td></tr></tbody></table>

**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 flow 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 can yield 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 the recovered parameter. For example, the libration amplitude is an inverse function of shell thickness.

In order to develop traceable measurement requirements, we used the Markov Chain Monte Carlo (MCMC) approach to develop 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 prior probability distribution for all model parameters, including the currently unconstrained (due to a lack of Love number measurements) viscoelastic moduli. The form of the prior probability distribution is summarized in Table 2.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Corner plot visualizing the posterior distribution of internal structure model parameters. Histograms for layer thicknesses ($h_i$) and shear modulus of the icy shell ($\mu_{shell}$) are shown. Note that $\mu_{shell}$ is shown on a logarithmic scale. $\mu_{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. $Z_l$ stands for the gravity topography admittance at degree $l$. The table at the top shows the measurement uncertainties for the corresponding contours. The infinity sign in the table indicates that the quantity was not included in the inversion.</figcaption></figure>

To generate synthetic observations, we compute a multilayer hydrostatic equilibrium model using Tricarico (2014). Viscoelastic moduli and densities are assumed constant within each layer. This gives us the hydrostatic shape and gravity spherical harmonic coefficients. Degree 2 complex Love numbers are computed in a way similar to Kamata et al. (2015). Enceladus is assumed to have a small nonhydrostatic topography, which is compensated with an Airy isostasy mechanism (Watts 2001). This yields modeled gravity-topography admittance. The correlation between gravity and topography is assumed to be unity. The libration amplitude at the orbital frequency is computed for rigid shells (Van Hoolst et al. 2008), which is a simplification of our model. As the libration measurement accuracy is improved, it becomes sensitive to the rigidity of the shell. Čadek et al. (2016) estimated that varying the shell rigidity between 1 and 5 GPa leads to a libration amplitude difference of 50 m, which is just below current observational uncertainty. The variation of libration amplitude due to the shell viscosity profile is yet smaller: <20 m for bottom shell viscosity ranging from $10^{12}$ to $10^{15}$ 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 coefficients will not be worse than the current ones (Iess et al. 2014; Thomas et al. 2016). We explore the combined sensitivity of physical libration amplitude at the orbital period; tidal Love numbers $k_2$, $h_2$, and $l_2$; and the gravity-topography admittance spectrum $Z_n$, to the internal structure parameters. For each set of measurement errors, we derive the posterior 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 the posterior distributions for selected internal structure parameters in the form of a corner plot (Foreman-Mackey 2016). The top boxes show the 1D histograms of the parameters of interest. In
<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

addition, we show 95% confidence regions for selected parameter pairs. The black contours in Figure 2 represent the current state of knowledge. The other (colored) contours show how the confidence regions shrink as more data are included in the MCMC inversion.

We find 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⁻¹ up to degree 10 and libration 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 flux. The Love numbers provide sensitivity to viscoelastic moduli. If $k_2$ is measured to $10^{-2}$ or better, it provides a constraint on shell rigidity and improves the accuracy of the shell thickness determination. If $k_2$ is measured to $10^{-3}$, it becomes sensitive to the shell viscosity.

### 3.1. Estimating the Amplitude of Tidal Deformations

Having estimated Enceladus’ Love numbers using an MCMC inversion, we can statistically assess the magnitude of surface tidal deformation. The posterior distribution for the Love numbers is shown in Figure 3. The derived confidence intervals for the Love numbers allow an estimation of the tidal displacement ranges and surface gravity changes to be measured by future spacecraft missions. The tidal surface displacement is proportional to the displacement Love numbers $h_2$ and $l_2$. The maximum range of tidal displacement from degree 2 tides is given by Park et al. (2020b) as

$$ \max(\Delta R) = 12 \sqrt{\frac{2}{7} \frac{h_2 R^2 \omega^2 e}{g}}, \quad (1) $$

$$ \max(\Delta N) = 9 \frac{l_2 R^2 \omega^2 e}{g}, \quad (2) $$
<figure><img src="image_7.png" /><figcaption>**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 as constrained by the gravity (Iess et al. 2014), shape (Tajeddine et al. 2017), and libration data (Thomas et al. 2016).</figcaption></figure>

<figure><img src="image_8.png" /></figure>

$$ \max(\Delta E) = 12 \frac{l_2 R^2 \omega^2 e}{g}, \qquad (3) $$

for the radial, northerly, and easterly displacements, respectively. Here, R is the satellite's radius, e is its orbital eccentricity, ω is the orbital frequency, and g is the surface gravity. In addition, the maximum surface gravity variation range is given by

$$ \max(\Delta g) = 24 \sqrt{\frac{2}{7}} \left(1 - \frac{3}{2}k_2 + h_2\right) R\omega^2 e. \qquad (4) $$

Taking the median values for the Love numbers from Figure 3 ($k_2$ = 0.0167, $h_2$ = 0.0429, $l_2$ = 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 symmetric model, radial tidal displacements are maximized at the equator.

The maximum radial tidal deformation range is 0.5–6.5 m (2.5th–95th confidence interval) in the equatorial region but can be amplified in the vicinity of the Tiger Stripes (Běhounková et al. 2017; Marusiak et al. 2021). We note that tidal displacements are expected to be much smaller compared to 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><img src="image_9.png" /><figcaption>**Figure 5.** Tidal dissipation of Enceladus as a function of the imaginary part of the Love number.</figcaption></figure>

Given the observed heat flux (Howett et al. 2011), we can derive an estimate of the imaginary part of the potential Love number using the tidal dissipation formula (Peale et al. 1979):

$$ \dot{E} = \operatorname{Im}(k_2) \frac{21 GM_{\mathrm{Saturn}}^2 R_{\mathrm{Enceladus}}^5 \omega e^2}{2 a^6}, \qquad (5) $$

where G is the universal gravitational constant, and a is the semimajor axis of Enceladus' orbit. Thus, the measurement requirement on $\operatorname{Im}(k_2)$ should be tied to the requirement on the heat flux measurement. Figure 5 shows the total tidal
dissipation as a function of Im($k_2$), from which follows that measuring 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 $h_2$ and $l_2$ with the same accuracy as for $k_2$ can help mitigate the ambiguity of the ice rheology that arises when measuring $k_2$ only (Wahr et al. 2006).

# 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 with interspacecraft tracking (i.e., GRACE/GRAIL-like configuration, Tapley et al. 2004; Zuber et al. 2013) to map the static and temporally variable gravity field 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 high-inclination orbits need to be considered (Russell & Lara 2009; Massarweh & Cappuccio 2020). High-inclination orbits are also of interest for gravity mapping of the South Polar Terrain of Enceladus that lies below 60° S.

Tidal deformation causes the gravity field of Enceladus to vary in time, which affects the motion of spacecraft around it. When operating a formation of multiple spacecraft around Enceladus, the tidal investigation would rely both on the effect of tides on each spacecraft but also on the differential effect of the perturbation, which affects the spacecraft relative motion. The effect of tides on the motion of a spacecraft can be split into two components:

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 variations in the orbital elements or long-term changes in the range rate between two spacecraft. The indirect effect may be amplified by resonances between the spacecraft's mean-motion 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 optically observing tidal deformation and detecting tidally induced gravity variations as they would maximize the observable signal. 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 to order m, with degree 2 Love numbers differing by up to a factor of 2 (Běhounková et al. 2017). In order to separate the nondegenerate Love numbers, a nonequatorial orbit would be needed. Thus, the choice of an optimal orbit should depend on the expected internal structure and should be driven by a hypothesis to maximize the sensitivity to the phenomenon in question.

## 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

spacecraft remaining in orbit without escaping or impacting Enceladus. In practice, a certain orbit may be considered stable as long as its lifetime, although finite, is long enough to complete the scientific exploration phase according to the mission requirements.

Tidal perturbations due to Enceladus' tidal deformation are small compared to the perturbation from Saturn—even compared to the third-order terms of the static gravity field—and do not pose a risk to orbital stability. For this reason, it is relevant to find the orbit configurations 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 scientific 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 first-order expansion of the three-body potential leads to the quadrupole approximation, which predicts the maximum eccentricity during an LK cycle to be

$$e_{\max} = \sqrt{1 - \frac{5}{3} \cos^2 i}, \quad (6)$$

where i is the orbital inclination. This approximation is only valid for $\cos^2 i < 3/5$, or for inclinations between the Kozai angles $i_{min} = 39.2^\circ$ and $i_{max} = 140.8^\circ$. The eccentricity of initially circular orbits with inclinations between these critical values 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 to estimate the inclination at which the periapsis reaches the surface of Enceladus. However, Lara et al. (2010) proved that the quadrupole expansion is not sufficient to accurately model the dynamics at Enceladus and found that a higher-order expansion predicts stable near-circular orbits at inclinations as high as 50°.

The maximum altitude that low-inclination near-circular orbits 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 effect opposes 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 stable below 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 also reach 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 specific objectives. For example, the low-altitude near-circular equatorial orbits already discussed are members of the family of distant prograde orbits (DPOs) and distant retrograde orbits (DROs),
respectively. 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 and observing the Tiger Stripes, a configuration 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 the simplified 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 to observe the South Polar Terrain. Although such highly inclined orbits tend to be unstable in the full-ephemeris model, Davis et al. (2018) demonstrated that halo orbits flying over the South Polar Terrain at altitudes lower than 200 km can be controlled with an estimated $\Delta V$ cost of $20 \text{ m s}^{-1}$ per month. Nuclear electric propulsion can significantly 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 of formation flying concepts for deep-space exploration has been discussed in the past, including design and control considerations (Gurfil et al. 2003; Howell & Marchand 2005). Examples of proposed applications include close-proximity exploration of small bodies (Baresi et al. 2016; Lippe & D'Amico 2020) and formations along generic periodic orbits in the CR3BP (Gurfil & Kasdin 2004). Still, formation flying in close proximity to tidally locked satellites remains an understudied area.

For the case of Enceladus, the perturbation from Saturn's and Enceladus' nonuniform gravity can produce strong differential accelerations on the spacecraft and destabilize the formation. Passively mitigating the differential effect of the perturbations during the design stage is crucial for formulating an efficient control plan that minimizes the overall cost of station-keeping. Active control may still be required to maintain the formation within operational requirements. Precise autonomous guidance, navigation, and control are required to ensure that the spacecraft can react to and counteract perturbations efficiently. Constant maneuvering can interfere with the continuity of the radio-tracking data. Thus, station-keeping maneuvers should be kept at a minimum when possible.

### 4.3. Mission Simulations

Based on the aforementioned considerations, we conducted a series of mission simulations and performed detailed covariance analyses for two mission configurations: a single orbiter with radiometric tracking to Earth and dual spacecraft (GRAIL-like configuration) 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 for spacecraft radio-tracking Doppler measurements. These bands respectively yield ranging accuracies of $\approx 5 \times 10^{-7}$, $5 \times 10^{-8}$,

and $5 \times 10^{-9} \text{ km s}^{-1}$ 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 the two-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 distance encountering 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 the latter 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 single-spacecraft case.

We assumed an X-band tracking accuracy for the single-spacecraft case, which is typical for a deep-space mission. We only considered X-band for single spacecraft because the Ka-band uplink capability is currently available only at the Goldstone DSN station. In addition, the Ka-band uplink capability is costly to maintain as it requires a water vapor radiometer to get the full accuracy. If a Ka-band tracking of a single 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 and 200 km, respectively. Both poles can be covered with flybys prior to being placed into this stable orbit, depending on the science requirements, which was not modeled in our simulations. We have simulated how accurately we can recover Enceladus' geophysical parameters, such as Love numbers, gravity field, 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 features of interest.

Assuming a single orbiter with two-way X-band tracking capability and assuming the Doppler data accuracy of $10^{-7} \text{ km s}^{-1}$ at a 60 s count time with a 28 day data collection duration, a gravity field up to degree $\approx 10$ can be recovered and $k_2$ can be recovered with an accuracy of $\approx 10^{-2}$. A high-resolution context imager with accurate pointing knowledge and image-motion-compensation capability, such as the Advanced Pointing Imaging Camera (APIC; Park et al. 2020b), can measure $h_2$ and $l_2$ with accuracies $\approx 0.01$ and $\approx 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 \text{ m}$ at the equator (a relative uncertainty of $<0.2\%$). The accuracy of $h_2$, $l_2$, and libration amplitude can be improved with the optimal distribution of crossover points (Park et al. 2015) or lower orbit altitude.

Recovering $h_2$ 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-of-the-art laser altimeters can achieve vertical resolutions of 10 cm (Araki et al. 2019). However, the error budget for the
<figure><img src="image_10.png" /><figcaption>**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, $\lambda \approx 2\pi R/n$. The light blue region indicates the expected range of the gravity signal. The upper bound is given by a power law fitted to an uncompensated gravity-from-shape using the shape model by Tajeddine et al. (2017) with a shell density of 920 kg m⁻³. The lower bound corresponds to gravity from compensated topography using the same shape model with a shell density of 600 kg m⁻³. The large-scale depressions refer to the chain of basins identified by Tajeddine et al. (2017).</figcaption></figure>

<table><caption>**Table 3**
A Summary of the Potential Love Numbers and Gravity Recovery</caption><thead><tr><th></th><th>Range-rate Accuracy (km s<sup>-1</sup>)</th><th>k<sub>20</sub></th><th>k<sub>21</sub></th><th>k<sub>22</sub></th><th>Global Degree Strength</th></tr></thead><tbody><tr><td rowspan="2">Single spacecraft</td><td>10<sup>-7</sup></td><td>9 &times; 10<sup>-3</sup></td><td>5 &times; 10<sup>-1</sup></td><td>4 &times; 10<sup>-3</sup></td><td>9</td></tr><tr><td>10<sup>-8</sup></td><td>9 &times; 10<sup>-4</sup></td><td>2 &times; 10<sup>-2</sup></td><td>2 &times; 10<sup>-4</sup></td><td>17</td></tr><tr><td rowspan="2">Dual spacecraft</td><td>10<sup>-9</sup></td><td>2 &times; 10<sup>-4</sup></td><td>3 &times; 10<sup>-3</sup></td><td>2 &times; 10<sup>-5</sup></td><td>23</td></tr><tr><td>10<sup>-10</sup></td><td>2 &times; 10<sup>-5</sup></td><td>4 &times; 10<sup>-4</sup></td><td>3 &times; 10<sup>-6</sup></td><td>30</td></tr></tbody></table>

**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 defines the global resolution of the gravity model.

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 beneficial for the $h_2$ measurement as well. Typical accuracies for such a configuration would be on the order of 1 m. The accuracy of the $h_2$ inversion then depends on the number of orbits, $n_{orb}$, a spacecraft completes within its mission lifetime. Because the number of crossover points increases with $n_{orb}^2$, millions of crossover points can be collected over mission lifetimes of a few months.

Using a radar sounder instead of a dedicated altimeter requires 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 the center frequency. Ice is particularly transparent to frequencies between 1 and 300 MHz (Blankenship et al. 2009). To achieve a 3 m vertical resolution, a 50 MHz bandwidth is required, which is larger than the bandwidth used by typical radar sounders. However, by using stereo imaging in combination with radar sounders, resolutions up to a factor of 4–5 better

than the inherent range resolution can be achieved (Steinbrügge et al. 2018).

A dual-spacecraft architecture with intersatellite Ka-band tracking (i.e., GRAIL-like scenario) was also considered for the same orbit configuration described above. We assumed range-rate accuracies of $10^{-8}$, $10^{-9}$, and $10^{-10}$ km s⁻¹. These or better (down to $3 \times 10^{-11}$ km s⁻¹) ranging accuracies were achieved by the GRAIL mission (Konopliv et al. 2013). Such measurement configuration can significantly improve the accuracy of the $k_2$ determination yielding a $k_2$ error down to $\approx 2 \times 10^{-5}$ (or relative accuracy of $\approx 0.1\%$; Table 3). The gravity field can be determined to degrees 20–30 (see Figure 6), depending on the ranging accuracy. This level of accuracy would allow a definitive determination of the tidal phase lag and total tidal dissipation within Enceladus. The core viscoelastic moduli and ocean compressibility would remain virtually unconstrained for the studied mission configurations, although a lower bound on the core viscosity could be retrieved from the Love numbers if they are measured to an accuracy of 0.1%.
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⁻¹⁴ at 1000 s; Pugh et al. 2017) would enable a GRAIL-like scenario at little resource cost to the primary mission but with less accurate satellite-to-satellite ranging accuracy. The CubeSat could be carried to the destination and network with the mothership via a deep-space deployer providing power/thermal during cruise. This CubeSat would be enabled by new technologies in thermal management and propulsion subsystems.

## 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 in this paper, as well as other compelling science questions. Thus, we conclude with the following recommendations:

**Recommendation 1:** Geophysical measurements should be an 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 the variations 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 flexural profiles at various locations. Thermal IR mapping away from the South Polar Terrain is needed to estimate the global heat flux. Laser or radar altimetry could be used to derive a high-resolution shape model and measure tidal deformations.

**Recommendation 2:** A dedicated gravity mapping investigation should be considered for Enceladus. A GRAIL-like mission would be able to measure tidal phase lag and significantly improve the accuracy of the gravity-topography admittance. This would enable investigating the energetics of Enceladus, 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 by leveraging current and upcoming developments in CubeSat technologies and small radios.

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 affine-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 by sampling 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 with more stringent measurement errors. The affine-invariant ensemble sampler uses walkers from the previous step in the

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:

$$ \log L \propto -\frac{1}{2}(\mathbf{X} - \mathbf{Y})^T \Sigma^{-1} (\mathbf{X} - \mathbf{Y}), \quad (A1) $$

where X is the vector of observations and Y is the vector of model predictions, which includes degree 2 gravity and shape coefficients, Love numbers, and libration amplitude at the orbital frequency. Σ is the covariance matrix that contains contributions from the observational and model covariances: $Σ = Σ_{\text{model}} + Σ_{\text{obs}}$. The observation covariance is given by the gravity, shape, tidal Love numbers, and libration determination accuracies. The observations were assumed independent, thus yielding a diagonal $Σ_{\text{obs}}$. Each Markov chain was inspected visually and the initial burn-in steps were discarded. The Markov chains were run until convergence as informed by visual inspection and computed chain autocorrelation times.

### 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-0003-0400-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-0001-5617-207X

Vishnu Viswanathan https://orcid.org/0000-0002-9027-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 Scientific 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
Č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 of Saturn, ed. P. M. Schenk et al. (Tucson, AZ: Univ. Arizona Press), 39

Goodman, J., & Weare, J. 2010, *Communications in Applied Mathematics and Computational Science*, 5, 65

Gurfil, P., Idan, M., & Kasdin, N. 2003, *JGCD*, 26, 491

Gurfil, 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 and the 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, E03003

Hsu, 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, 213

Monteux, J., Collins, G., Tobie, G., & Choblet, G. 2016, *Icar*, 264, 300

Nadezhdina, I., Zubarev, A., Brusnikin, E., & Oberst, J. 2016, *ISPR*, 41, 459
National Research Council 2011, Vision and Voyages for Planetary Science in the 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., Buffington, B. B., et al. 2011, *GeoRL*, 38, L22802

Park, R. S., Bills, B., Buffington, 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 the Icy 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