Gravity Investigation to Characterize Enceladus's Ocean and Interior

A key objective for the future exploration of the icy moon Enceladus is the characterization of the habitable conditions in its internal ocean. Radio science instrumentation on board a spacecraft in a low-altitude orbit about Enceladus would enable gravity measurements that are fundamental to providing constraints on its internal structure. We present here the concept of operations and expected results of the gravity investigation for a New Frontiers–class mission. Numerical simulations are carried out to determine the gravity field in spherical harmonics to degree and order 30 and the Love number k 2 and its phase. By combining Enceladus’s shape measured by Cassini and the geophysical constraints obtained through the processing of the simulated radio science data, a Bayesian inference network is used for the interior model inversion. Our results indicate that the gravity investigation would enable tight constraints on core radius and density, ocean depth and density, and ice shell rigidity. By assuming a high core rigidity and a preliminary modeling of dissipation in the ice shell, our interior model inversion also yields information on the ice shell viscosity. Further data on the hydrosphere properties might be gathered through optical navigation data by accurately measuring Enceladus’s orientation model.


Introduction
The exploration of the icy moons is fundamental to addressing key scientific questions regarding the primordial sources of organics and the existence of habitability conditions elsewhere in the solar system (National Academies of Sciences, Engineering, and Medicine 2022).Enceladus has been identified as one of the best candidates that may currently host sites of abiotic and biological synthesis of organic molecules (e.g., Glein et al. 2015).
The Cassini mission provided multiple lines of evidence for a subsurface liquid water ocean rich in salts and organic compounds (Postberg et al. 2011(Postberg et al. , 2023)).The metabolic energy required to sustain extraterrestrial life underneath Enceladus's ice shell may be sourced from reduction-oxidation (redox) reactions (e.g., Cockell et al. 2016).Chemosynthesis is then expected to be dominant since it may be supported by hydrothermal activity within the seafloor (Choblet et al. 2017).Cassini in situ measurements of plume material support the presence of ongoing hydrothermal processes at the oceancore boundary.Nanometer-scale dust particles consisting of silica, which were detected by the Cassini Cosmic Dust Analyzer (CDA; Srama et al. 2006), may be formed as frozen droplets from the liquid ocean in contact with rock (Hsu et al. 2015).The proportions of dissolved volatiles measured in the plumes by the Ion and Neutral Mass Spectrometer (INMS; Waite et al. 2004) that are consistent with modeling suggest a gas source in Enceladus's ocean (Bouquet et al. 2015).INMS also identified molecular hydrogen (H 2 ) in the plumes during Cassini low-altitude flyby E21 (Waite et al. 2017).This finding further indicates that the subsurface ocean is in contact with the rocky core, leading to the production of this geochemical fuel through serpentinization of chondritic rocks (e.g., Glein et al. 2015Glein et al. , 2018)).The in situ sampling of Enceladus's plumes by CDA and INMS also enabled an accurate characterization of ice grains and water vapor, respectively.
To unambiguously detect and classify the parent organic molecules in Enceladus's plumes, sophisticated mission concepts have been defined by assuming large-or flagshiplike classes (e.g., National Academies of Sciences, Engineering, and Medicine 2022; Voyage 2050 Senior Committee 2021).Scientific instruments developed to search for life in this extraterrestrial environment are designed to analyze the plume materials for evidence of life and biosignature detection (e.g., Reh et al. 2016;MacKenzie et al. 2021).A comprehensive knowledge of the environmental conditions that would support life underneath the ice shell may be obtained through complementary geophysical measurements.An enhanced understanding of the moonʼs interior, including ocean density and the structure and thermal state of the core, is a key objective to investigate habitability (e.g., Vance et al. 2018).
A radio science instrument provides key science measurements that allow the internal structure of Enceladus to be constrained.Radio science observations acquired by the Cassini mission enabled the estimate of Enceladus's mass (e.g., Rappaport et al. 2007) and gravity quadrupole and hemispherical asymmetry (Iess et al. 2014).The mass anomaly identified in the southern hemisphere was initially interpreted as evidence of the presence of a local subsurface water ocean (Iess et al. 2014).By combining gravity and topography data (Nimmo et al. 2011), an estimation of the moment of inertia (MoI) of Enceladus was obtained with a first-order solution by accounting for hydrostatic and nonhydrostatic parts (Iess et al. 2014).A refined estimate of the MoI based on a fourth-order solution (Tricarico 2014) suggests a high degree of internal differentiation and a hydrosphere thickness of ∼50 km (McKinnon 2015).
The analysis of control-point measurements observed with the Cassini Imaging Science Subsystem (ISS; Porco et al. 2004) led to a precise determination of the amplitude of the forced physical libration.This measurement provided evidence of the existence of a global water ocean that decouples the ice shell from the rocky core (Thomas et al. 2016).On the other hand, the properties of the ocean, including density and relative depth from the ice shell, are still unconstrained.Different interpretations of the measured libration only yield consistent information on the ice shell thickness and rigidity (e.g., Thomas et al. 2016;Cadek et al. 2016;Van Hoolst et al. 2016;Hemingway & Mittal 2017).An in-depth gravity investigation would lead to the characterization of the structure and composition of Enceladus's ice shell and ocean through highly accurate measurements of Enceladus's static gravity field and tidal response (i.e., gravitational tide amplitude and phase).
A straightforward recipe is proposed in this study to fill in the knowledge gap on Enceladus's interior after Cassini's outstanding findings.We advocate that the determination of highly accurate long-and short-wavelength gravity anomalies and the Love number k 2 and its phase are fundamental pieces of information missing from a geophysical perspective.A recent analysis of Cassini limb profiles and control points significantly refined Enceladus's topography, leading to a resolution in spherical harmonics to degree 16 (Tajeddine et al. 2017).Highly accurate measurements of Enceladus's gravity field are then required for thorough correlation and admittance analyses with topography to constrain global and local properties of Enceladus's internal structure.
This paper is structured as follows.In Section 2, we describe the proposed orbital mission and radio science instrumentation that would allow us to address outstanding geophysical objectives.In Section 3, we present the results of our numerical simulations based on the concept of operations discussed in Section 2, focusing on the expected accuracy of the gravitational field, tides, and orientation based on the data analysis for precise orbit determination (POD).In Section 4, we discuss the model inversion of Enceladus's interior by using the geophysical constraints obtained through the gravity and radio science investigations.A summary of the mission and instrument concept is reported in Section 5.

