Abstract

This study investigates the water plumes of Saturn’s moon, Enceladus, using Direct Simulation Monte Carlo (DSMC) modeling to analyze venting dynamics and plume structures. Building on prior research, we employ a parametrized DSMC approach to model water vapor and ice particle flows, leveraging Cassini spacecraft data from instruments such as the Ion and Neutral Mass Spectrometer and the Ultraviolet Imaging Spectrograph. The study explores whether vent conditions, such as mass flow rates, mixture temperatures, and particle sizes, can be inferred from observational data. We develop a computational framework to expand plume simulations beyond 10 km altitudes, incorporating gravitational and inertial forces in an Enceladus-fixed reference frame. A sensitivity analysis correlates vent parameters with observed data, identifying critical contributors such as vent orientation and location, mass flow rate, exit temperature, and ice grain characteristics. This approach reduces the dimensionality of fitting procedures, enabling robust parameter constraints and a more detailed understanding of plume dynamics. Key findings include constrained values for mass flow rates, ice grain radii assuming single-size particles, and exit temperatures (∼44–61 K), consistent with theoretical predictions. Additionally, variations in vent orientation and positional parameters were refined from the work of Porco et al. (2014, https://doi.org/10.1088/0004-6256/148/3/45). These results highlight the importance of collision dynamics in shaping plume structures. This work establishes a computationally efficient methodology for analyzing cryovolcanic plumes applicable to future missions exploring icy moons such as Enceladus or Europa. By prioritizing sensitive parameters, the study offers insights for optimizing observational strategies to maximize scientific yield.

Plain Language Summary

This study examines the water plumes of Saturn’s moon Enceladus using computer simulations to understand how water vapor and ice particles escape to space. Data from NASA’s Cassini spacecraft are used together with the Direct Simulation Monte Carlo modeling technique to study the structure and behavior of these plumes. The study focuses on key factors such as the speed and temperature of the escaping material as well as the size of the ice grains—assuming single-size particles—to determine how well these could be inferred from Cassini’s observations. By refining previous models, this work extends the simulations to altitudes beyond 10 km, considering the effects of gravity and motion on the plume. The findings provide new estimates for the amount of material ejected, ice grain sizes (∼1.1 μm), and vent temperatures (∼44–61 K), improving previous models. This study could help with future missions, such as those to Jupiter’s moon Europa or ESA’s L-class future mission to Enceladus, by identifying the most important parameters for studying cryovolcanic activity and optimizing spacecraft observations.

Key Points

Updated values for the total mass flow rate outgassed from the Tiger Stripes using comprehensive Direct Simulation Monte Carlo simulations

New constraints on the ice grain radii and vent exit temperatures

Method to reduce the size of an under constrained problem

1 Introduction

The Enceladus geysers, located in its South Polar Region, were first reported in 2005 by Cassini. Various instruments onboard the spacecraft (Schenk et al., 2018; Science Special Issue, 2006) performed observations of the plumes, while many studies focused on their characterization, as they offer an unprecedented opportunity to study the interior of this icy moon of Saturn. Many studies have shown that the outgassed material is mainly composed of water vapor and water-ice particles, together with other trace gases such as nitrogen-bearing and oxygen-bearing, and aromatic compounds (Hansen et al., 2011; Khawaja et al., 2019; Magee & Waite, 2017; Postberg, Khawaja, et al., 2018; Waite et al., 2006, 2017).

The geysers originate from many tens-of-kilometer-long crevasses, named the Tiger Stripes, and the plumes extend far into space: the ice grains ultimately form the Saturn E-ring. Ground temperature in the Tiger Stripe region was measured by the Composite Infrared Spectrometer (CIRS) (Howett et al., 2011; Spencer et al., 2006) and showed that the crevasses in the South Polar region were dramatically warmer than the surrounding regions (∼135 K) and the rest of the satellite surface (<80 K). Using the Imaging Science Subsystem (ISS) imaging and triangulating the visible particulate streams, Porco et al. (2014) reported the location and orientation of 98 vents at the surface.

Several studies have focused on the origin of the geyser phenomena: near-surface explosive boiling liquid water (Porco et al., 2006), clathrate decomposition (Kieffer et al., 2006), and evaporation from an underground salty ocean (Porco et al., 2014; Postberg et al., 2009, 2011; Schmidt et al., 2008). Several authors studied the flow inside the vent using a gas dynamics model (Ingersoll & Pankine, 2010; Schmidt et al., 2008; Yeoh et al., 2015; Zaharias et al., 2023). Others focused on the formation of particles, which were explained to be the results of sprays from bursting bubbles, the latter stemming from rising exsolved gases such as CO2 (Matson et al., 2018). Other studies have described how fluids move up the conduits (Ingersoll & Nakajima, 2016; Kite & Rubin, 2016; Nakajima & Ingersoll, 2016). Many studies have concentrated on constraining the outgassed mass flow rates (Burger et al., 2007; Dong et al., 2011; Hansen et al., 2006, 2008, 2011; Saur et al., 2008; Smith et al., 2010; Tenishev et al., 2010, 2014; Tian et al., 2007; Waite et al., 2006; Yeoh et al., 2015, 2017) using various physical models of the geysers. Examples include ballistic geysers from a point source based on three source parameters (the number density at the source, the flow speed, and the source temperature) and the orientation and location of the source (Dong et al., 2011), analytic solutions of the geyser particle velocity fields (Teolis et al., 2017) using the number density expression from Dong et al. (2011), and simple Direct Simulation Monte Carlo (DSMC) models (Portyankina et al., 2016, 2022; Spitale et al., 2015; Tucker et al., 2015).

The current study is an application of previous work by Mahieux et al. (2019), which presents state-of-the-art calculations and parametrizations of axisymmetric geyser flow fields using fully spatially and temporally resolved comprehensive DSMC simulations at steady-state. Mahieux et al. (2019) considered a single geyser, outgassing water vapor and single particle size pure water-ice grains, expanding into vacuum at supersonic speed in a steady state regime up to an altitude of 10 km, where the flow is essentially non-collisional. They parametrized the flow field characteristics, that is, number density, velocity, and temperature, at the 10-km upper boundary as a function of the vent conditions: mass flow rate, vent radius, water vapor exit speed, ice grain exit speed, mixture temperature, mass ratio between the two phases of water, size of the ice grains, and spreading angle. They then built functional forms to robustly and quickly reconstruct the flow field characteristics for any vent condition within a given range, defined based on physically plausible conditions.

In the present work, we aim to study whether the considered ground conditions of the vents can be inferred from the in situ Ion and Neutral Mass Spectrometer (INMS) and remote-sensing Ultraviolet Imaging Spectrograph (UVIS) observations, considering the location and orientation of the 98 vents described in Porco et al. (2014). INMS reported measurements of the local water vapor content sampled along the orbit of Cassini as it flew through the plumes. We uploaded the INMS Cassini trajectory data from the PDS online archive and used the water vapor local number density from Perry et al. (2015). UVIS performed solar (and stellar, not considered in this work) occultations of the plumes, from which the team derived integrated column densities as a function of the spacecraft’s position and a look direction (Portyankina et al., 2018). We use the aforementioned functional forms as lower boundary conditions to simulate the flow at altitudes above 10 km, assuming collisionless flow conditions and considering the main external forces in an Enceladus-fixed reference frame. We model the INMS and UVIS observations of the Enceladus plumes. We fitted the INMS and UVIS data to the synthetic model we developed using the Levenberg-Marquardt least-squares algorithm. An unsuccessful attempt to fit the data using Markov Chain Monte Carlo is discussed in Appendix E in Supporting Information S1.

In this work, we show that in addition to the mass flow rates of various vents, some information regarding the ground conditions (mixture exit temperature and the size of the ice particles) can be constrained using the UVIS and INMS data for several vents, providing insights into underground physics and chemistry. A subsidiary result of this work is a fast and reliable method for determining which of the many physical parameters of a physically realistic plume model can be constrained by limited observational data, considering only the geometry of the observation.

Section 2 describes the method we use to expand the flow into space. We present a model Monte Carlo study of the sensitivity of the measuring instrument signals to several vent parameters in Section 3 and rank these parameters by their sensitivity to the observation geometry. Section 4 discusses the results obtained using a least-squares algorithm of the most sensitive vent parameters for the E3, E5, E7, E14, E17, and E18 INMS observations and the UVIS solar occultation. Finally, we present our conclusions in Section 5.

2 Flow Field Calculation

2.1 Parametrized Radial Profiles at 10 km Altitude

We first summarize the work presented by Mahieux et al. (2019). They performed detailed state-of-the-art simulations of the entire flow field of axisymmetric plumes and derived functional forms, which provide the number density, the radial and vertical velocities, and the kinetic temperature for water vapor and single particle size pure ice grains at 10 km altitude above a vent. These 10 km flow field conditions are used in this work as the lower boundary conditions of the simulations described in Section 2.2. Their analytical expressions of the parameterized radial profiles were calculated from comprehensive and fully spatially and temporally resolved DSMC simulations obtained using the code PLANET, which assumes that (a) the mass flow is composed of a mixture of water vapor and uniform-radius, water-ice particles exiting the ground at two independent supersonic speeds for each water phase, at a given mass ratio between water vapor and water-ice grains and uniform temperature, (b) the flow does not experience any force beside Enceladus’s gravity g E ${g}_{E}$ , (c) the vent is circular, (d) the flow exits the moon’s surface with cylindrical symmetry relative to the axis normal to the vent center, (e) Enceladus has no atmosphere, (f) there are only neutral particles and the influence of electromagnetic fields is neglected, (g) as the composite flow exits through the vent, it spreads over an angle (“spreading angle”) that varies uniformly from the vertical vent centerline to a maximum value at the vent edge to represent the initial flow, and (h) a single particle size—or mono-dispersed population—is an adequate approximation of true grain size distribution.

Assumption (f) is valid for large ice particles but might not always be true for particles with radii smaller than 1 μm, where charging can significantly affect their trajectories. However, the importance of charged sub-micron particles depends on their fraction in the plume. If it is small, they may be ignored as the neutrals would dominate the flow. Small charged particles would be swept away immediately by the swiftly rotating magnetosphere, leaving the neutral plume free of their influence. Additionally, for charged particles to play a significant role in the plume, their flux must be high enough to charge sub-micron neutrals there. At Enceladus’ distance from Saturn, the flux of trapped charged particles is low. Moreover, the rings act as a sink for charged particles by sweeping out large volumes. Therefore, there is reasonable justification for considering the neutrals-only case.

As already noted in Mahieux et al. (2019), our plume-flow framework adopts, for reasons of tractability, assumption (h): a single particle size or monodisperse population of water-ice grains is used at the vent exit, that is, all particles are assigned a single representative radius. From a purely physical standpoint, this is, of course, an approximation. Direct and indirect evidence from Cassini Plasma Spectrometer (CAPS), Cosmic Dust Analyzer (CDA), Visible and Infrared Mapping Spectrometer (VIMS), and ISS investigations, as well as the condensation and wall-erosion scenarios outlined in Dong et al. (2015), show that Enceladus actually produces multiple grain populations whose characteristic sizes depend on their formative pathway: nanometer-scale condensates nucleated in the rapidly cooling vapor, sub-micron grains formed by homogeneous nucleation deeper in the conduit, and larger-than-micron fragments ablated from warm fracture walls (Ershova et al., 2024; Kempf et al., 2008; Waite et al., 2017). Collisions between these grains and the surrounding water molecules in the dense near-vent region do, in principle, affect momentum and energy exchange: the collision rate scales with the number density of grains, which itself scales inversely with the cube of grain radius at fixed mass flux, and thus modifies the early acceleration of both phases.

Yet, several pragmatic considerations argue for retaining a single effective grain size in the present study. First, the data sets we exploit, namely INMS local water-vapor number densities and UVIS column densities, are not directly sensitive to the microphysical grain-size spectrum. What these instruments do sense is the large-scale shape—and integrated density for UVIS—of the vapor plume. Numerical computations in Mahieux et al. (2019) demonstrate that once the flow has expanded beyond the first few meters where most collisions occur, the gross vapor field that determines the INMS/UVIS signal varies by only a few percent across a realistic range of grain radii. Put differently, the vapor phase dominates the observables at Cassini altitudes, and the details of the size distribution imprint only second-order perturbations.

Second, allowing each vent to possess its own grain-size distribution (modal radius, width, multiple modes) would dramatically multiply the parameter count. Such a high-dimensional search space is severely under-constrained by relatively few measurement points available per fly-by, leading to unstable fits or indeterminate solutions.

Finally, our objective is to identify which gas-phase vent conditions can be constrained solely from vapor observations. Instruments that are directly sensitive to grain size, such as CDA, are outside the scope of the current work. Incorporating those measurements in future studies would indeed warrant revisiting the full multi-population description.

In summary, while a multi-modal grain-size distribution is undoubtedly more realistic, adopting a single representative radius is a justified compromise. It preserves the dominant gas–particle coupling that shapes the macroscopic plume, keeps the inversion problem computationally manageable, and avoids introducing poorly constrained variables that the INMS and UVIS data simply cannot resolve. We therefore retain the single particle size assumption here, but we explicitly acknowledge that future analyses integrating CDA or next-generation in situ dust detectors could and should treat the full size spectrum explicitly (see Section 3 for further discussion).

The vertical domain of the DSMC computation extends from the surface up to 10 km altitude, where the gas mixture was shown to be collisionless, that is, the gradient-based Knudsen number, K n ${K}_{n}$ , is much larger than 1 for both gas and grains. This is the case for the entire range of vent parameters considered here and in Mahieux et al. (2019). Only a 1° azimuthal cylindrical slice of the flow was computed with respect to the axis normal to the surface passing through the center of the vent under the assumptions listed above.

Mahieux et al. (2019) considered the vent parameters in Table 1 as the “default” or nominal vent parameters. They were obtained from previous work conducted by our group at The University of Texas at Austin and summarized by Yeoh et al. (2017). The analytical expressions of the parametrizations are described in detail in Appendix A in Supporting Information S1. It should be noted that the computation of a single DSMC simulation from surface to 10 km altitude takes up to 48 hr using 4 CPUs (on Intel E5-2690 v3 12-core machines), while the reconstruction of the radial profile at 10 km using the parametric fits only takes a fraction of a second on the same machine—this is the key benefit of creating the parametric fit. As discussed in Mahieux et al. (2019), the significant computation time saved is obtained at the cost of reduced precision of the density and velocity profiles at 10 km altitude when compared to the DSMC calculations for random vent parameter values. The fields obtained using the above-described parametrizations generally agree with the corresponding full DSMC simulations to within 5%. In the worst case, they differ by 20%.

Table 1. List of the Confidence Intervals for the Vent Parameters, Location, and Orientation

Vent parameter	Symbol	Units	Scale	Relative/absolute value	Default value	Min value	Max value
Mass flow rate	m ˙ $\dot{m}$	[kg/s]	Linear	Absolute value	4.81	0	100
Water vapor speed	v gas , ini ${v}_{\text{gas},\text{ini}}$	[m/s]	Linear	Absolute value	902	271	1,082
Water-ice speed	v ice , ini ${v}_{\text{ice},\text{ini}}$	[m/s]	Linear	Absolute value	902	271	1,082
Vent radius	r vent ${r}_{\text{vent}}$	[m]	Linear	Absolute value	1.4	0.2	6.0
Ice grain mass ratio	f ice ${f}_{\text{ice}}$	[-]	Linear	Absolute value	5	0	0.5
Ice grain radius	r ice ${r}_{\text{ice}}$	[μm]	Log	Absolute value	1	0.01	5
Mixture temperature	T ini ${T}_{\text{ini}}$	[K]	Linear	Absolute value	53	33	143
Vent spreading angle	α v ${\alpha }_{v}$	[°]	Linear	Absolute value	0	0	9
Vent latitude	L $L$	[°]	Linear	Relative value	0	−0.32	+0.32
Vent longitude	l $l$	[°]	Linear	Relative value	0	−0.32	+0.32
Vent zenith angle	θ Z ${\theta }_{Z}$	[°]	Linear	Relative value	0	−0.8	+0.8
Vent azimuth angle	θ A ${\theta }_{A}$	[°]	Linear	Relative value	0	−0.8	+0.8

• Note. The fourth and fifth columns refer to the scales used to calculate the random values: linear or logarithmic scales and absolute or relative values. The relative values are calculated relative to the location and orientation of the vents, as reported by Porco et al. (2014).

2.2 Expansion of the Plume at Higher Altitude

We compute the further expansion of the plumes toward space using the analytical curve fit expressions for the flow fields at an altitude of 10 km as the lower boundary condition. We define a rotating reference frame with the origin at the center of Enceladus and the x $x$ -axis pointing toward Saturn—because Enceladus is tidally locked in its orbit around Saturn, the z $z$ -axis pointing toward Enceladus’ North, and the y $y$ -axis such that the axis system is right-handed. We assume here that (a) Enceladus’s orbital plane and Saturn’s equatorial plane coincide (the inclination angle between the two planes is negligible at 0.019°), (b) the orbit of Enceladus around Saturn is circular—the orbital eccentricity is extremely small (0.0047), (c) the rotation axis of Enceladus is normal to its orbital plane, and (d) Saturn is a spherical body.

In the above-defined Enceladus rotating reference frame, the particles are subject to different accelerations: g E ${\boldsymbol{g}}_{\boldsymbol{E}}$ and g S ${\boldsymbol{g}}_{\boldsymbol{S}}$ , relative to the Enceladus and Saturn gravity fields, a E ${\boldsymbol{a}}_{\boldsymbol{E}}$ the centrifugal acceleration due to the orbit of Enceladus around Saturn, and the Coriolis and centrifugal accelerations due to Enceladus’ self-rotation ω E ${\boldsymbol{\omega }}_{\boldsymbol{E}}$ .

If r = ( x , y , z ) $\boldsymbol{r}=(x,y,z)$ denotes the position of a particle, e x , e y , e z $\left({\boldsymbol{e}}_{\boldsymbol{x}},{\boldsymbol{e}}_{\boldsymbol{y}},{\boldsymbol{e}}_{\boldsymbol{z}}\right)$ are the unit vectors of the right-handed coordinate system, and e r ${\boldsymbol{e}}_{\boldsymbol{r}}$ the unit radial vector, we can write the equation of motion of a particle (Kempf et al., 2010; Yeoh et al., 2017):

r ¨ = g E + g S + a E − 2 ω E × r ˙ − ω E × ω E × r $\ddot{\boldsymbol{r}}={\boldsymbol{g}}_{\boldsymbol{E}}+{\boldsymbol{g}}_{\boldsymbol{S}}+{\boldsymbol{a}}_{\boldsymbol{E}}-2,\left({\boldsymbol{\omega }}_{\boldsymbol{E}}\times \dot{\boldsymbol{r}}\right)-{\boldsymbol{\omega }}_{\boldsymbol{E}}\times \left({\boldsymbol{\omega }}_{\boldsymbol{E}}\times \boldsymbol{r}\right)$ (1)

where g E = − μ E r 2 e r ${\boldsymbol{g}}_{\boldsymbol{E}}=-\frac{{\mu }_{E}}{{r}^{{2}},{\boldsymbol{e}}_{\boldsymbol{r}}$ is the gravitational acceleration from Enceladus, with μ E ${\mu }_{E}$ the product of the universal gravitational constant and the Enceladus mass ( μ E = 7.21 × 10 9 m 3 / s 2 ${\mu }_{E}=7.21\times {10}}{9},{\mathrm{m}}^{{3}/{\mathrm{s}}}{2}$ ), g S = − μ S ‖ r − d E S e x ‖ 2 r − d E S e x ‖ r − d E S e x ‖ ${\boldsymbol{g}}_{\boldsymbol{S}}=-\frac{{\mu }_{S}}{{\Vert \left(\boldsymbol{r}-{d}_{ES},{\boldsymbol{e}}_{\boldsymbol{x}}\right)\Vert }^{{2}},\frac{\left(\boldsymbol{r}-{d}_{ES},{\boldsymbol{e}}_{\boldsymbol{x}}\right)}{\Vert \left(\boldsymbol{r}-{d}_{ES},{\boldsymbol{e}}_{\boldsymbol{x}}\right)\Vert }$ is the gravitational acceleration due to Saturn, where d E S ${d}_{ES}$ is the distance between Enceladus and Saturn ( d E S = 238 , 020 km ${d}_{ES}=238,020,\text{km}$ and μ S = 3.798 × 10 16 m 3 / s 2 ${\mu }_{S}=3.798\times {10}}{16},{\mathrm{m}}^{{3}/{\mathrm{s}}}{2}$ is the product of the Enceladus mass and the universal gravitational constant), a E = − 2 π τ E 2 d b E S e x ${\boldsymbol{a}}_{\boldsymbol{E}}=-{\left(\frac{2\pi }{{\tau }_{E}}\right)}^{{2},{d}_{{b}_{ES}},{\boldsymbol{e}}_{\boldsymbol{x}}$ is the centrifugal acceleration due to the orbit of Enceladus around Saturn, where τ E = 1.18 × 10 5 s ${\tau }_{E}=1.18\times {10}}{5},\mathrm{s}$ is the Enceladus orbital period and d b E S ${d}_{{b}_{ES}}$ is the distance between the Saturn-Enceladus system barycenter and the Enceladus mass center, which is nearly equal to d E S ${d}_{ES}$ (relative difference of the order of 10−7), 2 ω E × r ˙ $2,\left({\boldsymbol{\omega }}_{\boldsymbol{E}}\times \dot{\boldsymbol{r}}\right)$ is the Coriolis acceleration, where ω E = ω E e z ${\boldsymbol{\omega }}_{\boldsymbol{E}}={\omega }_{E},{\boldsymbol{e}}_{\boldsymbol{z}}$ and ω E = 2 π τ E ${\omega }_{E}=\frac{2\pi }{{\tau }_{E}}$ , and ω E × ω E × r ${\boldsymbol{\omega }}_{\boldsymbol{E}}\times \left({\boldsymbol{\omega }}_{\boldsymbol{E}}\times \boldsymbol{r}\right)$ is the centrifugal acceleration due to Enceladus’ self-rotation. Even though we include all force terms in the simulations described below, the terms g E ${\boldsymbol{g}}_{\boldsymbol{E}}$ and 2 ω E × r ˙ $2,\left({\boldsymbol{\omega }}_{\boldsymbol{E}}\times \dot{\boldsymbol{r}}\right)$ are the more important ones.

Appendix B in Supporting Information S1 provides details on the analytical and free-molecular methods. It also provides two applications to simple cases, one for INMS observation configuration and one for UVIS, and a comparison between the analytical solution and the solution from a free molecular code.

3 Vent Parameter Sensitivity Study

3.1 Introduction

In this work, we focus exclusively on the 98 vent positions and orientations reported by Porco et al. (2014). These vents were identified using high-resolution Cassini ISS imaging, with their locations and directions triangulated from resolved particulate jets. This catalogue remains the most spatially complete and geometrically constrained data set available for the Enceladus South Polar Terrain. Alternative indicators of vent activity—such as thermal anomalies from CIRS (Howett et al., 2011; Spencer et al., 2006) or compositional signatures from INMS or CDA (Magee and Waite., 2017; Waite et al., 2006, 2017)—offer complementary insights but typically lack either the spatial coverage, positional accuracy, or directional information required for the present modeling approach. Using this uniform and well-characterized vent set allows us to consistently apply our sensitivity and inversion framework across all observations.

Our final goal is to obtain, for as many of the 98 geysers as possible, the vent parameters listed in Table 1 to which an observation geometry is sensitive, that is, that could be further constrained. We note that the logarithm of the ice grain radius is fitted, as this form was also used in the parametrization of the flow fields presented in Mahieux et al. (2019). Also, if sensitive to these variables, the water vapor and water ice speeds are fitted separately, as their respective speeds were reported to be different (Ershova et al., 2024; Kempf et al., 2010; Postberg et al., 2011; Schmidt et al., 2008), even though their respective boundaries in Table 1 are equal. We also consider potential adjustments to the latitude/longitude of the vents and the azimuth and zenith angles of their pointing directions listed in Porco et al. (2014), according to the uncertainties provided in Table 2 of the reference: we allow the vent latitude and longitude to vary by ±0.32° and the azimuth and zenith angles by ±0.8° from the reported values. In the remainder of the text, we include the location and orientation of each vent under the designation of “vent parameter.” That is, we have a list of 12 parameters to consider per vent, or 1,176 free variables in total for the 98 vents. Each INMS or UVIS observation has between 10 and 125 measurement points sampled at different locations along the Cassini flight path, and the UVIS occultation has ∼100 column-integrated values, leading to under-constrained fits.

Table 2. Number of Vent Parameters by Category for Each Observation Geometry

Instrument Observation Number of observation points Mass flow rate Gas speed Ice grain speed Vent radius Ice grain diameter Mass ratio Spreading angle Temperature Latitude Longitude Azimuth Zenith Total number of variables INMS E3 17 17 – – – – – – – – – – – 17 E5 10 10 – – – – – – – – – – – 10 E7 37 37 – – – – – – – – – – – 37 E14 125 96 – 1 6 4 1 – 3 2 6 3 3 125 E17 105 85 – 1 5 3 1 – 3 2 2 1 2 105 E18 119 90 – 1 6 4 1 2 3 1 7 2 2 119 UVIS SO 105 98 – – – – – – 2 – 2 2 1 105

The latter motivates the implementation of specific methods to determine the vent parameters that likely contribute to a particular observation and should be considered in a subsequent fitting procedure. To reduce the size of the problem, since a single plume computation only takes a few seconds, we use a Monte Carlo statistical method (Higdon et al., 2017) and calculate correlation probabilities using the Pearson algorithm, available in MATLAB and thoroughly described in Appendix C in Supporting Information S1, to a given Quantity of Interest (QoI). For the QoI, we use the simulated water vapor number density at each INMS measurement point along the trajectory or the integrated water vapor column density at each UVIS measurement point. We note that this work does not include any contribution from material originating from the E-ring, which adds a diffuse continuum to the measured signal of the various instruments (see Postberg et al. (2009), Tenishev et al. (2024), and Yeoh et al. (2015) e.g.).

In this section, we apply the Pearson correlation coefficient method to the vent parameters without considering the physical dependence of these parameters on the flow conditions. An essential advantage of this method, as compared to methods that vary one parameter at a time while considering all other parameters constant, is that one can isolate the relative importance of a given parameter within a complex system driven by many parameters while effectively canceling the effects of all the others—since all parameters are varied simultaneously and randomly. We vary each parameter over the ranges given in Table 1 and arbitrarily set a threshold value to 0.01 to determine whether an observation is sensitive to a given vent parameter. The ranges were designed to encompass the broadest possible physically plausible values for each parameter as discussed by Mahieux et al. (2019). We illustrate our approach, discuss its validity, and examine the sample size needed to correctly calculate the correlations in a simple case, considering simplified observation geometry in Appendix C in Supporting Information S1.

Given the high dimensionality of the problem and the limited number of observational constraints, multiple parameter reduction strategies are possible. In this work, we chose to fit only those vent parameters that showed meaningful sensitivity to the specific observation geometry, while fixing the remaining parameters to values derived from Yeoh et al. (2015). This choice was motivated by numerical stability considerations for the fitting algorithm presented in Section 4 as well as the desire to preserve local variability in parameters where sufficient sensitivity exists. Alternative strategies, such as assuming a single global value for certain vent parameters (e.g., grain size or exit speed), could provide a more constrained framework but may mask localized physical differences or temporal variability. Since previous studies have shown time variability in mass flow rate associated with tidal forcing, it is reasonable to expect other vent properties to vary as well. Our approach aims to capture these potential differences without overfitting the available data. Future work may compare alternative parameter reduction strategies to assess the robustness of the retrieved plume characteristics.

3.2 Sensitivity to the 98 Geyser Vent Parameters

We consider the actual orbit configurations of six INMS observations (E3, E5, E7, E14, E17, and E18) and the UVIS solar occultation, together with the vent positions and orientations given by Porco et al. (2014).

We present the vent parameter correlation coefficients obtained using the abovementioned Pearson method for these observations. The total number of vent parameters for each observation geometry corresponds to the number of non-zero observation points to avoid under-constraining the problems, as reflected by the equal numbers provided in columns “Number of observation points” and “Total number of variables” for each observation in Table 2.

The mean correlation values for the INMS observations are depicted graphically in Figure 1 for E14 and in Appendix C in Supporting Information S1 for the other observation geometries; see Figures C2 to C6. The mean correlation values are also tabulated in Appendix C in Supporting Information S1 for all observation geometries; see Tables C2 to C8. The number of correlated parameters for each observation—those whose mean correlation value exceeds 0.01—is summarized in Table 2.

South polar view of the mean correlation of the vent parameters to the Ion and Neutral Mass Spectrometer E14 observation geometry, where each symbol represents the position of the vents reported by Porco et al. (2014), and the color code is the mean correlation value (see color bar). The black line represents the projection on the surface of the Cassini spacecraft trajectory. The parameters correlated with values greater than 0.05 are triangles, lower are circles, and the uncorrelated parameters are black points.

As an example, we consider the sensitivity to observation E14, presented in Figure 1, which had its closest approach at an altitude of 95 km and a latitude of −89°, with the spacecraft path nearly parallel to the main direction of the Tiger Stripes in a Southwest-Northeast configuration. Ninety-six vent parameters show a confident correlation (with a correlation value larger than 0.01) to mass flow rate, one to ice grain speed, six to vent radii, four to ice grain diameter, one to the water phase mass ratio, three to the exit temperature, and 14 to either location or orientation components. Thus, we can now rank the vent parameters that contribute the most to mass flow rate, vent radius, ice grain exit speed, etc., for each observation to the measured INMS local water vapor number density.

The mean correlation coefficients for mass flux decrease with lateral distance from the spacecraft ground track, as expected from the rapid radial fall-off of plume density: vents nearest the line of sight contribute most strongly to the measured signal. In contrast, parameters such as the ice-to-vapor mass ratio or grain exit speed show isolated non-zero correlations only where the trajectory skims the edge of a vent’s plume cone; these kinds of hot spots mark locations of locally enhanced sensitivity, whereas most vents remain unconstrained. Finally, comparing flybys E14 (Figure 1) and E7 (Figure C4 in Supporting Information S1) illustrates the role of trajectory orientation: E14’s ground track lies nearly parallel to the Tiger Stripes and therefore samples multiple vents along their principal axis, resulting in broad sensitivity to both mass flux and vent radius. In contrast, E7 traverses the stripes nearly perpendicular to their strike, sampling each plume in cross-section and producing a more localized and parameter-specific correlation pattern. These differences underscore how the spacecraft’s path relative to the fracture geometry controls which vent parameters can be constrained.

Regarding the other INMS observations, similar conclusions can be drawn. The spacecraft flew far from the satellite for the E3 (Figure C2 in Appendix C in Supporting Information S1) and E5 (Figure C3 in Supporting Information S1) observation geometries. The closest approach occurred at an altitude of 49.1 km and latitude of −20° for E3 and 2525 km and −2828° for E5; thus, both were far from the South Polar Region. The minimum distance between the spacecraft and the South Pole was 225 km for E3 and 162 km for E5 due to a large angle relative to the Enceladus equatorial plane of the Cassini orbit. In addition, the mean speed of the spacecraft with respect to Enceladus for E3 and E5 was also larger than during E7: 14.4 and 17.7 km/s, respectively, which limited the number of sampling points. Consequently, INMS sampled the geysers far away from the moon and at a relatively large speed, resulting in a poorer data set. Consequently, the mass flow rate of 17 vents is correlated with the observation geometry for E3 and of 10 vents for E5. INMS E7 presents similar results. It had its closest approach at an altitude of 96 km and a latitude of −87°, with the spacecraft path following a South West–North East ground track that intersected all four Tiger Stripes at a steep crossing angle, rather than running parallel to them. The mean speed of the spacecraft with respect to Enceladus during the measurement period was 7.7 km/s. The observation geometry is sensitive to the mass flow rates of 37 vents. We cannot infer information about the other vent parameters from these three observations.

The E17 and E18 observation geometries are similar to that of E14, as the spacecraft was flying on an orbit nearly parallel to the Enceladus Equatorial plane, at a minimum distance from the satellite of 72 km for E17 and 71 km for E18, located close to the South Pole (−87° and −88° respectively) and at low relative mean speed (∼7 km/s). They also occurred on a flight path nearly parallel to the main direction of the Tiger Stripes. Similar conclusions can be drawn for E14: correlations to the mass flow rate of many vent parameters are obtained (85 and 90, respectively), and to some other vent parameters, 19 for E17 and 29 for E18. Besides being sensitive to the mass flow rates, we see that some information relative to nearly all other vent parameters (besides the water vapor exit speed) could theoretically be obtained from remote sensing observations, retrieving only information on the water vapor number density.

Finally, we consider the UVIS solar occultation observation, which sampled the whole plume system, measuring column-integrated water vapor densities. The results related to the sensitivity study of the vent parameters are presented in Figure C7 in Appendix C in Supporting Information S1. The Sun was located to the Northwest, while the spacecraft was orbiting in the Southeast. The sensitivity analysis returns correlations to all 98 vent mass flow rates. This is related to the fact that during the solar occultation, due to the sampling technique, UVIS accounted for the contribution of all geysers, in contrast to INMS observation, which sampled local densities. In addition, the observation was sensitive to the exit temperature of two vents and the location and orientation of five vents.

It is important to note that given the large number of vent parameters and the limited number of measurement points per observation, the system is inherently under-constrained. While the Pearson correlation method allows us to identify potentially sensitive parameters, weak correlations—especially those limited to a few vents far from the observation path—may not always reflect robust physical constraints. Rather, they may arise from statistical variability or complex interactions within the high-dimensional parameter space. To mitigate this, we apply a minimum threshold of 0.01 to the mean correlation coefficient and consider only parameters whose correlations are statistically significant (p-value <0.05). We interpret strong and spatially coherent correlations as likely indicators of physical sensitivity, particularly when they appear consistently across neighboring vents or multiple observations. In contrast, isolated weak correlations, especially for parameters other than mass flow rate, are presented as tentative indications of potential sensitivity and should be interpreted with caution. Nevertheless, for some observations that have more observation points, notably E14, E17, E18, and UVIS, we find meaningful sensitivity to a range of parameters beyond the mass flow rate, such as vent radius, exit temperature, and orientation, suggesting that these properties can, in favorable configurations, be partially constrained by remote sensing or local sampling of the water vapor number density fields.

The results presented in this section are based on Cassini observation geometries, which were limited in spatial and temporal coverage. In the following section, we explore how the same methodology could be used prospectively to support observation planning for future missions, where flyby geometries and measurement strategies can be optimized.

3.3 How This Approach May Serve as an Operational Tool for Future Space Missions That Involve Observations of Cryovolcanic Activity

The approach presented in this section provides the ability to compute the number density and velocity fields of cryovolcanic plumes rapidly and precisely up to hundreds of kilometers above the surface and to estimate the resulting sensitivity to the vent parameters. A reduction in the computation time by four orders of magnitude compared to DSMC and to free molecular simulations has been demonstrated, while the accuracy of the simulations degraded generally by less than 5%—and up to 20% in the worst-case scenarios—relative to full DSMC and free molecular simulations.

Future missions to Enceladus or Europa (Dayton-Oxland et al., 2023) or other icy satellites showing geysers will need to plan flyby and orbital trajectories that best sample and characterize those geysers based on presumed or known vent locations. The approach presented herein could be used during mission planning to predict what information could be constrained based on particular observation geometries and measurement techniques. This would include in situ observations with instruments similar to INMS or remote sensing observations such as occultations as observed by the UVIS instrument. If a particular vent radius, mass flow rate, or other property were of specific interest, the sensitivity approach described above may determine suitable trajectories, pointing directions, and sampling rates for the mission. Similarly, a related sensitivity study employing the same techniques could be undertaken if vehicle safety were an issue for a space mission to a satellite producing a substantial gas or particulate flux (e.g., Io). Reliable sensitivity maps could be quickly created for possible observation geometries (including spacecraft path, speed, and sampling rates) to maximize the mission’s scientific outcome. We note that the approach provided in this section is limited to plumes on airless bodies where the gravity is too low to produce canopy shocks. Thus, it would also apply to the putative plumes that may exist on Europa (Dayton-Oxland et al., 2023).

4 Constraining the Vent Parameters

Here, we consider specific retrieval fitting techniques for particular vent conditions. We present the results of the sensitive vent parameters that fit the data measured by INMS and UVIS. To do so, we follow the procedure described in Appendix D in Supporting Information S1, which uses the Levenberg-Marquardt least-squares algorithm (Dentler, 2024). This procedure returns the maximum number of fitted parameters that do not under-constrain the problem to be solved and a robust evaluation of the uncertainty on each fitted parameter. We note that our original objective was to constrain the vent parameters using the Markov Chain Monte Carlo method, but this did not yield satisfactory results. This effort is summarized in Appendix E in Supporting Information S1.

We order the free vent parameters from five to the number of independent measurement points based on their sensitivity to the problem geometry. We perform a series of parametric fits while increasing the number of free vent parameters. We control the fit’s root mean square (RMS) error and the total mass flow rate. Initially, these two quantities do not vary much from one set of vent parameters to the next; see Figure D1 in Appendix D in Supporting Information S1. Then, at some point, both quantities start to differ while the RMS error increases. The third column of Table 3 gives the number of successfully fitted parameters. Next, for each set of fitted variables, we evaluate the dependence on the a priori values used for the fit and the effect of the non-fitted variables on the quality of the fit. First, we check that the fit results are independent of the a priori values. Second, we compute one hundred fits in which the vent parameters that are not fitted for the active vents, that is, for which we fit the mass flow rate, are set to random values. We evaluate the variance in the resulting fitted parameters and add it to the uncertainty of the fitted free vent parameters. The way these effects are accounted for in the resulting fits is described in Appendix D in Supporting Information S1. This approach ensures that (a) the problem is not under-constrained, (b) that the total mass flow rate is consistent and robust since it varies little from one set of free vent parameters to the next one, (c) the computations are done based on robust and extensive DSMC computations that were parametrized in Mahieux et al. (2019), and (d) that the possible effect of the vent parameters that were not fit because of their weak influence on the observation (reflected in the sensitivity values), is accounted for in the procedure.

Table 3. Total Fitted Mass Flow Rate for the Six Ion and Neutral Mass Spectrometer Observations and the Ultraviolet Imaging Spectrograph Solar Occultation and Number of Fitted Vents’ Mass Flow Rates

Observation	Total mass flow rate	Number of free vent parameters	Number of fitted active vents
INMS E3	135 ± 34 kg/s	17	17
INMS E5	204 ± 37 kg/s	10	10
INMS E7	109 ± 15 kg/s	37	37
INMS E14	426 ± 54 kg/s	82	76
INMS E17	486 ± 50 kg/s	69	66
INMS E18	403 ± 13 kg/s	51	51
UVIS Solar Occ.	217 ± 24 kg/s	103	98

The results of the simulations using the fitted p