A Recipe for Geophysical Exploration of Enceladus

Orbital geophysical investigations of Enceladus are critical to understanding its energy balance. 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 a 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.


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 , 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.
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 2013McKinnon , 2015Beuthe 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;Hemingway & Mittal 2019;Čadek et al. 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 spatio-temporal 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 (Spitale et al. 2015;Běhounková 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 Fig. 1.
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 2007Bě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:

Quantity Current Knowledge Information Provided
Gravity field Up to degree 3 field measured by Cassini (Iess et al. 2014).
Constraint on the internal structure and compensation mechanism in the ice shell.
Constraint on the internal structure and ice shell thickness variations.
Constraint on compensation state of topography, total shell thickness, elastic shell thickness and shell density and shell-ocean density contrast. Libration amplitude is strongly sensitive to the shell thickness. Observed large amplitude requires a decoupled ice shell and hence implies a global subsurface ocean.
Precession and nutation Not currently measured. Precession rate estimated 2.6 rad/yr ).
In combination with degree 2 gravity, enables moment of inertia determination without the need for hydrostatic equilibrium assumption.
Provides a constraint on 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 shell and its rigidity. Imaginary parts are sensitive to the viscosity profile and, thus, provide a constraint on total tidal dissipation within Enceladus.
Radar sub-surface mapping Not currently achieved. Enables independent shell thickness determination.

Magnetic induction
Not currently measured. The magnitude of the forcing field is ≈10 nT.
Sensitive to the sub-surface ocean conductivity and thickness.
Constraint on the total heat flux. † Gravity-topography admittance is not an independent geophysical observable. It is derived as the ratio of the gravity amplitude to topography amplitude. 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 (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). 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 non-hydrostaticity 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 (Thomas et al. 2016;Nadezhdina 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 meters (Park et al. 2020a)-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-meter 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 2 O, salt, organics, and other volatiles indicating a connection to a reservoir of ocean-derived material . Regional gravity and topography mapping would allow to characterize the crustal structure in the South Polar Terrain. 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 be enable assessing the shell viscosity structure (Akiba et al. 2021) 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 non-hydrostaticity between the shell and the core -Separating the non-hydrostatic 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 to identify tectonics and convection zones allowing to separate the signals originating within the ice shell. A direct detection of the ice-ocean interface would allow subtracting the respective gravity signal of the ice shell and infer 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 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 gravitytopography 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), where n is degree and m is order. For a spherically symmetric body, the Love numbers are degenerate with respect to m. Thus, separating tidal signal corresponding to different orders within the same degree would indicate deviations of the body 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 two due to the nonuniform shell thickness and faults in the South Polar Terrain.

Where is the heat generated and how is it transported?
Tidal heating is thought to provide 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) onboard 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)  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). 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 ( 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 iceocean interface and return the shell thickness with an accuracy better than 100 m. This would constrain the regional dynamics and the 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 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.

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 as small as ≈2,000 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 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 high 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 Eq. 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 years 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.

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 Gy, Hansen et al. 2006) and also offers the opportunity to sample material from Enceladus' ocean. The rate of surface-interior exchange is important for sustaining habitability (e.g., Soderlund et al. 2020). 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), nor are the resurfacing mechanisms or rates of resurfacing and surface-interior exchange well constrained (Spencer et al. 2018;Bland et al. 2015). 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 the 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 since a significant advection of heat occurs in the fractures. Better constraints on the spatio-temporal pattern of eruptions 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 leave dielectric signatures leading to strong reflectors in the 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. Vertical fractures could be identified by point scattering usually leaving characteristic hyperbolas in the (non-focused) radar return. Due to the complex surface topography of Enceladus, the radar sounding data has 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 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) measurementsanalogous 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).

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 3.
To generate synthetic observations, we compute a multi-layer hydrostatic equilibrium model using Tricarico (2014). Viscoelastic moduli and densities are assumed constant within each layer. This gives us the hydrostatic shape and  Table 2. Prior probability distribution of Enceladus internal structure parameters. 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.
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 non-hydrostatic 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 10 9 and 5 · 10 9 Pa leads to a libration amplitude difference of 50 meters, which is just below current observational uncertainty. The variation of libration amplitude due to the shell viscosity profile is yet smaller: < 20 meters for bottom shell viscosity ranging from 10 12 to 10 15 Pa s ). More details on the MCMC geophysical inversion are given in Appendix A.
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 , l 2 ; and the gravitytopography admittance spectrum Z l , where l is the spherical harmonic degree, 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. Fig. 2 shows the posterior distributions for selected internal structure parameters in form of a corner plot (Foreman-Mackey 2016). The top boxes show the 1D histograms of parameters of interest. In addition, we show 95% confidence regions for selected parameter pairs. The black contours in Fig. 2 represent the current state of knowledge. The other (colored) contours show how the confidence regions shrink as more data is 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 meters 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 Figure 2. Corner plot visualizing the posterior distribution of internals 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 95% 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. 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.

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 Fig. 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. (2020a) as: 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.5%-97.5% percentile ranges are given at the top. max(∆R) = 12 2 7 h 2 R 2 ω 2 e g (1) max(∆N ) = 12 l 2 R 2 ω 2 e g (2) max(∆E) = 9 l 2 R 2 ω 2 e g for the radial, northerly and easterly displacement, respectively. In addition, the maximum surface gravity variation range is given by: Taking the median values for the Love numbers from Fig. 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.9 and 0.6 meters in the radial, northerly and easterly directions. The map of tidal displacement ranges for these median values of Love numbers is shown in Fig. 4. In addition, Fig. 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 tidal displacements are maximized at the equator. Figure 4. Radial, northerly and easterly tidal displacement ranges as well as 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).
The maximum radial tidal deformation range is 0.5-6.5 meters (2.5-95% 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 meters 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 sub-meter accuracy on the Enceladus orientation.
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): Thus, the measurement requirement on the Im(k 2 ) should be tied to the requirement of on the heat flux measurement. Fig. 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 constrain 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).