Mission Scenario
Enceladus is one of the most coveted places to search for extraterrestrial life in the outer solar system, despite its modest size.The planetary science community is keen to define mission architectures and scientific objectives that will help identify potential biosignatures on the Saturnian moon.A significant advantage that characterizes Enceladus's environment is the presence of plumes that are continuously being released into space and could be directly accessible from an orbiter.In a recent paper, Villanueva et al. (2023) confirmed that the outgassing rate of the plumes that feed the E ring have remained constant since Cassini observations, confirming the potential for future missions to sample material ejected directly from the subsurface water ocean (Reh et al. 2016).A radio science investigation optimizes constraints on the spacecraft trajectory evolution, leading to enhanced mission assurance.An optical navigation camera would provide complementary measurements for a precise onboard localization of the spacecraft during periods that are not covered by radio tracking passes.
An in-depth characterization of Enceladus's internal structure with gravity science data requires a dedicated orbital mission phase with a low-eccentricity orbit around the moon that would be navigable for a required time span.Uniform spatial coverage and temporal coverage are key scientific requirements for this geophysical investigation.
In this work we consider a New Frontiers-class mission with a nominal launch date in 2032 and arrival in the Saturnian system in the early 2040s, after about 9 yr of interplanetary cruise.After spending a few years orbiting Saturn, the solarpowered spacecraft (Reh et al. 2016) is transferred to an orbit about Enceladus in 2048, and the gravity science phase would last for a minimum of 1 month.The orbital parameters of the spacecraft trajectory around Enceladus are modified to target a quasi-circular near-polar orbit, with altitudes between 100 and 150 km.This class of orbits allows for a rather uniform coverage of the underlying gravity field of Enceladus (Figure 1), thanks to the low radial distance, low orbital speed, and steady precession of its orbital plane.Throughout the gravity science phase, Enceladus will complete several revolutions around Saturn (the orbital period is ∼33 hr), enabling the sampling of the satelliteʼs tidal bulge at a wide range of true anomalies.The small eccentricity of Enceladus's orbit gives rise to eccentricity tides that can be described by its tidal Love number k 2 , whose precise measurement (magnitude and phase) is conducive to the characterization of the ocean.
While a near-polar circular orbit is ideal for recovering gravity science data, it is decidedly within a highly unstable orbital regime.Without active orbit maintenance maneuvers (OMMs) to control the orbit, a spacecraft will crash or escape within as little as 1 day (Russell & Lara 2009).Propulsive maneuvers are thus required a few times per Earth day.A preliminary trajectory analysis was carried out to define this orbit profile that enables a trade-off between a uniform spatial coverage and the minimization of OMMs.The designed orbit configuration fulfills the expected engineering and scientific requirements.
The near-circular near-polar orbit analyzed here relies on an onboard autonomous navigation capability.The system would design a correction maneuver for a predefined future opportunity and then provide the maneuver parameters to the flight system to execute at the prescribed time.Autonomous navigation systems can perform the orbit determination and maneuver design in tens of minutes, as opposed to the tens of hours required for a ground-based navigation team.Such systems have significant flight heritage, including Deep Space 1 (Riedel et al. 2000), Deep Impact (Mastrodemos et al. 2005), and OSIRIS-REx (Norman et al. 2022).

Radio Science Instrumentation and Requirements
The spacecraft is equipped with onboard radio science instrumentation that consists of elements that are based on the solid heritage of previous interplanetary missions (e.g., Ciarcia et al. 2013;Asmar et al. 2017), including a transponder that allows for establishing coherent two-way radio links with the stations of NASAʼs Deep Space Network (DSN) and support Telemetry, Tracking and Command (TT&C) activities.The onboard instrument must support dual-frequency links in X band (7.2-8.4GHz) and Ka band (32.5-34.0GHz), for both the uplink and downlink legs.The spacecraft is tracked by the DSN while orbiting Enceladus for 8 hr day −1 , throughout the gravity mapping phase.This schedule is also well suited to meet the spacecraft requirements for energy consumption.The radio links are established through the spacecraftʼs High Gain Antenna (HGA), to guarantee the required signal-to-noise ratio (S/N).We set the overall performance requirement of the onboard radio system at an Allan deviation of 6 × 10 −15 at 1000 s integration time (e.g., Asmar et al. 2005Asmar et al. , 2017)).The requirements on the uplink and downlink S/N are set at 40 dBHz, preventing thermal noise from becoming the leading noise source on the radiometric observables.The link budget will be likely dominated by either plasma or troposphere noise.
A dual-link configuration would lead to the compensation of 75% of the noise associated with charged particles (e.g., Bertotti et al. 1993) in proximity of superior solar conjunctions.The radio science instrumentation will be used during different mission phases that can be near these events when an enhanced characterization of the solar plasma is required for a precise reconstruction of the spacecraft trajectory.Furthermore, Ka band is about 16 times more resistant to plasma noise than X band, making the use of the Ka/Ka link the baseline configuration for science operations.Troposphere noise is strongly dependent on local weather conditions.Sophisticated instrumentation at the ground station with one or multiple Advanced Water Vapor Radiometers (AWVRs) would be capable of calibrating the wet component of the troposphere.The Juno gravity science experiment has shown that troposphere noise can be reduced by as much as 70% with the use of AWVRs at 60 s integration time (Buccino et al. 2021).Currently, only one of the 34 m deep space antennas of the Goldstone DSN complex (DSS-25) is equipped with an AWVR in its proximity.For other stations, the Tracking System Analytical Calibrations (TSACs) are available.Due to Saturnʼs season in the late 2040s, ground stations located in the northern hemisphere are characterized by shorter tracking periods than the ones in the southern hemisphere.Therefore, the preferred daily tracking pass from Goldstone needs to be complemented by additional tracking time at Canberra to meet the 8 hr requirement.We exclude observations that occur at a station elevation angle below 10°, in order to curtail the effect on the radio signal traversing elongated paths into the troposphere.
The predicted noise level on the Goldstone Doppler observables is ∼1.1 mHz for Ka-band uplink and downlink (at 60 s integration time), which corresponds approximately to 0.010 mm s −1 two-way range-rate accuracy.For non-Goldstone passes, the noise increases to ∼2.3 mHz (Ka band, 60 s) or ∼0.021 mm s −1 to reflect the available statistics (e.g., Buccino et al. 2021).For Ka-band range measurements, the utilized noise level is about 20 cm for an integration time of 1000 s (e.g., Genova et al. 2021).Depending on the time of the observations, we expect the presence of spacecraft occultations by Enceladus as seen from Earth.The presence or absence of occultation events is related to the angle between the line-ofsight direction and the spacecraft orbital plane (β Earth angle).A requirement of the orbit profile posed by the gravity investigation is an initial edge-on configuration (β Earth ∼ 0°) that enables the acquisition of Doppler measurements well suited for an accurate determination of the spacecraft trajectory and geophysical parameters, but with significant occultation events.The orbital perturbations, however, cause a precession of the R.A. of the ascending node, leading to β Earth ∼ 30°at the end of the gravity mapping phase.No Saturn occultations are detected during this mission phase.The data percentage lost because of Enceladus occultations is ∼20%.

Scientific Objectives
The gravity and radio science investigation provides information that supplements and contextualizes data on ocean habitability obtained from chemical analysis of plume material.An accurate estimation of the icy moonʼs gravity field, orientation, and tides allows study of Enceladus's properties from the deep interior to the ice shell.Key objectives are the determination of the ocean density and the heat source of the hydrothermal system.

