Gravity Investigation to Characterize Enceladus's Ocean and Interior

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 solar-powered 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₂, whose precise measurement (magnitude and phase) is conducive to the characterization of the ocean.

Zoom In Zoom Out Reset image size

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

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.4 GHz) and Ka band (32.5–34.0 GHz), for both the uplink and downlink legs. The spacecraft is tracked by the DSN while orbiting Enceladus for 8 hr day⁻¹, 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⁻¹⁵ at 1000 s integration time (e.g., Asmar et al. 2005, 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⁻¹ 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⁻¹ 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-of-sight 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%.

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.

2.3.1. Deep Interior and Ocean

2.3.2. 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 short-wavelength 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).

Zoom In Zoom Out Reset image size

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₂. 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₂ 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¹²–10¹⁶ 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¹⁶–10²² 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⁷ and 10⁹ Pa, and an effective viscosity, for which the interval (10¹⁰ and 10¹⁶ 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¹³ 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¹⁴ 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₂ 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.

Zoom In Zoom Out Reset image size

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.

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₂, C₂₂, and J₃; orientation parameters; and tides. The gravity field coefficients J₂ and C₂₂ 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:

$$\begin{eqnarray}&&\begin{array}{l}{J}_{2}^{\mathrm{obs}}={J}_{2}^{\mathrm{hyd}}+{J}_{2}^{\mathrm{nh}}\\ {C}_{22}^{\mathrm{obs}}={C}_{22}^{\mathrm{hyd}}+{C}_{22}^{\mathrm{nh}}\\ {H}_{20}^{\mathrm{obs}}={H}_{20}^{\mathrm{hyd}}+{H}_{20}^{\mathrm{nh}}\\ {H}_{22}^{\mathrm{obs}}={H}_{22}^{\mathrm{hyd}}+{H}_{22}^{\mathrm{nh}}.\end{array}\end{eqnarray} \tag{ 1 }$$

Table 1. Current Knowledge and 1-standard-deviation Predicted Formal Uncertainties of the Low-degree Spherical Harmonic Unnormalized Coefficients (J₂, C_22, J₃), R.A. $\left({\alpha }_{0}\right)$, and decl. $\left({\delta }_{0}\right)$ of the Pole, Spin Rate $\left(\dot{W}\right)$, Amplitude of the Physical Librations (ϕ), and Love Number k₂ Real and Imaginary Parts

Parameter	Current Knowledge		Radio Science	Radio Science and Optical
	Value ± 1σ	References
J2 (×106)	5526.1 ± 35.5	Hemingway et al. (2018)	1.95 ×10−2	1.94 ×10−2
C22 (×106)	1575.7 ± 15.9		7.41 ×10−3	7.10 ×10−3
J3 (×106)	−118.2 ± 23.5		2.39 ×10−2	2.35 ×10−2
α0 (deg)	40.66	Archinal et al. (2018) a	2.4 ×10−3	1.2 ×10−3
δ0 (deg)	83.52		1.5 ×10−4	1.1 ×10−4
$\dot{W}\ (\deg \ {s}^{-1})$	262.732		4.5 ×10−10	4.5 ×10−10
ϕ (deg)	0.120 ± 0.007	Thomas et al. (2016)	⋯	4.2 ×10−4
Love number	${Re}({k}_{2})$	⋯	1.2 × 10−4	9 × 10−5
	$\mathrm{Im}({k}_{2})$	⋯	1.7 × 10−4	1.2 × 10−4

Note.

^a Formal uncertainties of the pole coordinates are not included in the IAU report.

Download table as:  ASCIITypeset image

For synchronous moons significant deviations from the first-order solution are observed for the ratios of the hydrostatic gravity terms and of the semiaxes of the hydrostatic equilibrium ellipsoidal figure, i.e., ${J}_{2}^{\mathrm{hyd}}/{C}_{22}^{\mathrm{Hyd}}\simeq 10/3$ 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),

$$\begin{eqnarray}&&{C}_{\mathrm{lm}}^{\mathrm{hyd}}={\sum }_{i=1}^{2}{\left(\tfrac{{a}_{i}}{{a}_{1}}\right)}^{l}\left(\tfrac{{M}_{i}}{M}\right){C}_{\mathrm{lm},{\rm{i}}}^{\mathrm{hyd}},\end{eqnarray} \tag{ 2 }$$

where the mass fraction of each layer is