Gravity and tides sensitivity
We have studied the sensitivity of single spacecraft with Earth-based Doppler tracking and dual spacecraft with inter-spacecraft 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. 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 relative motion depends on the relative position of the spacecraft. Assuming that a formation of two spacecraft co-orbits Enceladus at the same mean altitude along circular equatorial orbits, the design space reduces to choosing their angular separation and orbit altitude. 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 spacecraft's mean-motion and tidal harmonics.
The amplitude of the tidal deformation and tidal gravity perturbation is at the maximum at the equator (Fig. 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.
If the body is spherical, the Love numbers are degenerate with respect to spherical harmonic order m. Thus, for such a body: k 20 = k 21 = k 22 . Deviations from spherical symmetry would break this degeneracy. Běhounková et al. (2017) estimated that Enceladus' degree 2 potential Love numbers could be different by a factor of two. In order to separate the effects of the degree 2 Love numbers, a non-equatorial 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.

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 (Lidov 1962;Kozai 1962) produces a long-period exchange between inclination and eccentricity, effectively increasing the eccentricity of the orbit while reducing its inclination and vise 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: 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°and i max = 140.8°. 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, Eq. (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 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. 2020b). 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 ∆V cost of 20 m/s 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 ensuring 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 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.

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 the Earth and dual spacecraft (GRAIL-like configuration) with inter-satellite range-rate measurements. The covariance analysis is based on a least-squares principal and is a powerful tool for assessing the expected uncertainties in the estimated parameters (Park et al. 2011(Park et al. , 2012(Park et al. , 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, 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 ≈ 5 · 10 −7 , 5 · 10 −8 and 5 · 10 −9 over a 10-s integration time (Bills & Ermakov 2019;Asmar et al. 2005). The impact of the frequency used for spacecraft Doppler tracking depends on the signal-to-noise ratio 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 the 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 km. 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 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 since 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 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 to 10 measurement accuracy 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 Fig. 6. In addition, Fig. 6 shows the expected amplitude of the gravity signal as well as the spatial scales of several geologic features of interest. Figure 6. Summary of the gravity anomaly signals due to various geologic landmarks. The vertical axis shows the Root-Mean-Square (RMS) magnitude of the gravity signal. The top horizontal axis shows 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 fitted to an uncompensated gravity-from-shape using the shape model by Tajeddine 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 a shell density of 600 kg/m 3 . The large scale depressions refer to the chain of basins identified by Tajeddine et al. (2017).
Assuming a single orbiter with two-way X-band tracking capability, and assuming the Doppler data accuracy of 10 −7 km/s at a 60-s count time with a 28-day data collection duration, a gravity field up to degree ≈10 can be recovered and k 2 can be recovered with an accuracy of ≈ 10 −2 . A high-resolution context imager with accurate pointing knowledge and image-motion-compensation capability, such as the Advanced Pointing Imaging Camera (Park et al. 2020a) (APIC), can measure h 2 and l 2 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 h 2 , l 2 , and libration amplitude can be improved with 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 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 measurement of radial tides is dominated by interpolation errors at cross-over 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. Since the number of cross-over points increases with n 2 orb , millions of cross-over 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 4-5 better than the inherent range resolution can be achieved (Steinbrügge et al. 2018).
A dual spacecraft architecture with inter-satellite 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·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 k 2 error down to ≈ 2 · 10 −5 (or relative accuracy of ≈0.1%, Table 3). The gravity field can be determined to degree 20 to 30 (see Fig. 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 −14 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.

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 GRAILlike 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.