OPS4
Session assets
Introduction and Abstract
The RADAR onboard Cassini probe (2004-2017), operating in active and passive (or radiometry) mode at 2.2-cm wavelength recorded the backscatter from the surface (through normalized radar backscatter cross section σ0 ) in its active mode, while it measured brightness temperatures (Tb) in its passive mode. Both resolved and unresolved observations of Enceladus have concluded on the extremely radar-brightness of Enceladus, the largest in the Solar system [1,2]. Such radar-brightness can be partially explained by the presence of ultra-clean water ice particles at Enceladus’s surface, in particular in the SPT (South Polar Terrain), which offers favorable medium for scattering. Nevertheless, so far none of the purely random wave scattering models had succeeded to reproduce the measured σ0 or radar albedos. Furthermore, during Cassini’s unique close flyby of Enceladus’ SPT (E16 flyby) swathing an area of few tens of kilometers, Cassini RADAR used as a radiometer revealed thermal anomalies that had not been detected in the infrared [3]. However, the amplitude of the internal heat flux remained to be constrained.
To understand better both the scattering and thermal anomalies of Enceladus’ SPT, we developed a model able to jointly simulate E16 Cassini active and passive observations. For the first time, our concurrent active and passive radar simulations considering a rough, undulating surface reproduces the observed radar brightness and thermal emission measurements at microwave wavelengths. This unified approach put constraints on the SPT’s subsurface properties, refining estimates of ice grain size, subsurface porosity, and endogenic heat flux.
Method
To predict backscatter (σ0) and thermal emission in the microwave domain we combine two models: (1) a thermal model providing depth profiles of the physical temperature beneath the surface at E16 flyby epoch (2) a radiative transfer model to simulate both active and passive observations.
Thermal Model
We adapted a multi-layer thermal model called MultIHeaTS [4] to the case of Enceladus. We account for Solar flux and radiative flux equilibrium at surface and a constant temperature at the bottom. The subsurface of Enceladus is modeled as bi-layer medium with an icy porous regolith overlying a denser water ice substrate. Main parameters of model are porosities of the top and bottom layers ϕ1 and ϕ2 (primarily control effective thermal properties), thickness d of top layer and the ocean level i.e ice shell thickness, D at the SPT region. Based on plume deposition rate modeling [5] a value of up to 90% was assumed for ϕ1. The thickness d is unknown, but [5] suggest it could be up to 700 meters and at least few decimeters based on [6]. We thus vary d from 1 m to 500 meters. Lastly, ice shell thickness, D is assumed to be in the range 2 to 5 km as the average ice shell thickness in the SPT region falls within this range [10]. At these depths, we assume the presence of an ocean, where the temperature corresponds to the melting point of water ice. Once thermal equilibrium is reached, the heat flux driven by temperature gradient is calculated providing the endogenic heat flux for the given subsurface properties.
Fig. 1a displays simulated temperature profiles obtained for D = 2.0 km at the SPT for different values top layer thickness. Temperature profiles clearly show a discontinuity at transition from top to the bottom layer, top layer acting like an insulating layer (especially if very porous). Fig. 1b shows the endogenic heat flux for same porosities, top layer thickness, d = 10 meters and for different ocean levels.
Radiative Transfer Model
We use Snow Microwave Radiative Transfer model (SMRT), a multi-layer RT model initially designed for snow or sea-ice [8]. Permittivity of water ice is assumed constant with temperature [11], effective permittivity depends on its porosity and includes possible contaminants fraction, assumed as organic dust. The parameters of RT model are thus dust fraction and water ice grain radius size (r1, r2) in the top and bottom layers. It is consistent with the assumptions on ϕ1, ϕ2 and d considered in thermal model. We vary r from 100 microns to 1 mm [9]. We consider two types of surface- smooth where reflections and transmissions are specular in nature given by Fresnel coefficients and rough, gently undulating surface parameterized statistically by root mean squared slope (s), which applies geometrical optics solution with Kirchhoff’s approximation [12].
Results, conclusion and discussion
Figure 2 compares modeled radar backscatter and brightness temperatures for a bi-layer medium at one E16 swath location, using the thermal profile from our thermal model. Results shows a clear anti-correlation between σ0 and Tb. Smooth surface assumption fail to reproduce the observed high σ0, while a gently undulating rough surface of moderate RMS slope reproduces both high σ0 and Tb for a highly porous (more than 80%) layer of clean water ice having a thickness of around 50 meters and containing large ice grains (500–600 µm). This provides an insight to an outstanding puzzle, as no prior random scattering model could explain Enceladus’ extreme radar brightness. Rough surface causes diffuse reflection and transmission, consequently volume scattering from the diffused transmitted waves in the porous icy regolith favors the radar back-scatter. Surprisingly for these high σ0 cases, electrical skin depth is much shallower (around 3 m) than the expected for a homogeneous pure water ice substrate. Moreover, it is consistent with high similarity between the optical and SAR images of the SPT region observed during E16.
Figure 3 demonstrates that σ0 remains stable and as high as 2 across incidence angles (up to 50◦) for highly porous medium, which may explain Enceladus’ high disk-integrated radar albedo [2]. Finally, from the best deduced parameters (ϕ1, ϕ2, d and D), thermal model constrains endogenic heat flux to 4–6 mW/m2 (Figure 1b).
How to cite: Raza, M. S., Le Gall, A., Schmidt, F., Bonnefoy, L., Picard, G., Mergny, C., and Leyrat, C.: Enceladus’ SPT surface and subsurface properties as constrained by Cassini RADAR observations, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-797, https://doi.org/10.5194/epsc-dps2025-797, 2025.
Almost all the known ocean-bearing moons of our solar system orbit within the magnetosphere of their host planet. Given these oceans are expected to contain salt, any fluid flow therein will constitute a moving conductor within the ambient magnetic field, permitting ‘motional induction’, a mechanism where the deflection of moving ions and resulting electrical shorting gives rise to secondary magnetic fields. At Earth, the motionally-induced magnetic signature of the ocean tides has been successfully extracted from magnetometer data obtained at satellite altitude. The signature of Earth’s time-mean ocean circulation is yet to be detected at altitude but has been modelled numerically and extensively. In contrast, motional induction at icy moons has received little attention and, in the case of Enceladus, has not been explored.
At Enceladus, ocean salinity has been demonstrated as a fundamental control upon the stratification and time-mean ocean circulation, uncertainty in which leading to orders of magnitude of uncertainty in the bottom-to-top transit timescale of particulates from the ocean depths to the south polar plumes (Ames et al., 2025). Considering this, we set out to explore whether motional induction within Enceladus could be large enough to be detectable and, if so, whether it could be used to constrain the ocean salinity and, by extension, timescales of transport within the ocean.
To this end, we conduct fully global, 3D simulations of Enceladus’ ocean circulation using a general circulation model (MITgcm) at varying ocean salinity. The model is non-hydrostatic, and accounts for geothermal heating, non-linearity in the equation of state for water density, as well as variation in the pressure-dependant freezing temperature at the ocean top. Simulated time-mean ocean flow and conductivity fields are then inserted into GEMMIE (https://gitlab.com/m.kruglyakov/gemmie), a global 3D electromagnetic solver, to obtain the ocean induced magnetic field (OIMF) at Enceladus’ ice surface.
Our results show that different salinity oceans generate different magnetic signatures. For example, the polarity of the dipole in the radial OIMF component reverses at low vs high salinity (see figure 1). The OIMF is dominated by the meridional (north-south) and radial ocean flows, with weak contributions from zonal (east-west) flows owing to the near perfect alignment of Saturn’s rotation and magnetic dipole axis. We find 3D variability of the flow field to be important. Longitudinal variation in the ocean flows permit higher order components of the OIMF that can significantly raise the magnitude of the induced signature. However, we find that the magnitude of Enceladus’ OIMF is very likely beyond the detection capability of modern fluxgate magnetometers located on its ice surface, especially where the assumed ocean salinity is low. Reasons for a weaker OIMF at Enceladus vs Earth include a much weaker (~ two orders of magnitude) ambient magnetic field, a more sluggish ocean circulation and a weaker ocean conductivity (owing to lower temperatures). Our results provide a first estimate of the ocean induced magnetic field at Enceladus and inform future developments in the field.
Fig.1 Left: Simulated colatitudinal ocean flow velocity (m/s), here shown at the ocean top. Centre: Longitudinal component of the electrical current (amperes; also a function of the radial ocean flows) at the ocean top. Right: Radial component of the ocean induced magnetic field (nT; OIMF) at Enceladus’ ice surface. Solutions are shown at 5 (top) and 35 (bottom) g/kg ocean mean salinity. Note colour bars are saturated throughout.
References:
Ames, F., Ferreira, D., Czaja, A., Masters, A., 2025. Ocean stratification impedes particulate transport to the plumes of Enceladus. Communications: Earth and Environment. DOI 10.1038/s43247-025-02036-3
How to cite: Ames, F., Ferreira, D., Czaja, A., and Masters, A.: Simulating the motional induction of 3D time-mean ocean flows within Enceladus, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-1885, https://doi.org/10.5194/epsc-dps2025-1885, 2025.
At first glance, Saturn’s moon Mimas might seem like an unlikely candidate ocean world. Its small size, low bulk density, heavily cratered surface, and fairly eccentric orbit all suggest a solid interior primarily composed of H2O ice. This view was first challenged by measurements of Mimas’ physical librations from Cassini ISS data [1], whose amplitude exceeded what is expected from a completely frozen, hydrostatic interior—instead hinting at the presence of either an elongated core or a subsurface ocean. Recent work has further suggested that Mimas’ periapsis drift and current high eccentricity can be explained by a young, non-steady state ocean. This paints a picture in which Mimas’ ice shell has a thickness of 20–30 km and is currently thinning, resulting in an ocean that is only 10–25 Myrs old [2][3]. While this scenario is mathematically viable, more work is needed to reconcile it with surface observations, particularly the absence of compressional features [4] and the apparent lack of relaxation of large crater topography, both of which would typically be expected if an ocean were present.
Topographic relaxation, the process by which surface relief is lost as topographic stresses decay, has been widely used to probe the thermal history of icy moons [e.g., 5]. The rate of relaxation depends on the composition and thermal structure of the moon’s interior, as well as the size of the topographic feature in question. The existence of a subsurface ocean implies a warm, weak interior that would promote relaxation, seemingly contradicting the minimally relaxed state (< 10% [6]) of Mimas’ largest crater, Herschel. Given that bigger topography generates larger stresses and thus relaxes more rapidly, Herschel crater’s preserved relief offers a unique opportunity to test for the existence and longevity of the hypothesized ocean. Determining the conditions under which Herschel’s relief could be maintained can constrain whether its present-day topography is compatible with a young subsurface ocean, and if so, how long this ocean may have existed for.
In this work, we use the finite element software COMSOL Multiphysics to forward model the relaxation of Herschel crater over 25 Myr, the maximum proposed lifespan for the young ocean. We test 3 different ice shell temperature structures, varying our model’s effective surface temperature from 80 K (Mimas’ nominal surface temperature) to 100 and 120 K, corresponding to 35–50% porosity in the uppermost ~1 km of the ice shell. Our domain consists of a 2D axisymmetric shell of radius 198.2 km, with Herschel’s crater topography at the top (Figure 1A) and gravitational acceleration decreasing radially inwards from 0.064 ms-2 at the surface. Our initial crater has a diameter of 140 km, depth of 11 km, and a 6 km tall central peak, based on measurements by [7]. We test ice shell thicknesses ranging from 15–35 km, corresponding to surface heat fluxes between ~18–50 mWm-2. We assume a pure water ice, viscoelastic shell with no impurities, where viscous behavior is driven by dislocation creep and grain boundary sliding [8]. In terms of the temperature profile, our models assume the shell to be purely conductive, with a surface temperature between 80–120 K and a base temperature of 273 K, at which point the remainder of the domain becomes isothermal (Figure 1B). We use a temperature-dependent thermal conductivity for ice of 651/T.
We find that in models with a surface temperature of 80 K, Herschel retains most of its original relief after 25 Myr across all the ice shell thicknesses tested (15–35 km). During this time, the rim-to-floor depth decreases only by 210-947 m, corresponding to relaxation fractions between ~2–9% (Figure 2). These results suggest that in the absence of porosity, minimal relaxation of Herschel crater is consistent with the existence of a subsurface ocean. However, the addition of a near-surface porous layer— whose insulating effect is simulated by increasing the effective surface temperature of our models— significantly speeds up the relaxation process. For an effective surface temperature of 100 K (corresponding to a moderately porous, 1-km thick insulating layer) Herschel relaxes by >15% for the thinnest ice shells tested (Figure 2), far exceeding the currently observed relaxation fraction. Under these conditions, maintaining Herschel’s present-day relief over 25 Myr requires an ice shell thicker than 25 km. With an even higher near-surface porosity ( = 120 K), Herschel shallows by over 10% in merely 1.2 Myr to 2.4 kyr, strongly constraining the longevity of the proposed ocean.
Our results suggest that the relaxation state of Herschel could be consistent with the existence of a young ocean on Mimas only if its ice shell has negligible to no porosity. Our models highlight the significance of near-surface temperature variations and suggest that moderate-to-high porosity would be difficult to reconcile with Herschel’s observed morphology if the proposed ocean does indeed exist. However, we note that our results alone cannot confirm or reject the existence of an ocean, but rather place possible constraints on its longevity. Ongoing work aims to incorporate dynamic processes in the ice shell, particularly the thinning proposed by [3]. Ultimately, this will refine the range of ice shell thicknesses and structures that could be consistent with Herschel’s minimally relaxed state.
Figure 1. Model set up for an example 30 km thick ice shell. A) Close-up of the geometry defined for Herschel, B) Ice sheet temperature gradient.
Figure 2. Modeled relaxation fraction of Herschel after 25 Myr as a function of ice shell thickness for a surface temperature of 80 and 100 K.
References:
[1] Tajeddine et al., Science, 346, 6207, 322–324 (2014). [2] Lainey et al., Nature, 626(7998), 280–282 (2024). [3] Rhoden et al., Earth and Planetary Science Letters, 635, 118689 (2024). [4] McKinnon et al., 56th Lunar Planet. Sci. Conf., Abstract #2897 (2025) [5] Bland et al., Geophysical Research Letters, 39(17) (2012) [6] Schenk et al., 56th Lunar Planet. Sci. Conf., Abstract #2435 (2025) [7] Moore et al., Icarus, 171, 2, 421–443 (2004) [8] Goldsby and Kohlstedt, JGR: Solid Earth, 106, 11017–11030 (2001).
How to cite: Blanco-Rojas, M. and Sori, M.: Does Mimas wear an ocean on its sleeve? Testing for a young ocean with crater topography and ice shell porosity, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-433, https://doi.org/10.5194/epsc-dps2025-433, 2025.
Dione, the fourth largest moon of Saturn, is primarily composed of water ice [1] with a surface that shows abundant evidence of tectonic resurfacing, both contractional and extensional [e.g., 2–4]. There appear to be at least two eras of tectonism on Dione: the more recent formation of the Wispy Terrains on the trailing hemisphere—some of the most recent tectonic landforms—and the leading hemisphere, dominated by more ancient tectonic features, including ridges (Figure 1) [5]. The pattern of fractures preserved in icy planetary shells forms in response to stresses, and those patterns can be diagnostic of the stress fields in which they formed [e.g., 3,6]. The long geological history of tectonism preserved on Dione can thus provide a temporal record of the stress history, which can be linked to changes in orbital dynamics and interior environments (e.g., formation or freezing of subsurface oceans) [7].
Figure 1. Global distribution of ancient (above) and recent (below) tectonic structures across Dione. Scarps are interpreted as normal fault scarps.
Here, we use SatStress [8,9] to create model stress fields to compare with Dione’s observed fracture pattern. We divide each mapped feature into 1-km segments and compare each segment to the predicted orientation from our models at that location on Dione’s surface. Preliminary results show that both diurnal stresses and stresses due to non-synchronous rotation (NSR) can replicate the fracture pattern of more recent tectonism observed in the Wispy Terrain. Predicted orientations of diurnal stresses match 63% and 82% of our observed orientations within 15° and 30°, respectively. Similarly, the predicted orientations of stresses due to NSR match 68% and 80% of observed orientations within 15° and 30°, respectively. In both cases, our modeled global stress mechanisms reproduce the observed orientations more accurately than an isotropic population with random orientations like what would be expected from ice shell thickening alone (Figure 2). While the predicted fractures from diurnal tides and NSR match the observations at a similar level, it is important to note that diurnal stresses are on the order of tens of kPa while NSR produces much higher stresses on the order of MPa.
Figure 2. Preliminary SatStressGUI modeling results for late tectonism. The angle between predicted and modeled fracture orientations, or the orientation misfit, for each 1-km segment of a mapped graben, scarp, or trough is shown for different stress models.
Future work will include the analysis of ancient tectonic features and including additional stress mechanisms such as orbital recession and despinning. This will enable us to determine if there has been a shift in the dominant stress mechanisms throughout Dione’s geologic history. Changes in stress magnitudes and mechanisms may indicate alterations in planetary interiors, such as variations in shell thickness or ocean depth, or modifications to orbital parameters like eccentricity or obliquity.
References: [1] Thomas, P. C. (2010) Icarus, 208, 395–401. [2] Kirchoff, M. R. and Schenk, P. (2015) Icarus, 256, 78–89. [3] Collins, G. C. et al. Planetary Tectonics, Cambridge University Press (2009). [4] Collins, G. C. (2010) AGU 2010, P24A-08. [5] Martin, E. S. et al. (2024) AGU 2024, P31A-09. [6] Kattenhorn, S. A. and Hurford, T. Europa, (2009), p. 199. [7] Martin, E. S. et al. (2015) 46th LPSC, Abstract No. 1620. [8] Wahr, J. et al. (2009) Icarus, 200, 188–206. [9] Patthoff, D. A. et al. (2019) Icarus, 321, 445–457.
How to cite: Nichols-Fleming, F., Martin, E. S., Patthoff, D. A., and Watters, T. R.: Constraining the Evolution of Dione from Changes to the Dominant Stress Mechanisms and Tectonic Record, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-879, https://doi.org/10.5194/epsc-dps2025-879, 2025.
Introduction:
Digital terrain models (DTMs) form the basis of many science investigations, including geologic mapping, thermal modeling, and photometric modeling. Previous DTMs were constructed while the Cassini mission was ongoing, and contain limited images. We present global and regional DTMs (i.e., shape models) constructed using the entire Cassini spacecraft imaging campaign and the Stereophotoclinometry (SPC) software suite (Gaskell 2008, 2023). The efficacy of SPC was demonstrated during the OSIRIS-REx mission, where global and regional DTMs successfully guided the spacecraft to the surface of asteroid Bennu 77 cm from the center of the target landing site (Lauretta et al., 2021; Olds et al., 2022). These models are for the Saturnian moons Dione, Mimas, Phoebe, Rhea, and Tethys. Enceladus and Iapetus are forthcoming. These models are freely available on the NASA Planetary Data System (PDS), or will be upon completion.
In addition to the topography and albedo, we also provide detailed photometric data for each vertex of the regional models. For each image, we can provide the reflectance value, phase angle, emission angle, and incidence angle. This photometric data can be used to determine photometric parameters such as geometric albedo, single scattering albedo, and surface roughness for each DTM vertex.
Methods:
SPC uses two proven techniques, stereophotogrammetry and photoclinometry, and simultaneously solves for a single product. The stereophotogrammetry (stereo) solution is accomplished by triangulation of the spacecraft images to a single point on the surface. The multi-image photoclinometric solution uses each pixel from an image along with the image’s incidence and emission angles to calculate an x-slope, y-slope, and albedo for each DTM vertex. Every image contributes to the determination of the slopes and albedo value for each point on the surface. Stereo gives sparse absolute heights, while photoclinometry gives slopes that are then integrated to determine heights in between the stereo positions. The result can be a pixel-to-pixel DTM for the surface, with accurate data for both height and albedo.
Products:
The resulting DTMs can be used to generate views that are nearly indistinguishable from spacecraft images (Fig. 1).
Figure 1. Left shows a regional SPC DTM of Dione centered at -0.3 Lat 293.9 ELon rendered with the same viewing and lighting conditions as the spacecraft image on the right. Spacecraft image is a cropped view of Cassini image N1556123061 (726 m/px resolution).
To aid data users, we provide assessment data such as best image resolution and number of images (Fig. 2 and 3). The data in Figs. 2 and 3 indicate where the topography is over- or under-sampled relative to the images, and can be used to determine topographic quality. Additional assessments, such as internal sigma values (Fig. 4), give another check on the quality of the topography. Internal sigma values worse than the best image resolution (a rare occurrence), indicate the sigma values more accurately describe the topographic quality.
Figure 2. Assessment data for regional DTM in Fig. 1. Shows the best image resolution for each vertex of the DTM.
Figure 3. Assessment data for regional DTM in Fig. 1. Shows the number of images used by SPC for each vertex of the DTM.
Figure 4. Assessment data for regional DTM in Fig. 1. Shows the SPC internal sigmas for each vertex of the DTM.
As part of SPC processing, each image will be aligned to the topography within 1 pixel. With such an accurate alignment and topography, the photometric data can be calculated on a pixel-to-pixel scale. This data is given in the United State Geological Survey (USGS) Integrate Software for Imagers and Spectrometers (ISIS) cube format. An example of the local emission angle is given in Fig. 5.
Figure 5. Graphical representation of the local emission angle of image N1556123061 for each vertex of the DTM. Black corresponds to 5 deg while white corresponds to 72 deg.
Conclusion:
An SPC global DTM, along with assessment data, of Phoebe (0.3 km/vertex) is currently available from the PDS (Weirich, Gaskell, et al., 2023). Global DTMs and assessment data of Dione (1.6 km/vertex), Mimas (0.6 km/vertex), Rhea (2.2 km/vertex), and Tethys (1.5 km/vertex) are under review at the PDS and will likely be available at the time of this conference. Regional DTMs, assessment data, and photometric data cubes of Dione (~0.4 km/vertex) and Tethys (0.3 to 0.4 km/vertex) are also under review. Similar products for Enceladus and Iapetus are forthcoming.
Though not part of the Saturn system, many other SPC products like these are available, including those for asteroid Bennu, Earth’s Moon (Weirich, Palmer, et al., 2023), Mercury (Weirich, Domingue, et al., 2024), and many other objects. These products are key in the assessment of surface morphology and regolith structure, each providing insight into surface processing of these icy worlds.
References:
Gaskell, R.W., Barnouin-Jha, O.S., Scheeres, et al., (2008). Characterizing and Navigating Small Bodies with Imaging Data. Meteoritics and Planetary Science 43 (September): 1049–61.
Gaskell, R.W., Barnouin, O.S., Daly, M.G., et al., (2023). Stereophotoclinometry on the OSIRIS-REx Mission: Mathematics and Methods. Planetary Science Journal 4: id 63, 15pp.
Lauretta, D. S., Enos, H. L., Polit, A. T., et al., 2021, in Sample Return Missions, ed. A. Longobardo (Amsterdam: Elsevier), 163, doi:10.1016/B978-0-12-818330-4.00008-2
Olds R.D., Miller C.J., Norman, C.D. et al., (2022). The Use of Digital Terrain Models for Natural Feature Tracking at Bennu. Planetary Science Journal 3: id 100, 11pp.
Weirich J.R., Gaskell R., et al., (2023). Phoebe SPC Shape Model and Assessment Products Bundle V1.0. urn:nasa:pds:satellite-phoebe.cassini.shape-models-maps::1.0. NASA Planetary Data System, doi:10.26033/3k3c-5713.
Weirich J.R., Palmer E.E., et al., (2023). Lunar Topography and Photometric Cube Data V1.0. urn:nasa:pds:lunar_lro_lroc_topography_domingue_2023::1.0. NASA Planetary Data System, doi:10.17189/jdpv-z697.
Weirich J.R., Domingue D.L., et al., (2024). Nathair Facula Topography and Photometric Cube Data V1.0. urn:nasa:pds:nathairfacula_messenger_mdis_topography_domingue_2022::1.0. NASA Planetary Data System, doi:10.17189/31AX-FT48.
How to cite: Weirich, J., Domingue, D., Palmer, E., and Gaskell, R.: Digital Terrain Models and Photometric Data of Five Large Saturnian Moons, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-115, https://doi.org/10.5194/epsc-dps2025-115, 2025.
Indubitably, the most dramatic albedo hemispheric dichotomy in the Solar system occurs at Iapetus. Its leading (L) side is covered by an optically low-albedo material, contrasting with its bright trailing (T) side. This dichotomy is also visible at microwave wavelengths, which sense not only the surface but also the subsurface and can thus bring insights into the vertical extent of the dark L side deposit layer and, in general, on the variations with depth of the composition and structure of Iapetus’ subsurface. We here present all available radio maps of Iapetus and the thermal emission model we have develop to interpret them. This model combines a multi-layer thermal and radiative transfer model and its comparison to VLA and Cassini RADAR observations provide new constraints on the physico-chemical properties of both the T and L hemispheres of Iapetus.
Introduction
The origin of Iapetus two-tone coloration has long been controversial; there now seems to be a consensus toward an exogenic deposition of low-albedo materials originating from Phoebe’s vast debris ring that crosses Iapetus’ orbit [e.g., 1, 2, 3]. However, questions remain about the composition and vertical extent of Iapetus’ L hemisphere dark deposit. Measuring the thermal emission from Iapetus at different wavelengths, thus probing different depths, can help answering these questions. In particular, probing the L side of Iapetus at multiple wavelengths in the microwave domain can put a firmer number on the thickness of the dark layer thus bringing key constraints for dynamical models aiming at reproducing dust distribution on Iapetus [e.g., 4].
Iapetus radio maps
Since 2018, the Very Large Array telescope (VLA) has been used at several instances to map the thermal emission from Iapetus at wavelengths ranging from 0.9 to 6.0 cm i.e., from the Ka-band to C-band. A total of ~30 hours of observations was collected resulting in 13 radio maps (see 6 of them in Fig. 1). VLA observations only probe Iapetus dayside but with different portions of the L and T sides in the visible disk. To complete the analysis, we also consider 6 data segments collected by the Cassini RADAR acting as a radiometer during 2 targeted flybys of Iapetus (Fig. 1). These data were acquired at a wavelength of 2.2 cm, sometimes at night, and sample either the L or T hemisphere.
Both VLA and Cassini spatially-resolved maps show a clear asymmetry in the emitted flux with more flux coming from the L face of Iapetus than from the T face (e.g., May 28, 2018). Yet, the radio maps do no perfectly mimic the optical maps thus bringing complementary information
Fig. 1: Radio maps of Iapetus collected by VLA and the Cassini RADAR (Cassini deconvolved maps are from [5]). For each map, the size of the primary beam and the albedo map of the visible disk are indicated.
Thermal emission model
The interpretation of Iapetus’ radio maps requires a thermal emission model combining a thermal model and a transfer radiative model, to produce brightness temperature (in K) maps of Iapetus that can be converted into radiance (erg.cm−2.s−1/sr−1).
Model parameters: In bright terrains, Iapetus subsurface is modeled as a homogeneous column of porous water ice. On the L side, we assume a two-layer subsurface with a porous dark layer overlying the ice substrate whose characteristics are the same as those of the bright terrains. More specifically, the model parameters are: the porosity of the dark layer and of the icy substrate , the thickness of the dark layer (which should not exceed a few meters [6]), the composition of the dark material (organics, silicates or a mixture) and the degree of volume scattering in the icy substrate .
Thermal model: We adapted the “Multi-layered Implicit Heat Transfer Solver” (MultIHeaTS) [7] to the case of Iapetus. The thermal properties of the subsurface are primarily controlled by and Because the sensed subsurface temperatures are a combination of diurnal and seasonal components, the thermal model is set to provide temperature profiles down to tens of meter depths.
RT model: We used the model proposed by [8] for RT in stratified agricultural soils adding a contribution from volume scattering () in the water ice substrate. Composition, and dictate the vertical variations of the subsurface electrical properties. Silicates are more radio absorbing /emitting than organics.
Results and discussion
Each radio observation is compared to model predictions looking for the parameter combination that best reproduces it (e.g., Fig. 2). One main result is that all observations including the T side point to the need for a significant degree of volume scattering in the water ice substrate (=20-40%). In addition, in most cases, assuming a silicate dark composition provides better fits than organic one. Yet, best-fit models for observations on the L dark side generally underestimate and/or mislocate the peak emission which suggests variations in the thickness of the dust layer (Fig. 3). This will be further investigated. Future work will also focus on adjusting all data together to the model to find the best combination and in particular put a firm constrain on the thickness of the dark layer. Lastly, we will test the hypothesis that some detected cold spots could be associated with geological features such as impact craters.
Fig. 2: VLA map of May 28 2018 (left), best-fit model (=1.55) (middle) and residual map. In the model, a uniform value is assumed on the L side.
Fig. 3: Visible disk in May 28, 2018 (left) with best-fit model now assuming a thicker dust layer in the centre of the L side (right).
References:
1. Buratti et al., 2002, Icarus 155, 375–381
2. Verbiscer et al., 2009, Nature 461, 1098.
3. Dalle Ore et al., 2012, Icarus 221 (2), 735–743
4. Tamayo et al. 2011, Icarus 215, 260
5. Bonnefoy, 2020, PhD Thesis, Paris Observatory
6. Ostro et al., 2006, Icarus 183 (2), 479–490.
7. Mergny and Schmidt, 2024, Planet. Sci. J. 521
How to cite: Le Gall, A., Lellouch, E., Butler, B., Raza, S., Bonnefoy, L., Harrar, L., Leyrat, C., Schmidt, F., Mergny, C., Robin, C., and Meyer, F.-X.: Investigating Iapetus' dichotomy with multi-wavelength microwave observations, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-810, https://doi.org/10.5194/epsc-dps2025-810, 2025.
Saturn’s narrow, clumpy F ring is a region disturbed by chaotic orbital dynamics. We therefore model it as a stochastic process (specifically, a finite Markov chain). The ring appears dominated by dust in camera images, but the main mass of this ring resides in a core of elongated clumps called kittens, observed by ring occultations. Cassini UVIS sees such features about 1/3 of the time, the same frequency as the radio occultation detections of the F ring core; both have similar size distribution. We model the F ring core as kittens (transient aggregates with size100 m < 𝚫r < 3 km), including perturbations due to Prometheus encounters, resonance confinement, and mutual collisions. We solve for the stationary state. Without confinement, the probability of detecting the kittens is uniformly distributed. With corotation resonance confinement [1], the stationary state is sharply peaked, consistent with the longitudinal distribution of detections of the F ring by radio occultation. Considering shepherding alone, the Cassini radio observations are nonetheless better fit by the non-confinement stationary distribution, and even better by just 20% confined. We find acceptable fits for fractions up to 70% of the clumps shepherded in the 109:110 Prometheus CER. This alternative combines with the explanation of [1] to conclude that some fraction of the population, or some fraction of the time, the F ring is shepherded by Prometheus. We argue that the persistence of the F ring due to negative diffusion [2][3], where the ring is confined by Prometheus aligning particle when they are driven to collide when streamlines cross. We include the negative diffusion in the Markov chain using an Ehrenfest diffusion model [4]. A small asymmetry explains the distribution in resonant argument of the radio occultation detections. In all cases, Prometheus is the agent for confinement. The F ring is thus shepherded by a combination of a Prometheus corotation and a Lindblad resonance. When the center of mass of the material in the F ring is ever located at the Lindblad resonance with Prometheus, perturbations will drive negative diffusion to maintain that location. For our combined model: Shepherded fraction has the range fshep < 0.2; Diffusion asymmetry factor Pneg has the range 0.48 - 0.50. The negative diffusion thus can maintain the longitudinal distribution either alone, or in combination with shepherding, if the mean motion resonance coincides with the true core of the F ring
Comparing 3 models to the kittens detected in the Cassini radio occultations. Red line showing the purely negative diffusion model (no shepherding, only negative diffusion with Pneg = 0.40). The quality of the fit is similar to that of the Cuzzi et al. (2024) hypothesis (green line).
- Cuzzi, JN al. (2024). Saturn’s F ring is intermittently shepherded by Prometheus.
Science, 10 May 2024. - Lewis, MC al. (2012). Negative diffusion in planetary rings with a nearby moon. Icarus 213, 201.
- Sickafoose, AA and Lewis, MC (2024). Numerical simulations of Chariklo’s rings with a resonant perturber. Planetary Science Journal, February 2024.
- Feller, W (1968). An Introduction to Probability. Wiley.
How to cite: Esposito, L. W. and AlRebdi, A.: Saturn’s F Ring is Confined by Prometheus and Negative Diffusion, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-1052, https://doi.org/10.5194/epsc-dps2025-1052, 2025.
Saturn’s rings have been considered youthful—no more than a few × 108 years—because dark, non-icy micrometeoroids should steadily accumulate on the predominantly water-ice ring particles. Cassini measurements, however, show the ring is still remarkably clean. Hyodo, Genda & Madeira (2025) Nature Geoscience revisit this apparent contradiction with a three-stage, end-to-end numerical investigation that follows (i) hypervelocity impact vaporisation, (ii) post-impact vapour condensation, and (iii) electromagnetic removal of the impactor debris. Their work demonstrates that a pollution-resistance mechanism—rather than a recent origin—is sufficient to explain the rings’ pristine appearance.
1. Hypervelocity impact simulations
Typical interplanetary micrometeoroids 1–100 μm in size collide with ring particles at ~30 km s-1. Three-dimensional smoothed-particle-hydrodynamics (SPH) calculations using five-phase H₂O and SiO₂ M-ANEOS equations of state reveal peak temperatures ≥10,000 K and pressures ≥100 GPa throughout the projectile, causing complete vaporisation of its non-icy material. Only a micrometeoroid-sized volume of the icy target is vaporised; the rest is launched as slow, icy ejecta that remains in the ring plane. Consequently, no pristine, solid non-ice grains are emplaced into the rings at the moment of impact.
2. Vapour expansion and condensation
The mixed vapour cloud expands ballistically. Adiabatic cooling drives silicate gas to the liquid–vapour boundary, where homogeneous nucleation produces nanometre-sized condensates while ~60 % of the vapour remains atomic or molecular. In contrast, H₂O vapour almost never reaches the densities required for nucleation and therefore stays gaseous.
3. Charging and dynamical evolution
Both the residual vapour and nano-condensates become ionised in Saturn’s magnetosphere. Hybrid N-body–Lorentz integrations track test particles across the full radial extent of the C, B and A rings while varying charge-to-mass ratios q/m = 10-7–10-6 e amu-1. Lorentz coupling deflects low-mass, highly charged particles out of the Keplerian mid-plane; most are either (a) precipitated into Saturn’s atmosphere (“ring rain”) or (b) accelerated along open field lines and ejected as high-speed dust streams. The fraction that re-impacts the rings—the accretion efficiency η—is only ~1–3 %, two orders of magnitude smaller than the ≥10 % assumed in earlier age estimates.
4. Implications for apparent ring youth
If η < 1 %, the timescale required to darken Saturn’s rings by micrometeoroid flux lengthens from ~102 Myr to several Gyr, consistent with the Solar System’s age. Thus the observed whiteness is not prima facie evidence of recent formation; it can be maintained by continuous self-cleaning of exogenic silicates. The same mechanism naturally produces (i) the nanometre-sized “stream particles” discovered by Cassini and (ii) the water-rich ion rain measured in Saturn’s ionosphere, reconciling multiple Cassini data sets with a single process chain.
5. Broader consequences
Pollution resistance should operate anywhere hypervelocity impacts strike porous ice in a strong planetary magnetosphere. It may therefore influence the colour and inferred ages of Uranus’ and Neptune’s dark rings and of bright terrains on icy moons. Future work must include porosity, grain-scale granularity, mixed silicate-metal chemistries, and time-dependent stochastic charging to refine η, but the first-order conclusion remains: ring cleanliness does not mandate ring youth.
How to cite: Hyodo, R.: Micrometeoroid Pollution Resistance Sustains the Pristine—and Possibly Ancient—Rings of Saturn, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-36, https://doi.org/10.5194/epsc-dps2025-36, 2025.
Introduction and Scope of Work:
Imaging of Saturn’s rings by the Cassini spacecraft revealed fine detailed structure of its dense rings [e.g., Spitale & Porco 2010, Tiscareno et al. 2019], including the existence of i) irregular vertical structures (IVS’s) on the edge of the A- and B-rings, ii) propeller structures created by embedded moonlets, and iii) mountain-like wall structures at the edges of rings from fully gap-opening satellites (Fig. 1). We have been investigating these vertical structures in order to better understand their formation mechanisms, to provide better quantitative constraints on the physical properties of ring particles, and to construct scaling relationships of structure properties with ring particle and moonlet properties. Furthermore, these physical features of Saturn’s rings provide a natural laboratory to explore the relationship between dense rings, embedded moonlets, and exterior satellites.
Figure 1. (a) An example of an Irregular Vertical Structures (IVS) in the B-ring outer edge [Spitale & Porco 2010]. (b) An embedded moonlit forming a propeller (named Blériot) on the unlit side of the rings [Tiscareno et al. 2019]. (c) Cassini image of the Keeler gap, showing Daphnis and the vertical structure of density waves induced on the edge casting shadows on the ring.
We use the N-body code pkdgrav to directly simulate the formation of i) IVS’s and ii) propeller structures from embedded moonlets. pkdgrav has the ability to model dense granular regimes by incorporating detailed multi-contact, frictional, and cohesive forces between particles through a Soft-Sphere Discrete Element Method (SSDEM) [Schwartz et al. 2012, Zhang et al. 2018, Ballouz et al. 2017]. For modelling IVS’s, we also use the N-body code epi_int [Hahn et al. 2025] to generate initial conditions for modeling the edge of the B-ring that is influenced by the Mimas 2:1 Linblad resonance [Spitale & Proco 2010].
While there have been previous simulation studies of vertical structure using hydrodynamic codes, which capture the scale of the phenomenon, and N-body simulations, which directly model the ring particle interaction, few prior studies have been able to directly model vertical structure formation at the appropriate physical scale for direct comparison with observations. Here, we leverage the highly-parallel nature of pkdgrav to directly simulation the formation of vertical structures with physically-realistic ring particle sizes (1-10 m). For the presentation, we will present progress in modeling the formation of IVSs and propeller structures from embedded moonlets. Here, we highlight some progress in the modeling of propeller structures.
Simulation Setup and Preliminary Results
We simulate a patch of particles with a differential size frequency distribution defined by a power law with exponent -2, and a radius range of 0.5 to 1.5 meters. We vary the bulk density of individual particles, and set a nominal value of 0.6 g/cm3. The radius of the embedded moonlet is 50-m. We vary the normal, 𝜀n, and tangential, 𝜀t, coefficients of restitutions. The friction parameters are chosen such that particle shave an angle of friction of ~ 32°. The simulations are initialized by creating a uniform density ring patch with a prescribed dynamical optical depth (𝜏dyn = 1.0) set in the B-ring (semi-major axis = 100,000 km). The azimuthal and radial extents of the patch are approximately 6 km and 1 km, respectively, which is approximately 85 x 15 embedded moonlet hill radii. The total number of ring particles simulated is N = 3.6 million. The simulations run for approximately 50 orbital periods, which is sufficient to reach a steady state configuration.
Fig. 2 shows examples of two simulations with different 𝜀n and 𝜀t. In both cases, we see the formation of propeller structures that extend to 4-5 km along the azimuthal direction. For the case with more dissipative collisions, we noted the formation of irregular dense vertical structures at the edges of the gaps (magenta box in Fig. 2b) that have a saw-tooth pattern, which have been observed in different regions of Saturn’s rings. Fig. 3 shows a closer look at some of these structures.
Figure 2. End-stage of embedded moonlet simulations with patch size of 1.1 km (radial) by 6 km (azimuthal). The propeller structure is evident and contained within the confines of the simulation box. a, case with less dissipative inter-particle collisions, 𝜀n = 0.8, 𝜀t = 0.8. b, case with less dissipative inter-particle collisions, 𝜀n = 0.5, 𝜀t = 0.5.
Figure 3. We highlight one of the regions with irregular saw-tooth structures in the simulation results shown in Fig. 2b. These structures are more pronounced for the case with more dissipative collisions between ring particles.
We analyzed the case shown in Fig. 2b case further to study how the structure seen in the simulations might control the appearance of the ring as seen by Cassini. We sub-divided the ring patch into a grid made of 20 m × 20 m cells and calculate the maximum height, H, of each cell, normalized by the embedded moonlet’s radius, rm. We find that parts of the vertical structure have H/rm > 1. Fig. 4 shows that the height of the ring patch is largest adjacent to the accretion tubes.
Figure 4. Sub-dividing the simulation shown in Fig. 2b into 20 m × 20 m cells, we measured the height range, H, normalized by the embedded moonlet’s radius (rm) in each cell to better compare with Cassini observations of Saturn’s dense rings.
Acknowledgements
This work is supported by the NASA Cassini Data Analysis program through grant 80NSSC23K0220. This work also made use of a NASA High-End Computing allocation (SMD-24-30244393)
References
Ballouz et al. 2017 AJ 153, 146.
Ballouz et al. 2021 MNRAS 507, 5087.
Hahn, J.M., et al. 2025 ApJ, In Review.
Schwartz, S.R., et al. 2012 Granular Matter 14, 363
Spitale, J.N., & Porco, C.C. 2010 AJ 140, 1747.
Tiscareno, M., et al. 2019 Science 364, 6445.
Zhang, Y., et al. 2018 AJ 857, 15.
How to cite: Ballouz, R.-L., Spitale, J. N., Hahn, J. M., and Schwartz, S. R.: Investigating Vertical Structure Formation in Saturn’s Rings with N-body simulations, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-336, https://doi.org/10.5194/epsc-dps2025-336, 2025.
Stellar occultation measurements of Saturn’s rings from the Cassini Ultraviolet Imaging Spectrograph (UVIS) High Speed Photometer (HSP) reveal dozens of regions in the heart of the B ring that are opaque at the sensitivity level of the measurements (Colwell et al. 2021, 2024). We have identified 24 regions with an average transparency less than 0.25%, including 9 with an average transparency less than 0.1% and radial extents greater than 10 km. There are broad expanses within these opaque regions where the transparency is indistinguishable from zero.
Occultations of the star Beta Centauri with photon count rates of up to 600 per 1 ms integration period provide the strongest constraints on transparency. The elevation angle B (sin(B)=0.918) and brightness of Beta Centauri make it best suited among the UVIS stellar occultations to probe high optical depth regions. The radial resolution of these occultations is about 8 meters, while the azimuthal resolution in the frame co-moving with the ring particles is ~20-30 meters. Each measurement samples an area in the frame of the particles of ~400 m2. The average transparency of each designated opaque region is > 0 because of the presence of narrow, transient features we dub “phantoms” that increase the average transparency of the region. Phantoms are analogous to the ghosts observed by Baillie et al. (2013) and Green et al. (2024) in the C ring and Cassini Division. While ghosts were defined by being indistinguishable from 100% transparent in an otherwise moderate optical depth background, phantoms are defined as narrow, non-axisymmetric features with non-zero transparency in an otherwise opaque background (Figure 1).
Figure 1: Stacked and offset (by 5%) transparency profiles of the O1 opaque region in Saturn’s B2 ring region as measured by the Cassini UVIS HSP. The different curves are for different observations of the same star system: Beta Centauri. Narrow spikes in transparency are identified as phantoms. The vertical dashed lines indicate the boundaries of the opaque region.
We observed phantoms that are only a single HSP measurement, indicating a radial extent of < 10 m, but most span several measurements. Phantom transparencies range from a few percent up, but in some cases the transparency exceeds 10% and can approach 50%. There are 15 measurements of the B ring with Beta Centauri occultations, and the phantoms do not repeat from one occultation to the next, highlighting their non-axisymmetric and/or temporally variable nature. The Beta Centauri system has two stars that contribute to the HSP signal, and the projected separation of these two components in the ring plane enables us to place limits on the azimuthal extent of some phantoms, which in some cases is < 100 m.
At some locations there are clusters of phantom-like features that appear in all occultations at the same location. We identify these axisymmetric features “grassy regions” based on their appearance in plots of transparency (Figure 2). The average transparency of grassy regions is still low, typically less than 2%. The opaque regions are bordered by regions of higher optical depth (Figure 2).
Figure 2: Stacked and offset (by 5%) transparency profiles of an opaque region in Saturn’s B2 ring region as measured by the Cassini UVIS HSP. The different curves are for different observations of the same star system: Beta Centauri. This opaque region is home to phantoms (isolated spikes in transparency) as well as several grassy regions, identifiable by their presence at the same location in all occultation profiles and marked by blue dashed lines.
The measured optical depth of the opaque stretches between phantoms can exceed 6 in some cases (see for example the bottom curve in Figure 2 between 100,250-100,300 km). N-body simulations have struggled to reproduce such high optical depths. Individual ring particles in the simulations clump together due to their mutual gravity, leaving gaps which elevate the overal transparency of the region. We will present results of simulations with modified boundary conditions and particle properties to try to explain the existence of such high-optical-depth regions and the sub-structure of phantoms and grassy regions within them.
Baillié, K., J. E. Colwell, L. W. Esposito, and M. C. Lewis 2013. Meter-sized Moonlet Population in Saturn’s C Ring and Cassini Division. Astron J. 145, 171, doi:10.1088/0004-6256/145/6/171.
Colwell, J. E., M. Brooks, R. Jerousek, C. Coleman, M. S. Tiscareno, K.-M. Aye, M. Lewis, L. W. Esposito 2021. Irregular Structure in the Core of Saturn’s B Ring. 2021 AGU Fall Meeting P34A-01.
Colwell, J. E., A. Goforth, C. Coleman, R. Jerousek, L. W. Esposito 2024. Saturn’s B Ring Phantoms: Frequency, Properties, and Distribution. 2024 Meeting of the American Geophysical Union. P33F-2931.
Green, M. R., J. E. Colwell, L. W. Esposito, R. G. Jerousek 2024. Particle sizes in Saturn’s rings from UVIS stellar occultations 2. Outlier populatinos in the C ring and Cassini Division. Icarus 416, 116081.
How to cite: Colwell, J., Faulkner, E., Goforth, A., Jerousek, R., and Esposito, L.: Opaque Expanses in Saturn’s B ring, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-1101, https://doi.org/10.5194/epsc-dps2025-1101, 2025.
The presence of an ocean beneath the Enceladus’ ice shell makes this Saturnian moon a high priority target for future planetary exploration [1]. Water jets that have been observed at the south pole by NASA’s Cassini mission [2] are thought to originate from a global ocean and provide a direct window into the subsurface composition [3]. These jets generate a highly porous material that, due to its low thermal conductivity, affects the thermal state of the ice shell.
The analysis of pit chains on the surface of Enceladus indicates that locally the porous layer can be as thick as 700 m [4]. Such a thick porous layer can locally increase the temperature of the ice shell, leading to a low viscosity. This may promote solid-state convection in regions where the ice shell is covered by such a layer, whereas regions with thin porous layers could be characterized by conductive heat transport. Moreover, due to its effect on the ice shell temperature, the porous layer can strongly attenuate the signal of radar sounders that have been proposed to investigate the Enceladus’ subsurface [5, 6].
Here, we use the geodynamical code GAIA-v2 [7] to investigate the effects of a porous layer on the thermal state and dynamics of Enceladus’ ice shell. GAIA-v2 was originally built to model convection in the rocky mantle of terrestrial planets [8,9,10] but has recently been adapted to investigate large-scale dynamics in the shells of icy moons [11,12]. Our initial simulations use a 2D cylindrical geometry with an 18° arc computational domain, free-slip boundaries, and an Arrhenius rheology with diffusion creep, dislocation creep, as well as basal slip and grain-boundary sliding deformation mechanisms.
Using the resulting thermal state, we calculate the associated two-way radar attenuation at each location within the ice shell. We test different values of the ice shell thickness (5 – 35 km, [13]), porous layer thickness (d = 0 – 750 m), and its thermal conductivities (k = 0.1 – 0.001 Wm-1K-1 [14,15]). To account for chemical impurities within the ice shell we test a “low” loss scenario that considers a pure water ice shell and a “high” loss case that assumes a homogeneous mixture of water ice and chlorides in concentrations extrapolated from the particle composition of Enceladus’ plume [5].
Fig. 1: Modeled convection pattern in an ice shell with an ice shell thickness of 25 km, a porous layer thickness of 250 m, and porous layer thermal conductivity of 0.025 Wm-1K-1.
Fig. 2: Laterally averaged profiles for Fig. 1 showing a) temperature with depth b) 2-way radar attenuation under low-loss pure water ice assumption c) 2-way radar attenuation under high-loss 300 µmol chloride assumption. We assume here an attenuation limit of 100 dB [16].
Our results show that the porous layer thickness and its distribution have a first order effect on the thermal state and dynamics of the ice shell. Regions covered by a thick porous layer are characterized by a warm ice shell temperature and thus a lower viscosity, becoming more prone to undergo solid-state convection (Fig. 1). The vigor of convection depends on both the temperature-dependent ice shell viscosity and the temperature difference across the ice shell. While a thick porous layer would result in a low ice shell viscosity, thus increasing the convection vigor, such thick porous layers lead to an almost isothermal ice shell, due to their strong insulation, which, in turn, decreases the convection vigor. In ice shell regions covered by a thick porous layer, the penetration depth of a radar sounder is limited (Fig. 2), but as discussed in a recent study that only investigated a purely conductive ice shell [6], the high temperatures obtained in this case may lead to the formation of shallow brines detectable by radar measurements.
Next, we will increase the angular size of our models and allow the ice-ocean interface at the lower boundary to vary to capture spatial variations of the ice shell thickness. The upper part of the computational domain in our model will be treated as the solid ice shell, while the lower part will mimic a water ocean layer with a lower viscosity. The viscosity of the water ocean will be set to several orders of magnitude lower than that of the overlying ice shell. We will vary the viscosity contrast between the ice shell and the water layer to investigate whether ice flow can occur and possibly dampen the variations in the ice shell thickness. This will help us understand how convection can dynamically alter ice shell thickness, in contrast to previous modeling which assumes a purely conductive ice shell [6].
In future models, we will extend our models to include the effects of chemical impurities and tidal forces on the thermal state of Enceladus’ ice shell. We will track the redistribution of chemical anomalies due to solid-state convection in the ice shell. To this end, we will test different initial distributions of chemical heterogeneities and include tidal heating. Both the presence of chemical anomalies and tidal heating can affect the convection vigor and the convection pattern that characterizes large-scale dynamics in the ice shell. We will compare these more complex models with previous chemically homogeneous models in order to examine the effects of chemical anomalies and tidal heating on the thermal and dynamical state of the ice shell, and on the 2-way radar attenuation.
References:
[1] Choblet et al. (2021), Experimental Astronomy; [2] Porco et al. (2006), Science; [3] Postberg et al. (2009), Nature; [4] Martin et al. (2023), Icarus; [5] Souček et al. (2023), GRL; [6] Byrne et al. (2024), JGR:Planets; [7] Hüttig et al., (2013), PEPI; [8] Laneuville et al. (2013), JGR:Planets; [9] Tosi et al. (2015), GRL; [10] Plesa et al. (2016), JGR:Planets; [11] Rückriemen-Bez et al. (2023), Galilean Moons Workshop; [12] Plesa et al. (2023), Galilean Moons Workshop; [13] Hemingway & Mittal (2019), Icarus; [14] Seiferlin et al. (1996), PSS; [15] Ferrari et al. (2021), A&A; [16] Kalousova et al. (2017), JGR:Planets.
How to cite: DeMers, E., Byrne, W., Plesa, A.-C., Benedikter, A., and Hussmann, H.: Modeling Convection Dynamics in the Ice Shell of Enceladus, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-2027, https://doi.org/10.5194/epsc-dps2025-2027, 2025.
Introduction: One of the most compelling discoveries in planetary science is the potential existence of habitable subsurface oceans within icy moons. Among these, Saturn's moon Enceladus stands out due to its remarkable geological activity, characterized by the continuous venting of water vapor, ice particles, and organic compounds from its south polar region [1]. Cassini's comprehensive suite of measurements, including gravity data, confirmed the existence of a global subsurface ocean [2]. Accurately modeling the internal structure of Enceladus is crucial for assessing its habitability. Geophysical parameters, such as total mass and moment of inertia (MoI), are effective at constraining the total thickness of the hydrosphere, but they provide limited information on the separate contributions of the ice shell and the underlying liquid ocean. In contrast, measurements of physical librations in longitude are particularly sensitive to the internal structure of the hydrosphere, as they are influenced by the degree of mechanical decoupling provided by the liquid layer and the rigidity of the overlying ice shell. By combining mass and the moment of inertia measurements with libration amplitude estimate, a Bayesian inference approach enables tighter constraints on the deep interior and hydrosphere properties. This study presents an internal structure framework based on the Markov Chain Monte Carlo (MCMC) method that is well-suited for the estimation of Enceladus’ interior properties with rigorously quantified uncertainties by exploring the parameter space.
Methods: The mass, MoI and libration amplitude provide constraints that are integrated into a Bayesian inference framework, which is used to infer the internal properties of the investigated body. In particular, we implemented a MCMC algorithm to invert interior models by varying the free parameters associated with the structure of an icy moon. In accordance with the methodologies employed and tested in previous studies [3-4], our approach considers the body to be a multi-layered structure with free parameters of layer size, density, and rheology. These parameters are iteratively refined within the MCMC framework using the Metropolis–Hastings algorithm, which generates a diverse set of interior models. To ensure robust mapping of the parameter space, 20 independent chains are employed, each generating approximately 50,000 accepted models. Following the convergence of all chains, probability distributions for each parameter are derived, yielding constraints on the likely internal structure that are consistent with the observed geophysical data.
Enceladus Interior Model Inversion: The internal structure of Enceladus is constrained using a combination of geophysical parameters derived from Cassini mission data. The mass of Enceladus, determined from radio science data, is 1.08022 ± 0.00108 × 1020 kg [5], while the normalized MoI is estimated at 0.335 ± 0.002 [6]. These parameters place fundamental constraints on the total thickness of the hydrosphere, while remaining insensitive to the separate contributions of the ice shell and the subsurface ocean [6]. Physical librations, however, offer a crucial complementary constraint, as it is highly sensitive to the thickness and mechanical properties of the ice shell. The libration amplitude of Enceladus was determined from optical tracking of surface features by the Cassini spacecraft, revealing a physical libration of 0.120 ± 0.007°, which corresponds to an equatorial displacement of approximately 530 ± 30 m [7].
In this work, we implement the elastic libration model [8], which accounts for the gravitational torques acting on both the periodic and static tidal bulges and includes the complex gravitational interactions between the individual layers of the icy moon arising from their mutual misalignment and the pressure forces exerted by the liquid ocean on the surrounding solid layers. The model captures the contributions from both the direct external gravitational torque and the internal gravitational coupling between layers, as well as the pressure feedback from the ocean.
The moon is modeled with a three-layer structure comprising a rocky core, a subsurface ocean, and an ice shell. The thicknesses of the core and ice shell were treated as free parameters, while the ocean thickness was determined based on Enceladus’ known radius. The free parameters included the densities of the core and ocean, while the ice shell density was fixed at 917 kg m-3. The viscosities and shear moduli of the shell and core were included as free parameters. Within the MCMC framework, each proposed interior model starts from a simplified spherical approximation, characterized by an average radius and the densities of each layer. The hydrostatic shapes of the internal interfaces are then computed assuming hydrostatic equilibrium, calculating the polar and equatorial eccentricities of each layer using a fourth-order formulation [9].
The results suggest that incorporating libration amplitude can significantly improve the characterization of the internal differentiation withing Enceladus' hydrosphere, providing an estimate of the ice shell thickness with an uncertainty of approximately 1.5 km around a mean value of 21 km, consistent with previous studies.
Summary: The application of MCMC inversion with libration constraint offers a powerful approach to precisely determine the ice shell thickness of Enceladus, achieving a constraint accuracy of approximately 1.5 km. This framework can be readily adapted to other icy moons, providing a valuable tool for probing the internal structures of ocean worlds throughout the solar system.
References:
[1] Glein et al. (2015) GeCoA, 162, 202. [2] Iess et al. (2014), Sci, 344, 78. [3] Genova A. et al (2019) GRL 46(7), 3625–3633. [4] Petricca et al. (2023) GRL, 50, e2023GL104016. [5] Hemingway et al. (2018) in Enceladus and the Icy Moons of Saturn, Univ. Ariz. Press, 57. [6] Genova et al. (2024) PSJ, 5(2), 40. [7] Thomas et al. (2016) Icarus, 264: 37-47. [8] Van Hoolst et al. (2013) Icarus, 2013, 226.1: 299-315. [9] Tricarico (2014) ApJ, 782(2), 99.
How to cite: Ciambellini, M., Genova, A., Gargiulo, A. M., and Boccacci, G.: Constraining Enceladus' interior structure by using libration measurement in a Bayesian framework, EPSC-DPS Joint Meeting 2025, Helsinki, Finland, 7–13 Sep 2025, EPSC-DPS2025-1889, https://doi.org/10.5194/epsc-dps2025-1889, 2025.