Deep Interior and Ocean
The quadrupole gravity field is a fundamental quantity that leads to constraints on the level of internal differentiation through the measurement of the MoI.The assumption of hydrostatic equilibrium would enable the computation of the dimensionless polar MoI, C MR 2 , by using the Radau-Darwin formalism (e.g., Darwin 1899; Murray & Dermott 2000), where M and R are Enceladusʼs mass and radius, respectively, and C is Enceladus's polar MoI.The analysis of Cassini radio science data led to a separate estimate of the degree-2 spherical harmonic coefficients J 2 and C 22 , supporting significant deviations from hydrostatic equilibrium (Iess et al. 2014).The presence of a nonhydrostatic degree-2 was also detected for the topographic shape (Nimmo et al. 2011).Therefore, a straightforward computation of the MoI with the quadrupole gravity coefficients is not achievable for Enceladus.An alternative procedure was then proposed to interpret the geophysical measurements obtained by Cassini (e.g., Iess et al. 2014;Hemingway & Mittal 2019).The observed degree-2 coefficients for gravity (J 2 obs , ) C 22 obs and topography (H 20 obs , ) H 22 obs are assumed to be given by the linear sum of the hydrostatic and nonhydrostatic parts.By using a fourth-order expansion that is required for a synchronous body (Tricarico 2014), the polar and equatorial eccentricities of the hydrostatic equilibrium shape are computed to determine the hydrostatic parts of both gravity and topography (McKinnon 2015).These terms are then subtracted from the observed quadrupole coefficients to determine the gravity nh and topography (H 20 nh , ) H 22 nh nonhydrostatic components that are used to yield the degree-2 admittances (see Section 3.1).By assuming small lateral variations of the shell mechanical properties (Hemingway et al. 2018), the MoI is derived from the equality condition between the degree-2 admittances of order 0 and 2, Z 20 = Z 22 (for a more detailed discussion see Section 3.1).The resulting uncertainty of the MoI based on this approach depends on both gravity and topography contributions.The current knowledge of Enceladus's MoI is significantly limited by the large uncertainty of J 2 obs and C 22 obs measured after only three Cassini flybys.The proposed orbital gravity mapping configuration and radio science instrumentation are conceived to refine the estimate of the gravity quadrupole that would strongly enhance our knowledge of Enceladusʼs internal mass distribution.The shape model provided by Tajeddine et al. (2017) is well suited to constrain the MoI should the nonhydrostatic gravity term be confirmed.
An independent computation of the mean MoI would be obtained through the measurements of the pole coordinates (R.A., α 0 , and decl., δ 0 ).Enceladus's internal dissipation is expected to provide an equilibrium condition (e.g., Chen & Nimmo 2011), the Cassini state, that results from the coplanarity between the spin axis, orbit plane normal, and Laplace pole (Peale et al. 1980).The gravity phase orbit would allow us to confirm the fulfillment of this equilibrium state by refining Enceladus's ephemeris through radio science data analysis and pole orientation from the joint inversion of radio and optical data.These measurements with the degree-2 gravity field would then be used to constrain the mean MoI (e.g., Baland et al. 2016).This approach, which does not require any additional assumption in addition to the Cassini state, would enable a valid cross-check of the measured MoI based on the gravity degree-2 only.
The combination of Enceladus's MoI and bulk density based on Cassini data analysis enabled preliminary constraints on the rocky core density and radius (e.g., Iess et al. 2014;Beuthe et al 2016;Cadek et al. 2016).An enhanced knowledge of Enceladus's MoI would then allow us to better understand the properties of the rocky core, including the mean radius of the seafloor.To investigate density and thickness of both ice shell and ocean, an additional geophysical measurement is required in the interior model inversion.The Love number k 2 that measures Enceladus's gravitational tidal response to Saturnʼs gravity quadrupole term provides crucial information on the shell thickness and rigidity and constraints on the ocean depth and density.A highly accurate estimate of the amplitude of the gravitational tides is thus a fundamental scientific requirement for the proposed radio science instrumentation.
Radio science data acquired through an orbital mission would enable an accurate k 2 estimate.A joint inversion of interior models with enhanced measurements of Enceladus's bulk density, MoI, and k 2 would better constrain the corehydrosphere density contrast and the thickness of the shell.However, to fully characterize the structure of the hydrosphere, further geophysical information would be required in the inversion.The measurement of the physical librations based on Cassini optical data (Thomas et al. 2016) provides an independent constraint on the ice shell thickness and rigidity (e.g., Van Hoolst et al. 2016).By combining the measured amplitude of the librations with the quantities determined through the proposed investigation, preliminary constraints on the ocean density would be achieved.Furthermore, a navigation camera onboard the orbiter might acquire data during the mapping phase, leading to significant refinements of the orientation model, including the librations (e.g., Konopliv et al. 2018).
A highly accurate knowledge of the gravitational degree-2 and degree-3 would also allow us to use a technique to retrieve further constraints on the hydrosphere.By assuming a threelayer interior modeling (i.e., rocky core, ocean, and ice shell) and internal compensation from elastic forces and isostasy to support the topography nonhydrostatic components, the degree-2 coefficients J 2 and C 22 and the zonal harmonic J 3 might be used to retrieve the ocean and ice shell thicknesses (Hemingway & Mittal 2019).A combined inversion of these quantities determined from the gravity investigation (i.e., J 2 , C 22 , J 3 , bulk density, MoI, k 2 ) would result in a detailed characterization of the hydrosphere properties, including the ocean density (see Section 4.1).
This approach to derive Enceladus's internal properties is based on a Bayesian inference network with constraints on the static gravity and tides.To use these geophysical measurements in the interior model inversion, the core is assumed to be hydrostatic.However, departure from this equilibrium figure has been previously advocated in the literature (e.g., McKinnon 2013; Monteux et al. 2016).Highly accurate measurements of the low-degree gravity coefficients (l = 2-5) might provide evidence of the negligible contribution of the rocky core irregular shape (e.g., McKinnon 2013).If the presence of nonhydrostatic components of core topography is detected, an alternative interior model inversion is a joint analysis of tides and MoI computed with the pole obliquity measured by the gravity investigation and the amplitude of physical librations measured by Cassini (Thomas et al. 2016), or by the combination of radio and optical measurements.This technique is fully independent from assumptions on the core shape and would guarantee a backup solution for the ocean density.These two techniques to yield Enceladus's interior models would strengthen our findings on its core and hydrosphere.