$$\begin{eqnarray}&&\tfrac{{M}_{i}}{M}=\tfrac{{V}_{i}\left({\rho }_{i}-{\rho }_{i-1}\right)}{M}=\left(\tfrac{4\pi }{3}\right){a}_{i}^{3}\sqrt{1-{{e}_{p}}_{i}^{2}}\sqrt{1-{{e}_{q}}_{i}^{2}}\tfrac{\left({\rho }_{i}-{\rho }_{i-1}\right)}{M},\end{eqnarray} \tag{ 3 }$$

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}_{\mathrm{lm},{\rm{i}}}^{\mathrm{hyd}}$ of each homogeneous triaxial ellipsoid are retrieved according to Yoder (1995),

$$\begin{eqnarray}&&\begin{array}{c}{C}_{20,i}^{\mathrm{hyd}}=-{J}_{2,i}^{\mathrm{hyd}}\ =\tfrac{\left({e}_{{q}_{i}}^{2}-2{e}_{{p}_{i}}^{2}\right)}{10},\\ {C}_{22,i}^{\mathrm{hyd}}=\tfrac{{e}_{{q}_{i}}^{2}}{20}.\end{array}\end{eqnarray} \tag{ 4 }$$

The unnormalized topography harmonic coefficients ${H}_{\mathrm{lm},{\rm{i}}}^{\mathrm{hyd}}$ are obtained through the formulation reported in Pěč (1983),

$$\begin{eqnarray}&&\begin{array}{c}{H}_{20,i}^{\mathrm{hyd}}=-{a}_{i}\left(\tfrac{2}{3}{\alpha }_{i}+\tfrac{1}{7}{\alpha }_{i}^{2}-\tfrac{1}{3}{\alpha }_{1}-\tfrac{1}{7}{\alpha }_{i}{\alpha }_{{1}_{i}}-\tfrac{1}{3.7}{\alpha }_{i}^{3}\right),\\ {H}_{22,i}^{\mathrm{hyd}}=\tfrac{1}{6}{a}_{i}\left({\alpha }_{{1}_{i}}-\tfrac{3}{7}{\alpha }_{i}{\alpha }_{{1}_{i}}\right),\end{array}\end{eqnarray} \tag{ 5 }$$

where ${{\alpha }_{1}}_{i}=\tfrac{{a}_{i}-{b}_{i}}{{a}_{i}}$ is the equatorial flattening and ${\alpha }_{i}=\tfrac{{a}_{i}-{c}_{i}}{{a}_{i}}$ (a_i > b_i > c_i ) 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 spectral admittance,

$$\begin{eqnarray}&&\begin{array}{c}{Z}_{20}=-\tfrac{{J}_{2}^{\mathrm{nh}}}{{H}_{20}^{\mathrm{nh}}},\\ {Z}_{22}=\tfrac{{C}_{22}^{\mathrm{nh}}}{{H}_{22}^{\mathrm{nh}}},\end{array}\end{eqnarray} \tag{ 6 }$$

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₂₀ = Z₂₂ (e.g., Iess et al. 2014; McKinnon 2015; Hemingway et al. 2018).

The expected MoI for Enceladus is computed with a two-layer model that consists of an outer hydrosphere with density ρ_Shell = 925 kg m⁻³ and a rocky inner core. The density of the core ${\rho }_{\mathrm{core}}$ is varied between 2000 and 3000 kg m⁻³, and the condition Z₂₀ = Z₂₂ 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).

Zoom In Zoom Out Reset image size

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 α₀ and δ₀ 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 $\tfrac{C}{{{MR}}^{2}}$ can thus be inferred from the obliquity (e.g., Baland et al. 2016). By adjusting the coordinates of the pole α₀ and δ₀ in our POD solution (Table 1), the resulting uncertainty of the pole obliquity is 034 and 060 for the case with radio and optical data and radio only, respectively. Depending on the obliquity central value, e.g., = 00015 from Chen & Nimmo (2011) and = 0000356 from Baland et al. (2016), 1σ of C/MR² 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 (ϕ). 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₂, 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-standard-deviation formal uncertainty of k₂ equal to ∼1.2 × 10⁻⁴. By using our estimate of Enceladus's MoI (0.335, Figure 4) and our current knowledge of its mass (1.0802 × 10²⁰ kg; e.g., Nimmo et al. 2011), we ran Monte Carlo simulations to predict the Love number k₂ 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 χ² 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₂(1σ ∼ 1.7 × 10⁻⁴), which is used to better understand the mechanisms of internal power dissipation.

Zoom In Zoom Out Reset image size

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 (Tajeddine 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.

Zoom In Zoom Out Reset image size

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 degree-strength 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 degree-strength 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 $\tfrac{40\ \,\times \,{10}^{-5}}{{l}^{2}}$. 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).