Ice Shell Structure and Rheology
The proposed mission configuration and radio science instrumentation are designed to globally map the icy moonʼs gravitational anomalies at high spatial resolution (<100 km).A gravity mapping phase from a polar or high-inclination orbit is key to enable a uniform latitude and longitude coverage of Enceladus's surface.The high-quality radio tracking data acquired from the NASA DSN stations through a dual-link (X/X, Ka/Ka) configuration are very sensitive to the shortwavelength gravity anomalies (i.e., spherical harmonic degrees l > 10).
The gravity coefficients at degrees higher than 5-10 preserve crucial information on the ice shell structure, including mean thickness and density.Figure 2 shows the expected gravity power spectra that are computed by predicting the signal from the surface and subsurface topography with a Monte Carlo simulation of the parameters' space.By measuring the gravity field to degree and order 16 in spherical harmonics, which is the spectral resolution of the existing topography (Tajeddine et al. 2017), we would investigate the ice shell properties.A global gravity/topography correlation and admittance analysis would yield constraints on the average density of the shell (e.g., Genova et al. 2022).Lateral variations of the shell thickness would also be retrieved through the combination of topography and gravity data (e.g., Wieczorek 2015).
The rheology of the ice shell would be studied by adjusting in the POD process the tidal phase lag, which is related to the tidal Love number k 2 .An accurate estimation of this parameter might provide constraints on the ice rigidity and viscosity.The tidal stresses cause strains that dissipate heat at the bottom of the ice shell.The dissipation integral depends on the ice shell rheology modeling (e.g., Maxwell rheology), including rigidity and viscosity.The combined estimation of the Love number k 2 amplitude and phase lag would enable the computation of the total tidal dissipation within Enceladus (Peale & Cassen 1978).A comparison between the tidal dissipation measured from gravity and that based on the ice shell rheology would then enhance our knowledge of the ice shell rigidity and viscosity.
However, a significant tidal dissipation is also expected in the undifferentiated core even if the interior is frozen (Roberts 2015).Estimates of more than 10 GW have been predicted for the energy dissipated in the unconsolidated core (Choblet et al. 2017).
Figure 3 shows the total power dissipated in the icy moon interior.We computed this budget by using Equation (11) and assuming the Maxwell rheology for both ice shell and rocky core.The relevant rheological parameters were varied in intervals consistent with previous studies of Enceladus tidal dissipation.The ice shell rigidity was varied in the range 2.5-4.5 GPa, while its viscosity is assumed to be uniform throughout the shell and selected from the range 10 12 -10 16 Pa s (e.g., Roberts & Nimmo 2008).We computed the tidal dissipation assuming two different scenarios for the rocky core, including a solid rigid core and an unconsolidated porous core as in Choblet et al. (2017).In the first case, we assumed the rocky core rigidity and the viscosity in the intervals 50-70 GPa and 10 16 -10 22 Pa s (Roberts & Nimmo 2008), respectively.Choblet et al. (2017) characterized the dispersive behavior of the unconsolidated core by assuming the damping ratio measured in laboratory tests along with the shear modulus.For simplicity, here we assume the classical formulation based on the shear modulus, varied in the same range as Choblet et al. (2017), i.e., between 10 7 and 10 9 Pa, and an effective viscosity, for which the interval (10 10 and 10 16 Pa s) was selected to obtain dissipated power in the same range presented by Choblet et al. (2017).This simplified model allows us to preliminarily understand how an accurate measurement of the tidal phase lag may help constrain the internal heat source.High dissipation (∼30 GW) is obtained for extremely low core viscosity (<10 13 Pa s), providing a heat exchange regime that is consistent with models that account for the convection of interstitial water in the porous core (Choblet et al. 2017).This level of power dissipated in the interior is mainly generated by the unconsolidated porous core.A rocky core with viscosities higher than 10 14 Pa s would poorly contribute to the power dissipation, and our constraint on the tidal phase lag would allow us to yield an accurate measurement of the ice viscosity.For this reason, to enable a preliminary determination of the properties of the ice shell, we include in our interior model inversion the measured MoI and the Love number k 2 real and imaginary parts (see Section 4.2) by assuming a high rigidity for the rocky core.Further analyses will be carried out to study the ice rheological properties by accounting for interior models with a less rigid core.
The contribution of the ocean in the total tidal dissipation is mostly considered negligible (e.g., Beuthe 2016; Matsuyama et al. 2018;Rekier et al. 2019), but an independent study suggests the presence of resonantly forced ocean tides that would impact the internal heat budget (Tyler 2020).To investigate the contribution of the ocean, a different modeling of the tides should be applied, and this is beyond the scope of this study.

Gravity Measurements through a Dedicated Mapping Phase
A gravity mapping phase of a New Frontiers-class mission can acquire highly accurate radio science data from a DSN station by establishing X/X and Ka/Ka radio links with a spacecraft in a near-circular near-polar orbit (see Section 2.1).To better define the concept of operations (ConOps) and the scientific requirement of the gravity investigation, we carried out numerical simulations by accounting for the radio subsystem characteristics and environmental constraints.A first step consisted in the simulations of the radio science data by assuming the station DSS-25 at Goldstone as baseline and DSS-43 at Canberra as for additional tracking passes.Because of Saturnʼs season and relative geometry with Earth, the DSS-25 located in the northern hemisphere enables tracking passes on average of 6.5 hr day −1 .To fulfill the surface coverage and 8 hr day −1 tracking requirements, we account for additional data from DSS-43.The noise added to the simulated radio science data is consistent with the measurement requirements described in Section 2.2.
The simulated data are generated and then analyzed through the JPL software MONTE Project Edition (Evans et al. 2018).These measurements are processed through a least-squares filter that yields an adjustment of the parameters of interest by minimizing the discrepancies between simulated and computed data (e.g., Tapley et al. 2004).Our numerical simulations are based on a covariance analysis approach that assumes the same dynamical models for the generation of the simulated data and for the POD process.
The total time span of the gravity mapping phase is divided in our POD analysis into 1-day arcs.This temporal discretization allows us to estimate the state vector of the spacecraft every 24 hr, leading to an overparameterization that mitigates effects associated with propagation errors.To account for errors associated with mismodeling of the nonconservative forces (e.g., orbit correction maneuvers, solar radiation pressure), we include in our estimation process the adjustment of stochastic polynomial accelerations reinitialized every 8 hr with an a priori uncertainty of 10 −10 m s −2 .
In each POD arc we compute the partial derivatives of the local (i.e., state vectors and stochastic accelerations) and global parameters.The latter are the geophysical quantities that are the science requirements of the gravity investigation.We include in our analysis Enceladus's gravitational parameter ( ) GM E and gravity field coefficients in spherical harmonics to degree and order 30, pole coordinates (α 0 , δ 0 ), spin rate (ω), and the Love number k 2 real and imaginary parts.
The information gathered from each arc is then combined in a final inversion through a least-squares filter that provides the joint adjustment of local and global parameters.In this study we report as baseline the gravity inversion based on the analysis of radio science data only.
Preliminary numerical simulation campaigns were carried out by accounting for radio tracking passes from a single ground station during the gravity mapping phase.The results based on this mission scenario suggest a strong degradation (by a factor of ∼20) of the uncertainties of the degree-2 gravity field and tides compared to the case that includes both DSS-25 and DSS-43 enabling ∼8 hr day −1 tracking.For this reason, our nominal scenario is based on the DSS-25 as baseline ground station with an additional 2 hr tracking pass with DSS-43.
A further scenario accounts for a joint processing of radio and optical data acquired during the gravity mapping phase.In our numerical simulations we consider a wide-angle camera with an instantaneous field of view of 1 mrad.The camera is assumed to continuously collect data with an acquisition rate of 600 s.The noise included in our synthetic optical data is modeled as Gaussian with a standard deviation that includes errors associated with the center finding (σ = 0.5 pixels) and attitude (σ = 3 pixels).The optical data are simulated by using a network of known features on Enceladus's surface equally distributed in latitude and longitude.On average approximately four landmarks are detected on each synthetic image.Both radio and optical data are weighted by using the standard deviation of the noise added in the simulated data.

Constraints on Moment of Inertia and Tides
An accurate estimation of the low-degree gravity field is a fundamental objective of the radio science investigation.Table 1 reports the standard deviations of the gravity field coefficients J 2 , C 22 , and J 3 ; orientation parameters; and tides.The gravity field coefficients J 2 and C 22 are key quantities to determine Enceladus's MoI.Observations of the degree-2 of Enceladus's gravity and shape suggest that these signals are both affected by nonhydrostatic contributions, preventing a straightforward computation of the MoI through the Radau-Darwin equation (e.g., Darwin 1899; Murray & Dermott 2000).By assuming that the nonhydrostatic terms are small compared to the hydrostatic terms, the components of the degree-2 gravity and topography can be linearly separated as follows: Hyd  and (b − c)/(a − c) ; 1/4.Therefore, we adopted the methodology reported in Tricarico (2014) to yield the shapes of the coaxial ellipsoids of a two-layer interior model of Enceladus by numerically solving the condition for equipotential surfaces expanded to the fourth order at each interface (further details see McKinnon 2015).This method provides the polar e p and equatorial e q eccentricities of each layer for a given set of semimajor axes and density contrasts.The hydrostatic parts of the quadrupole coefficients of the gravity field are thus computed (Tricarico 2014), where the mass fraction of each layer is M is the body total mass, and a i denotes the semimajor axis of the ith layer.Our model includes two layers, where i = 1 is the hydrosphere and i = 2 is the core.The unnormalized gravity harmonic coefficients C lm,i hyd of each homogeneous triaxial ellipsoid are retrieved according to Yoder (1995), The unnormalized topography harmonic coefficients H lm,i hyd are obtained through the formulation reported in Pěč (1983), is the flattening of the meridional section of the triaxial ellipsoid in whose plane lies the semimajor axis of the ith layer.
The nonhydrostatic parts of the quadrupole coefficients of gravity and topography are computed subtracting the hydrostatic terms from the observed values.The ratio between the nonhydrostatic parts yields a component-wise version of the = -= which depends on the degree and depth of isostatic compensation.By assuming no significant lateral variations of the mechanical properties of the shell, this approach enables the determination of the MoI through the fulfillment of the condition Z 20 = Z 22 (e.g., Iess et al. 2014;McKinnon 2015;Hemingway et al. 2018).The expected MoI for Enceladus is computed with a twolayer model that consists of an outer hydrosphere with density ρ Shell = 925 kg m −3 and a rocky inner core.The density of the core core r is varied between 2000 and 3000 kg m −3 , and the condition Z 20 = Z 22 is checked for each model, similarly to McKinnon (2015).By using the coefficients of Enceladus's gravity field rescaled to a reference radius equal to 252.22 km and the shape coefficients estimated by Tajeddine et al. (2017), we obtain a mean MoI of 0.335 ± 0.00016 (Figure 4).This central value significantly differs from previous studies (e.g., Hemingway et al. 2018) because our computation is based on the latest topography model by Tajeddine et al. (2017) rather than Nimmo et al. (2011).A larger MoI value is more consistent with the first estimate provided by Iess et al. (2014), but our solution accounts for a fourth-order expansion of the equilibrium surface.The MoI formal uncertainty resulting from the adjustment of the degree-2 gravity coefficients in our numerical simulations enables an enhancement of more than an order of magnitude compared to the current knowledge.This level of accuracy is well suited to constrain the degree of differentiation of the internal structure (see Section 4.1).
The method used to determine the MoI by combining gravity and topography degree-2 nonhydrostatic components is based on the assumption that the rocky core is in hydrostatic equilibrium.Constraints on the MoI can also be obtained through the adjustment of the pole coordinates α 0 and δ 0 that are directly tied to the pole obliquity ò.Enceladus is expected to be in a Cassini state (e.g., Chen & Nimmo 2011), which is a necessary condition to determine the MoI from the measured obliquity.Radio science data acquired during mission phases before the gravity mapping would allow us to refine the moonʼs ephemeris, yielding a better understanding of the orbital plane inclination.The Cassini state configuration would then be confirmed by measuring the pole obliquity after the gravity mapping phase.The polar MoI C MR 2 can thus be inferred from the obliquity (e.g., Baland et al. 2016).By adjusting the coordinates of the pole α 0 and δ 0 in our POD solution (Table 1), the resulting uncertainty of the pole obliquity is 0 34 and 0 60 for the case with radio and optical data and radio only, respectively.Depending on the obliquity central value, e.g., ò = 0°.0015 from Chen & Nimmo (2011) and ò = 0°.000356from Baland et al. (2016), 1σ of C/MR 2 may vary between 0.02 and 0.08.This level of accuracy would only provide a cross-validation technique because of the extremely low central value predicted for Enceladus's obliquity.
Another significant contribution of the optical data to the geophysical objectives is the estimation of the amplitude of the physical librations (f).Radio science data are not well suited to fully decouple the effects of spin rate, degree-2 gravity, and librations, since high correlations between these parameters are observed in a global inversion.The combination of optical images that enable the tracking of surface equatorial features and radio science data would allow us to determine the libration amplitude with an accuracy (1σ) of 2 m (Table 1).This measurement would represent a 15 times enhancement compared to current knowledge (1σ ∼ 30 m from Thomas et al. 2016).A better understanding of the amplitude of the physical librations would yield constraints on the polar MoI of the ice shell, enabling accurate measurements of shell average thickness and rigidity (e.g., Van Hoolst et al. 2016).
A key objective of the gravity mapping phase is the Love number k 2 , which plays a crucial role in the interior model inversion to better determine the properties of the hydrosphere.The processing of radio science data results in a 1-standarddeviation formal uncertainty of k 2 equal to ∼1.2 × 10 −4 .By using our estimate of Enceladus's MoI (0.335, Figure 4) and our current knowledge of its mass (1.0802 × 10 20 kg; e.g., Nimmo et al. 2011), we ran Monte Carlo simulations to predict the Love number k 2 amplitude that is consistent with the known geophysical constraints.The internal structure of Enceladus is modeled with three layers that account for the rocky core, ocean, and ice shell.Figure 5 shows a histogram of samples that resembles a χ 2 probability density distribution with a mode of ∼0.02, which is consistent with previous works (e.g., Baland et al. 2016;Ermakov et al. 2021).The shaded area represents the level of uncertainty that would be obtained with the proposed gravity investigation.Furthermore, the analyzed radio science measurements are sensitive to the imaginary part of the Love number k 2 (1σ ∼ 1.7 × 10 −4 ), which is used to better understand the mechanisms of internal power dissipation.