Zoom In Zoom Out Reset image size

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₂). 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₂ 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)

$$\begin{eqnarray}&&{Z}_{T}(l)=4\pi G\left(\frac{l+1}{2l+1}\right){\rho }_{\mathrm{ice}}\left[1-{\left(1-\frac{{T}_{c}}{{R}_{{\rm{E}}}}\right)}^{l+2}\left(\frac{{g}_{t}}{{g}_{b}}\right){C}_{l}^{0}\ \right],\end{eqnarray} \tag{ 9 }$$

where ρ_ice = 925 kg m⁻³ 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) and (15) in Hemingway & Mittal 2019), and ${g}_{t}=\ {GM}/{R}_{{\rm{E}}}^{2}$ and ${g}_{b}=\ \tfrac{{GM}}{{\left({r}_{\mathrm{core}}+{d}_{\mathrm{ocean}}\right)}^{2}}\,-\tfrac{4}{3}\ \pi G{\rho }_{\mathrm{Shell}}\tfrac{{R}_{{\rm{E}}}^{3}-{\left({r}_{\mathrm{core}}+{d}_{\mathrm{ocean}}\right)}^{3}}{{\left({r}_{\mathrm{core}}+{d}_{\mathrm{ocean}}\right)}^{2}}$ are the gravity accelerations at the top and bottom of the shell, respectively. By using the nonhydrostatic topography coefficient ${H}_{30}^{\mathrm{nh}}$ computed from Enceladus's measured topography and the model-based admittance (Equation (9)), we derive the zonal harmonic coefficient in our MCMC algorithm,

$$\begin{eqnarray}&&{J}_{3}=\frac{{Z}_{T}(3)}{3+1}\frac{{R}_{{\rm{E}}}^{2}}{{GM}}{H}_{30}^{\mathrm{nh}}.\end{eqnarray} \tag{ 10 }$$

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. 2017), 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 $\left({C}_{l}^{0}\ne 1\right)$ 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) in the 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 plug-in 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¹²–10¹⁶ Pa s. The fluid behavior of the ocean is approximated by treating this layer as a low-viscosity material (10³ Pa s) with negligible rigidity. This approximation is valid as long as the Maxwell time of the layer (τ_M = η₀/μ₀) is far from the forcing period (Roberts & Nimmo 2008). The core of Enceladus is modeled as a viscoelastic body with shear modulus ${\mu }_{\mathrm{Core}}$ allowed to range between 50 and 70 GPa and viscosity ${\nu }_{\mathrm{Core}}$ in the range 10¹⁶–10²² Pa s.

The Bayesian inversion is thus carried out by using the following observations: M = 1.080 22 ± 0.00108 × 10²⁰ kg, MoI = 0.335 ± 0.00016, J₃ = −118.0 ± 0.24261 × 10⁻⁶, and Love number Re(k₂) = 0.02 ± 3.5 × 10⁻⁴ and Im(k₂) = 0.01 ± 1.7 × 10⁻⁴. 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₃and k₂ 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, 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₂, and J₃ suggest that our gravity investigation would significantly enhance our understanding of the core properties (i.e., core radius ${r}_{\mathrm{core}}={196}_{-2}^{+2}\$ km and density ${\rho }_{\mathrm{core}}={2336}_{-12}^{+16}\$ kg m⁻³ compared to Cassini observations. The joint analysis of k₂ and J₃, 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}_{\mathrm{ocean}}={37}_{-2}^{+2}\$ km, and the ocean density is estimated to be ${\rho }_{\mathrm{ocean}}={1119}_{-43}^{+16}$ kg m⁻³, which shows a level of uncertainty of ∼30 kg m⁻³.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The imaginary part of the Love number k₂ $\mathrm{Im}\left({k}_{2}\right)$ quantifies the lag of the tidal response with respect to the forcing, and it is related to the total energy dissipation rate $\dot{E}$ within the body (Segatz et al. 1988)

$$\begin{eqnarray}&&\dot{E}=-\tfrac{21}{2}\mathrm{Im}\left({k}_{2}\right)\tfrac{{\left({{nR}}_{{\rm{E}}}\right)}^{5}}{G}{e}^{2},\end{eqnarray} \tag{ 11 }$$

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₂ 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 $\mathrm{Re}\left({k}_{2}\right)$ 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 ${\mathrm{log}}_{10}({\nu }_{\mathrm{ice}})$ 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, ${\mu }_{\mathrm{core}})$, 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.

Zoom In Zoom Out Reset image size