Gravity Field Resolution
A better understanding of the ice shell properties requires gravity and topography maps at the same spatial resolution.The latest estimation of Enceladus's shape leads to a spectral resolution in spherical harmonics to degree and order 16 et al. 2017).As shown in Figure 2, the gravity power spectrum from spherical harmonic degrees l > 7-8 mainly depends on the signal from the surface topography.Our predictions of the gravity from the seafloor topography, which are computed by varying radius, amplitude, and density contrast at core-ocean boundary, suggest a weaker signal compared to the contribution of the surface topography.Gravity/topography admittance analysis to degree 16 would then provide constraints on the ice mean density on global and local scales (e.g., Wieczorek 2015).The combination of gravity and topography is also important to determine the lateral variations of the shell thickness.
Figure 6 shows the power spectra of the gravity from uncompensated and fully compensated topography (shaded gray area) and of the formal uncertainty obtained with radio and optical data (blue) and radio only (red).By using an unconstrained gravity inversion (i.e., no a priori information on the parameters), our results suggest a global resolution in spherical harmonics to degree and order 15-16.This demonstrates the opportunity to use gravity and topography at the same resolution, enabling highly accurate studies on the mean properties of the hydrosphere.
Because of the uneven coverage of Enceladus's surface due to occultation events during the gravity mapping phase, the gravity field of the icy moon would not be measured at the same spatial resolution in each location.By using the degreestrength technique (e.g., Konopliv et al. 1999), predicted gravity accelerations are locally compared to the formal uncertainties of the adjusted accelerations.The intersection between these two curves computed at each latitude-longitude point on the map allows us to determine the local spatial resolution in spherical harmonics.Figure 7 shows the degreestrength map based on the DSN scheduling described in Section 3.This computation is obtained by assuming that the predicted gravity acceleration is consistent with the gravity predicted from partially compensated topography, which resembles a power spectrum that follows a Kaula rule of ´-.The minimum spatial resolution of ∼110 km, which is limited by Enceladus's occultation events, is observed at the equator.The polar region, on the other hand, shows spatial resolution of ∼70 km, which would enable more accurate local gravity analysis in the south polar terrains (SPT).

Interior Model Inversion
The geophysical constraints that result from the gravity investigation are used in a Bayesian inference network to constrain the properties of the internal structure.A Markov Chain Monte Carlo (MCMC) algorithm was developed to invert interior models by varying the free parameters that are adjusted to retrieve the internal properties of the icy moon.This Bayesian inversion enables the exploration of the parameter space with the goal of matching the statistical distribution of the observations (e.g., mass, MoI, tidal Love number).Previous studies used this methodology to investigate the internal structure of celestial bodies (e.g., Khan et al. 2007;Matsuyama et al. 2016Genova et al. 2019;Drilleau et al. 2021;Petricca et al. 2023), enabling significant constraints from the deep interior to the outer shell.
By using random walkers for each free parameter, we sample the multidimensional parameter space by assuming Gaussian distributions of the observables (e.g., Mosegaard & Tarantola 1995).The probability function P(j) of each model j is determined as where ΔO i is the difference between the measured observation O i and its value computed with the parameters of model j, and s is the standard deviation of the observation O i .Our technique is based on the Metropolis-Hastings algorithm (Metropolis et al. 1953;Hastings 1970), which samples a large parameter space by producing interior models in chains.Each chain is started by drawing a sample from the prior distributions, assumed to be uniform for each of the parameters across the ranges that we consider (see Table 2-parameter space).The proposal distributions used to perturb each model and to generate the following model in the chain are assumed to be Gaussian with standard deviation (see Table 2-random walker step size) selected by tuning the algorithm on the observations that are accounted for in our inversion.
To properly map the parameter space and to reduce the dependence of the inversion results on the initial model, 20 independent and distinct chains are generated.The convergence of the chains is tested with the Gelman-Rubin criterion (Gelman & Rubin 1992) by selecting a threshold for the potential scale reduction factor (PSRF) of 1.01 (for details see Petricca et al. 2023).Each chain is iterated until ∼50,000 models are accepted.By discarding a first set of 10,000 models (i.e., burn-in) that allows each chain to reach the targets, we obtain a PSRF lower than the threshold for each observation.The results are presented in terms of the mode of the posterior distributions, and the associated credible intervals are estimated through the highest probability density method (Chen & Shao 1999).
A three-layer model is included in our MCMC algorithm by assuming the presence of the core, ocean, and ice shell.The properties of each layer vary in one dimension only.To test and validate our methodology, we carried out an interior model inversion by using the two observations obtained by Cassini.We consider Enceladus's mass M = The probability distribution resembles a χ 2 distribution with a mode of 0.02, which is consistent with previous studies (e.g., Baland et al. 2016;Ermakov et al. 2021).The shaded area shows the level of uncertainty that is expected with the gravity investigation (this study, see Table 1).
(1.08022 ± 0.00108) × 10 20 where the standard deviation represents a conservative confidence level, and MoI = 0.335 ± 0.002 computed by using gravity coefficients scaled on the new reference radius (Hemingway et al. 2018) and topography (Tajeddine et al. 2017), as described in Section 3.1.The free parameters investigated in our MCMC algorithm for this case are the core density and size, ocean density and size, and ice shell thickness.
In our MCMC algorithm we evaluate the icy moonʼs mass M and MoI by accounting for the ellipsoid of each layer (i = 1 core, i = 2 ocean, and i = 3 ice shell).By considering the core radius r core , the ocean depth d ocean and Enceladus's mean radius R E , the semimajor axes of the core, ocean, and ice shell are a r , and a 3 = R E , respectively.We then compute polar and equatorial eccentricities for each layer through the fourth-order solution of the equilibrium figure theory (Tricarico 2014).This allows us to determine the polar (c i ) and equatorial (b i ) semiminor axes of each ellipsoid.By summing the contribution of each layer, the mass and MoI = C MR E 2 are computed through the following analytical expressions: where ρ i is the mean density of the ith layer.The axes a 4 , b 4 , and c 4 included in the sum are null.Table 2 reports the random walker step size and range of the free parameters.Figure 8     & Mittal 2017).Our technique provides more conservative uncertainties compared to parametric analysis by accounting for the correlations among the observations.Cassini measurements of the icy moon mass and MoI also enable an enhanced knowledge of the hydrosphere thickness, but these data do not provide any constraint on the ocean properties.The uniform distribution of the ocean density in Figure 8(d) strongly suggests that further geophysical data are required to characterize the hydrosphere.
The gravity investigation proposed in this study is conceived to accurately determine the low-degree spherical harmonic coefficients (J 2 , C 22 , and J 3 ), and the Love number k 2 amplitude and phase lag to better understand the properties of the ocean (Section 4.1) and ice shell (Section 4.2).

Ocean Density
An accurate knowledge of the ocean properties is one of the main objectives of the proposed gravity investigation.To determine the density of the internal ocean, our model inversion is carried out by combining constraints on the internal mass distribution (mass and MoI) with key measurements related to isostatic/flexural support of nonhydrostatic topography (low-degree gravity harmonics) and tidal response (Love number k 2 ).The modeling of mass and MoI in the MCMC is based on Equation (8).We describe below the approach implemented in our algorithm to include the odd zonal harmonic gravity coefficient that relies on the topography support and the tidal Love number k 2 as observations of the Bayesian inversion.
The nonhydrostatic surface topography is expected to be supported isostatically through Airy or Pratt compensation, or by the elastic flexure, or by a combination of these end-member mechanisms (Hemingway & Mittal 2019).Our method is tailored to deal with an unconstrained degree of compensation.By assuming equal pressures, the gravity/topography admittance is given by (e.g., Turcotte 1981; Hemingway & Matsuyama 2017) where ρ ice = 925 kg m −3 is the ice shell density, T c is the mean shell thickness, R E is Enceladus's mean radius, C l 0 is the compensation factor that is a function of the effective elastic thickness T e (see Equations ( 14 are the gravity accelerations at the top and bottom of the shell, respectively.By using the nonhydrostatic topography coefficient H 30 nh computed from Enceladus's measured topography and the model-based admittance (Equation ( 9)), we derive the zonal harmonic coefficient in our MCMC algorithm, Since the proposed gravity investigation would enable the determination of the gravity field at the same resolution of the existing topography (e.g., Tajeddine et al. we expect that local and global gravity/topography admittance analyses to degree and order 16 would allow us to constrain the elastic thickness, T e .For this reason, in our preliminary interior model inversion the elastic thickness is not adjusted, but it is assumed to be T e = 0, which is consistent with an Airy compensation scenario (compensation factor C l 0 = 1).By assuming a different a priori degree of compensation ( ) C 1 l 0 ¹ in our analysis, we only note differences in the mean values of the posterior distributions of the adjusted parameters but no variations in the resulting uncertainties.
To determine the gravitational tidal response for each interior model, we integrated the ALMA3 software (Spada 2008;Melini et al. 2022) MCMC algorithm.ALMA3 has been already employed to compute the Love numbers of multiple celestial bodies, including Mercury (e.g., Goossens et al. 2022), Venus, and Mars (e.g., Petricca et al. 2022).This further plugin allows us to retrieve the real and imaginary parts of the Love number k 2. ALMA3 accounts for concentric spherical shells.A Maxwell rheology is used for both ice shell and core, whose average viscosity and rigidity are treated as free parameters of the MCMC.The ice shell is characterized by a rigidity μ Shell varying between 2.5 and 4.5 GPa and viscosity ν Shell in the range 10 12 -10 16 Pa s.The fluid behavior of the ocean is approximated by treating this layer as a low-viscosity material (10 3 Pa s) with negligible rigidity.This approximation is valid as long as the Maxwell time of the layer (τ M = η 0 /μ 0 ) is far from the forcing period (Roberts & Nimmo 2008).The core of Enceladus is modeled as a viscoelastic body with shear modulus Core m allowed to range between 50 and 70 GPa and viscosity Core n in the range 10 16 -10 22 Pa s.The Bayesian inversion is thus carried out by using the following observations: M = 1.080 22 ± 0.00108 × 10 20 kg, MoI = 0.335 ± 0.00016, J 3 = −118.0± 0.24261 × 10 −6 , and Love number Re(k 2 ) = 0.02 ± 3.5 × 10 −4 and Im(k 2 ) = 0.01 ± 1.7 × 10 −4 .Enceladus's mass and MoI constraints are consistent with the case presented in Section 4 by only updating the MoI standard deviation with the expected level of accuracy from our gravity investigation (Section 3.1).The central values of J 3 and k 2 are in agreement with estimates from Cassini data (Hemingway et al. 2018) and predictions from theoretical models (see Figure 5), respectively.The imaginary part of the Love number is consistent with a dissipated power of ∼15 GW and the heat flux observed at the SPT (Howett et al. 2011).This is equivalent to assuming that Enceladus is in thermal equilibrium, with the tidal dissipation balancing the power conductively transported through the ice shell and radiated to space.This equilibrium configuration is consistent with recent astrometric observations of the tidal dissipation within Saturn (Lainey et al. 2012(Lainey et al. , 2017)).
The standard deviations of the observations result from our expected gravity uncertainties.Figure 9 shows the distributions of the ensemble for the observation included in our interior model inversion.The histograms of the inferred characteristics of Enceladus's rocky core, including density and radius, and of the ocean thickness and density are reported in Figure 10.The results of the MCMC algorithm based on the mass, MoI, k 2 , and J 3 suggest that our gravity investigation would significantly enhance our understanding of the core properties (i.e., core radius r 196

Internal Energy Dissipation
The imaginary part of the Love number k 2 ( ) k Im 2 quantifies the lag of the tidal response with respect to the forcing, and it is related to the total energy dissipation rate E  within the body (Segatz et al. 1988) where n is the mean motion of Enceladus, e its orbital eccentricity (e = 0.0045, Roberts & Nimmo 2008), and R E its radius.Since the tidal lag depends on the material viscosity and rigidity, an accurate measurement of this parameter informs on the mechanisms that lead to the internal energy dissipation.
The sensitivity of the Love number k 2 real and imaginary parts on the density and rheological structures enables constraints on the properties of the ice shell.Figure 11 shows the distribution for the core and ice shell rigidities resulting from our inversion.Given the large sensitivity of ( ) k Re 2 on the rigidity of the ice layer and the highly accurate estimate of this parameter with the proposed investigation, we expect a constraint on the ice rigidity of ∼0.4 GPa.Our inversion also indicates that the geophysical observations are well suited to determine the viscosity of the ice shell with a ( ) log 10 ice n accuracy of 0.04.The core viscosity is unconstrained, suggesting that our results enable unambiguous information on the internal energy dissipation in the ice shell only.The ice viscosity uncertainty, however, obtained with this preliminary approach is based on two simplified assumptions.As described in Section 2.3.2 and reported in Table 2 (range of explored core rigidities, ) core m , our analysis is based on the hypothesis that the core is rigid.Our modeling also accounts for an internal dissipation in the whole ice shell, although heat is expected to be dissipated at the bottom of the shell only.

Summary
Gravity and radio science investigations conducted through a spacecraft in orbit about Enceladus can address outstanding questions regarding the icy moonʼs internal structure.In this study, we presented a gravity mapping phase of a New Frontiers-class mission with an onboard radio science instrumentation that supports a dual-frequency configuration.Highly accurate Doppler measurements enable a precise characterization of the long-wavelength gravitational anomalies and tides, providing key constraints on the icy moonʼs rocky core and hydrosphere.The dedicated radio science campaign through an orbital phase is fundamental to accurately estimate the gravity field coefficients to degree and order 16 that jointly with the known surface topography inform on the lateral variations of the ice shell thickness and density.
A crucial requirement of the radio science measurements is thus the uniform spatial and temporal coverage of Enceladus's surface, yielding constraints on its gravity field and tides.By assuming ∼8 hr day −1 radio tracking passes with DSN stations, we expect spectral resolutions of global and local gravity field in spherical harmonics to degree 16 (see Figure 6) and 22 7), respectively.This enhanced knowledge of the gravity anomalies would be fully consistent in terms of resolution with the latest estimate of Enceladus's topography (Tajeddine et al. 2017).A 1-month time span of the gravity mapping phase is well suited to precisely detect and measure the signal associated with the degree-2 gravitational tides.Our numerical simulations suggest high precisions for the Love number k 2 and its phase.
A preliminary inversion of Enceladus's interior models is presented in this study to better understand the constraints on the internal properties resulting from the measured geophysical quantities, including MoI, low-degree gravity coefficients, and k 2 .Our approach is based on a Bayesian inference network that accounts for 1D models of Enceladus's interior by assuming three concentric layers (i.e., rocky core, ocean, and ice shell).Our is embedded with ALMA3 (Melini et 2022) to determine the viscoelastic response of the internal layers that affect the gravitational tides.A joint analysis of the geophysical constraints enables an accurate characterization of the properties of the core, ice shell, and ocean.
Ice shell thickness, ocean density, core density, and source of tidal heating are all key parameters that complement direct sampling of Enceladus's plume as a record of the ocean composition.Together they provide the most detailed picture of Re k 2 , and ( ) k Im 2 resulting from the numerical simulations of this study (see Table 1).the habitability of liquid water ocean obtainable for any body other than Earth.

Figure 1 .
Figure 1.Body-fixed 3D (left) and polar (right) views of the spacecraft orbit evolution during the gravity mapping phase.

Figure 2 .
Figure 2. Power spectra of the gravity field computed from isostatically compensated and uncompensated topography (gray shaded area) and seafloor topography by varying its mean radius ( ) R sf , amplitude (Δh), and density contrast with the ocean (Δρ).Each panel shows with colors the power spectrum of the gravity field from the seafloor topography by assuming constant values for Δρ = 1500 kg m −3 and Δh = 5 km (left), R sf = 190 km and Δh = 5 km (middle), and Δρ = 1500 kg m −3 and R sf = 190 km (right).

Figure 3 .
Figure 3. Power dissipated in Enceladus's interior as a function of ice shell and rocky core viscosity by assuming high rigidity (left) and low rigidity (right) for the rocky core.
For synchronous moons significant deviations from the firstorder solution are observed for the ratios of the hydrostatic gravity terms and of the semiaxes of the hydrostatic equilibrium ellipsoidal figure, i.e.,

Figure 4 .
Figure 4. Degree-2 admittances Z 20 (blue) and Z 22 (red) as a function of the MoI.A two-layer model is assumed in this computation by accounting for the fourth-order expression for equilibrium figures (Tricarico 2014).Ice density is considered fixed (ρ Shell = 925 kg m −3 ), and we only vary the core density Core r between 2000 and 3000 kg m −3 .The intersection of the admittances (Z 20 = Z 22 ) identifies the expected value of MoI = 0.335 that provides selfconsistent results (e.g., Hemingway et al. 2018).The admittance uncertainties (blue and red shaded areas) are retrieved by assuming 3σ uncertainties for the shape (Tajeddine et al. 2017) and gravity models (this paper).Gray and yellow shaded areas show the level of accuracy for the MoI based on the gravity estimates obtained with Cassini (σ MoI = 0.002) and this study (σ MoI = 0.00016), respectively.

Figure 5 .
Figure5.Histograms for the Love number k 2 predicted through the MCMC interior model inversion based on the current knowledge of Enceladus's mass and MoI.The probability distribution resembles a χ 2 distribution with a mode of 0.02, which is consistent with previous studies (e.g.,Baland et al. 2016;Ermakov et al. 2021).The shaded area shows the level of uncertainty that is expected with the gravity investigation (this study, see Table1).
shows the results of our MCMC algorithm with Cassini observations of Enceladus's mass and MoI.The core radius r m −3 are well constrained and consistent with previous studies (e.g.,Iess et al. 2014;McKinnon 2015;

Figure 6 .
Figure 6.Power spectra of the predicted gravity field based on isostatically compensated (gray dashed line) and uncompensated (gray solid line) topography (Tajeddine et al. 2017) and of the formal uncertainties of the gravity field determined with radio and optical data (blue dashed line) and radio only (red solid line).The black line shows the Kaula rule 40 × 10 −5 /l 2 that resembles the power spectrum of a partially compensated topography.

Figure 7 .
Figure7.Degree-strength map of Enceladus's gravity field obtained by using the covariance matrix from the numerical simulations and the expected acceleration profile based on the gravity from uncompensated topography.

Figure 8 .
Figure 8. Histograms for core (a) radius and (c) density and ocean (b) depth and (d) density resulting from the MCMC carried out by using as observations Enceladus's mass and MoI measured by Cassini.

+
m −3 compared to Cassini observations.The joint analysis of k 2 and J 3 , which are sensitive to the thickness of the hydrosphere and to the ocean density, is fundamental to decouple icy shell and ocean properties.The ocean depth is constrained to d kg m −3 , which shows a level of uncertainty of ∼30 kg m −3 .

Figure 9 .
Figure 9. Observation (black) and mapped (red) distributions and histograms of the ensemble for Enceladus's mass, MoI, J 3 , and Love number Re(k 2 ) and Im(k 2 ).

Figure 11 .
Figure 11.Histograms for core and ice shell (a, b) rigidities and (c, d) viscosities resulting from the MCMC interior model inversion.

Table 1
Current Knowledge and 1-standard-deviation Predicted Formal Uncertainties of the Low-degree Spherical Harmonic Unnormalized Coefficients (J 2 , C 22, J 3 ), R.A. Amplitude of the Physical Librations (f), and Love Number k 2 Real and Imaginary Parts 0 d of the Pole, Spin Rate ( ) W  ,

Table 2
Random Walker Step Sizes and Boundaries of the Free Parameters Explored in the MCMC Interior Model Inversion