<figure><img src="image_1.png" /></figure>

# Variations in plume activity reveal the dynamics of water-filled faults on Enceladus

Received: 21 August 2023

Accepted: 13 August 2024

Published online: 28 August 2024

<figure><img src="image_2.png" /></figure>

Check for updates

Ondřej Souček Ⓔ<fnref n="1" />✉, Marie Běhounková Ⓔ², Martin Lanzendörfer Ⓔ³, Gabriel Tobie Ⓔ⁴ & Gaël Choblet Ⓔ⁴

After discovering a jet activity near the south pole of Saturn's moon Enceladus, the Cassini mission demonstrated the existence of a subsurface water ocean with a unique sampling opportunity through flybys. Diurnal variations in the observed brightness of the plume suggest a tidal control, although the existence and timing of two activity maxima seem to contradict stress analysis predictions. Here, we re-interpret the observed plume variability by combining a 3D global model of tidal deformation of the fractured ice shell with a 1D local model of transport processes within south-polar faults. Our model successfully predicts the observed plume's temporal variability by combining two independent vapour transport mechanisms: slip-controlled jet flow and normal-stress-controlled ambient flow. Moreover, it provides a possible explanation for the differences between the vapour and solid emission rates during the diurnal cycle and the observed fractionation of the various icy particle families. Our model prediction could be tested by future JWST observations targeted when Enceladus is at different positions on its orbit and could be used to determine the optimal strategy for plume material sampling by future space missions.

As Enceladus travels around Saturn on its eccentric orbit, the tidal forces raised by the giant planet change periodically. The variations of extensional tidal stress along the south-polar fault system have been expected to control the plume brightness<fnref n="1" />,². However, the observations (see Fig. 1a) revealed two distinct activity maxima around 40° and 200° mean anomalies (MA), with the latter prominent mid-period peak shifted by several hours compared to the predicted maximum of extensional stress, e.g. refs. 3–5 (Fig. 1b, gold colour, with maxima corresponding to ~120° MA or 11h). Complementary to plume brightness observations, the measurements by Cassini's magnetometer⁶, Ultraviolet Spectrometer (UVIS)⁷, and Ion and Neutral Mass Spectrometer⁸, JWST's NIRSpec⁹ and Visual and Infra-red Mapping Spectrometer¹⁰ provided estimates of the vapour mass flux (200–1000 kg/s) with comparatively smaller temporal variations than for the brightness, but these

findings may still be inconclusive due to the scarcity of measurements.

To interpret the observed plume brightness variability, three mechanisms influencing the activity and location of a notable mid-period maximum have been proposed<fnref n="1" />: (i) diurnal eccentricity tides with an apparent more than 5-hour time delay, (ii) eccentricity tides with an in-phase physical libration of 0.8°, and (iii) right-lateral strike-slip displacement at the faults. According to the first mechanism, the viscoelastic behaviour of the shell could account for the timing of the main plume activity. However, this holds true only for a thick and low-viscosity ice shell¹¹. In the case of a realistic shell shape with significant south-polar thinning¹²⁻¹⁵, the delay associated with a viscoelastic response becomes inadequate to explain more than 5-hour observed lag between the timing of the extensional stress maximum and the main activity peak¹⁶,¹⁷. The second proposed explanation of the lag in

<footnote n="1">¹Mathematical Institute, Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic. ²Department of Geophysics, Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic. ³Institute of Hydrogeology, Engineering Geology and Applied Geophysics, Faculty of Science, Charles University, Prague, Czech Republic. ⁴Laboratoire de Planétologie et Géosciences, UMR 6112, Nantes Université, Univ Angers, Le Mans Université, CNRS, Nantes, France. ✉e-mail: ondrej.soucek@mff.cuni.cz</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | The diurnal variations of Enceladus' plume activity compared to modelled variations of stress and slip. **a** The diurnal variations of Enceladus' plume activity, with the variations on a longer time scale filtered out (see Methods and Supplementary Fig. 2). The symbols represent data groups associated with measurements obtained during one observational day, while colour-coding reflects the years of data acquisition. The main features of the activity curve involve two peaks (-40° MA and -200° MA, respectively) as well as a ramp (-80°-140° MA). We speculate the presence of a late-period peak for specific, more recent, observations. **b** Shows, for comparison, the modelled¹⁹ evolution of the normal stress σ<sub>n</sub> (gold) and the slip u<sub>slip</sub> (red), i.e., relative horizontal displacement of the opposing faults'</figcaption></figure>

<figure><img src="image_4.png" /></figure>

walls, averaged along the faults as a function of Enceladus' position on its orbit. The timing of slip's maxima and minima is correlated with the positions of both early-period and mid-period peaks, while the timing of maximal tension stress aligns with the position of the ramp in the plume activity between the two maxima. **c** Shows the spatial distribution of the slip along the faults at the specific orbital positions associated with the activity peaks correlating with maxima and minima of the averaged slip (depicted as stars). The distribution of normal stress is illustrated for mean anomalies corresponding to the maximum and minimum averaged normal stress (represented by triangles). Source data are provided as a Source Data file.

terms of libration was also dismissed as it was shown to require an abnormally large libration amplitude compared to the measured one¹⁸. The last option, strike-slip motions, was first rejected in ref. 4 due to a reported lack of fit to the observations. Pleiner Sládková et al.¹⁹ noted a correlation between the timing of the two activity maxima and the occurrence of slip maxima and minima. This time correlation was further corroborated by Berne et al.²⁰. In this study, we encompass the influence of slip-induced aperture as well, exploiting advancements in our understanding of Enceladus' structure and advanced modelling approaches reached during the last decade.

The above analyses based on the tidal stress field have also been complemented by models of more localised processes describing the faults' dynamics, such as thermodynamic conditions within the water-filled and vapour parts of the faults, e.g. ²¹⁻²⁴, or local displacement at the faults¹⁹. While all these previous models shed light on various important aspects of local physical interactions within the fault system, none of them offers a fully comprehensive explanation for the observed temporal pattern in plume activity.

In this work, we combine the global and local viewpoints. We incorporate models of localised vapour and water transport processes within the fault zones, forcing them by displacements and normal stresses computed using a global three-dimensional tidal deformation model of Enceladian ice shell¹⁹. Observing that the maxima and minima of lateral slip at the faults are strongly correlated with the plume activity peaks (Fig. 1b, c), we demonstrate in the following that a slip-controlled mechanism offers the primary control on the plume diurnal activity. By exploring a range of parameters, we derive a more complete description of the physical processes leading to the activity

modulation and provide physical constraints on the dynamics of the water-filled faults. Our main objective is to develop a consistent physical model explaining the main characteristics of the activity curve shown in Fig. 1a: the position of the two main peaks and their relative amplitude, the ramp between the two peaks, and the late-period peak.

## Results

Our modelling approach is schematized in Fig. 2: we first compute the kinematic and dynamic quantities using the 3D model of Enceladus' shell with the south-polar faults following ref. 5 and employing the shape model of ref. 14 (Fig. 2a). We then average the computed quantities over the fault areas (see Supplementary Fig. 1), and impose these averages as forcing for a 1D model of transport processes within the faults. This 1D model comprises the following three components (see Fig. 2b): (I) a bottom water channel connected to the internal ocean (inspired by ref. 23), (II) an upper vapour channel above (inspired by ref. 24) that is connected to the surface via (III) an effective surface aperture model combining two vapour/grains transport mechanisms (Fig. 2c, d)).

Consistent with previous studies²³,²⁴, we consider two characteristic widths for the water and vapour channel regions, of the order of metres and centimetres, respectively. This implies a narrowing of the fault zones towards the surface, as classically assumed²⁵. The hydraulic model for the evolution of the water table (I) forced by tidal motion of the channel walls (Equation (12) in Methods) is based on the application of hydraulic Darcy-Weisbach equation (a turbulent variant of Darcy's law). The analytical vapour-transport model (II) is based on the model by Nakajima & Ingersoll²⁴ and provides the integrated surface
<figure><img src="image_5.png" /></figure>

**Fig. 2 | Our activity model for diurnal variations combines a global and local component.** The 3D global model (a) describes the tidal deformation of a fractured ice shell with realistic thickness variations¹⁴ and surface geometry of the faults based on ref. 5. Kinematic and dynamic quantities are averaged along the faults and over the depth across the shell and used as input for the local model. The local 1D model (b) combines a model for liquid water hydraulic motion (I, white), and a model for the transport of vapour and solid grains above the water table (II, blue) with a near-surface aperture model (III, brown) comprising two mechanisms: a slip-controlled mechanism corresponding to macroscopic jet flow (c) and a normal-stress-controlled mechanism corresponding to ambient (more diffuse) vapour flow

(d). Numbers indicate the typical magnitude of the various length scales. The red dashed line (c) and the gold line (d) correspond to the periodic component of slip ($u_{\text{slip}}(t)$) and the normal stress ($\sigma_n(t)$), respectively, both outputs from the global 3D model (a). The red solid line (c) corresponds to the total slip ($u_{\text{slip}}(t) + u_{\text{slip}}^0$). Also indicated are the main controlling parameters in the model: the global scale effective friction at faults ($\mu_f$), slip offset at faults characterising the first aperture mechanism ($u_{\text{slip}}^0$), the effective Young's modulus characterising the second aperture mechanism ($E_{\text{eff}}$) and the activity scaling factors for the two mechanisms ($A_u^0$ and $A_\sigma^0$).

mass flux, solving for mass, momentum, and energy balance of a 1D vapour column originating at the water table and expanding to a vacuum at the surface of Enceladus. Finally, our effective surface aperture model (III) characterises the vapour/grains transmissivity of the fault zones considering two mechanisms.

The first mechanism (Fig. 2c), which is referred to as *slip-controlled* results from a combination of tidally driven slip and geometrical incompatibilities at the faults²⁶. By creating macroscopic voids that evolve during the tidal cycle, it provides a relatively direct pathway from the water table to the surface, accounting for eruptive jets of vapour/grains at the surface. Simultaneously, the perpetual strike-slip motion prevents the closure of active cracks by continuously breaking any ice freshly deposited on the crack walls. We characterise this mechanism in our model by an effective aperture ($\delta_u$). While primarily controlled by slip at the faults, this aperture is also partially modulated by normal stresses, as outlined in Equation (4) in the Methods section.

Furthermore, we consider an additional vapour/grains transport mechanism, referred to as *normal-stress-controlled* that may reflect transport through side cracks, isolated vents, and/or porous or damaged material within the fault zone and accounts for a possibly more diffuse eruptive component of the vapour flow (Fig. 2d). We assume that this second mechanism can be characterised by an independent aperture ($\delta_\sigma$), entirely controlled by the normal stress within the fault zone and the mechanical properties of the material in terms of an effective Young modulus ($E_{\text{eff}}$), as described in Equation (7) in the Methods section.

The combined effect of the two mechanisms mentioned above (see Equations (9–10) in the Methods section) is represented by the total integrated surface vapour/grains mass flux. In our analysis, we assume that the ratio between the solid particles and vapour remains constant in time for each transport mechanism. However, since we presume that the normal-stress-controlled mechanism represents a less direct and more tortuous pathway, we anticipate that it may

potentially exhibit a smaller ice grain-bearing capacity compared to the slip-controlled one. Thus we assume that the scaling factors between the activities (i.e., solid particle fluxes) and the vapour mass fluxes differ for the two mechanisms.

The constructed plume activity model involves five free parameters: an offset of the slip-controlled aperture mechanism $u_{\text{slip}}^0$, the effective Young's modulus $E_{\text{eff}}$ (characterising the normal-stress-controlled aperture mechanism), along with amplitudes of the two mechanisms $A_u^0$ and $A_\sigma^0$. The fifth parameter is the effective friction coefficient $\mu_f$, used in the 3D tidal deformation model, which characterises the friction at faults. Its value influences the kinematic and dynamic quantities (e.g., slip and normal stress) used as inputs in our 1D model analysis. We consider values $\mu_f = 0.01, 0.05, 0.1, 0.2, 0.4$, and $0.8$, whose limits have been confirmed¹⁹ to approach the frictionless case (for the lowermost value of $\mu_f$) and the case without mobile faults (the uppermost value). As the variability in solid particle flux is much better constrained from brightness curves, we use only these data as a proxy of eruption activity. We fit data compiled by Ingersoll and coworkers²⁷, which corresponds to the longest published data set covering the entire diurnal tidal cycle. We re-scale the data to suppress the long-period variability in plume activity to focus only on the diurnal variability (Fig. 1, see Methods, subsection Data processing, and Supplementary Fig. 2). Once the model parameters are adjusted to reproduce the activity curve as depicted in Fig. 3, we predict the vapour flux based on this calibrated model, as explained below.

## Plume activity fit

For a given value of $\mu_f$, we first calculated the stress and displacement along the faults using the 3D model¹⁹. Subsequently, these values have been spatially averaged over the fault zone. Utilising these averaged quantities as input, the four free parameters governing the 1D activity model were determined through a Bayesian inversion process (employing Markov chain Monte Carlo (MCMC) method). Please refer to the Methods section and Supplementary Figs. 3–8 for more details.
<figure><img src="image_6.png" /></figure>

**Fig. 3** | The activity model prediction (black solid line) compared to the rescaled data shown in Fig. 1 (here displayed as grey symbols (same as in Fig. 1a)) reveals details of the Enceladean "plumbing system". The activity signal is composed of the primary slip-controlled (red dashed) and the secondary normal-stress-controlled (gold dash-dotted) components. The former can explain the two peaks at around 40–50° MA and 200° MA, while the latter can explain the ramp ~80–140° and the increase of activity beyond 300° MA. The relative contribution of the two transport mechanisms varies over time, at the early and mid-period activity maxima, the ratio between the slip-controlled and normal-stress-controlled components is roughly 1:2 and 6:1, respectively). Source data are provided as a Source Data file.

To highlight the relative contributions of the various mechanisms in our model, we initially exclude the hydraulic response; this corresponds to a scenario with no change in the water table elevation ($\Delta h = 0$ instead of given by the solution to Equation (12)). The optimal fit was found for the following parameter values: $\mu_f = 0.05$, $u_{\text{slip}}^0 = -0.323$ m, $E_{\text{eff}} = 9.86 \times 10^4$ Pa, $A_u^0 = 2.21 \times 10^6$ m$^{q-p}$, and $A_o^0 = 160$ m$^q$. Figure 3 displays the predicted activity curve for this best-fitting model (black line) together with the processed (see Methods, subsection Data processing) activity data (grey). Two components of the activity that correspond to the distinct aperture mechanisms are also highlighted: the slip-controlled component (red dashed line) and the normal-stress-controlled component (gold dash-dotted line).

**Two-peak activity due to slip-controlled jet flow.** The first characteristic feature of the plume activity resides in the observation of two maxima: a major mid-period maximum ~200° MA and a less pronounced early-period maximum located ~40° MA. Our model shows that these are determined by the primary slip-controlled component of the aperture model (characterised by $\delta_u$). Their timing depends

<figure><img src="image_7.png" /></figure>

**Fig. 4** | The sensitivity of modelled activity to the main model parameters. a The influence of the primary slip-controlled aperture mechanism, responsible for macroscopic jet flow is twofold: (1) a smaller friction coefficient ($\mu_f$) results in a delayed early-period peak and an advanced mid-period peak (green); (2) the role of slip offset ($u_{\text{slip}}^0$) is strictly geometric and affects the relative amplitude of both peaks but not their timing (red). b The influence of the secondary aperture mechanism, controlled by normal stress and responsible for ambient flow, is particularly manifested around 120° MA. The dash-dotted lines indicate the contribution of this secondary mechanism while, as in other panels, the solid lines include the sum of both contributions. c The hydraulic amplification (blue line) results in changes mostly between ~300° and 50° MA, with a local maximum at around 324° MA possibly corresponding to the enhanced activity witnessed for specific observation days (Fig. 1a). Source data are provided as a Source Data file.

primarily on the effective coefficient of friction $\mu_f$, while their relative amplitude is governed by the slip offset $u_{\text{slip}}^0$ (see Fig. 4a and also Supplementary Fig. 9).

The best match of the model prediction with the observation is obtained for very low values of the friction coefficient $\mu_f = 0.05$ (see Supplementary Figs. 3–8 and Supplementary Table 1), consistent with the expected lubrication of the fault by liquid and water vapour$^{28}$. As the friction coefficient $\mu_f$ increases, the early-period activity peak is delayed, while the mid-period peak advances due to the escalating asymmetry in the slip between the compressional and extensional phases of the tidal period. When $\mu_f$ exceeds 0.2, the model reproduces the main activity peak at ~200° MA but fails to capture the first peak (refer to Supplementary Figs. 3–8). This is consistent with Berne et al.$^{20}$ who favour high values of $\mu_f$ larger than 0.3 to provide the best correlation of faults' slip motions with the plume activity but fail to reproduce the shape and timing of the early-period peak. This discrepancy arises from the diminishing significance of the slip-controlled contribution to the overall activity in these cases.

Variations in $u_{\text{slip}}^0$ primarily influence the relative amplitudes of the two peaks, with their positions remaining essentially unchanged. This geometric effect can be explained (in view of Equation (4)) through the dependence of the function $|A \sin(t + \phi) + B|$ on the parameter $B$, with $B$ representing the slip offset $u_{\text{slip}}^0$, while the term $A \sin(t + \phi)$ approximates the time-periodic component of the slip, the amplitude $A$ and the phase of the periodic slip $\phi$ result from the 3D deformation model. The optimal value of $u_{\text{slip}}^0$ (~0.323 m) corresponds to a right-lateral offset of ~1/3 of the slip amplitude (~1 m) for the optimal $\mu_f$. A mechanism hinting at the possible origin of such a slip offset was suggested and modelled in ref. 19, namely a compressive stress
background developed in the SPT of Enceladus due to frictional response of a periodically loaded fault with Coulomb-type friction. Our results are consistent with such an analysis in terms of sign (right-lateral slip), but the amplitude is roughly one order of magnitude larger than what is proposed by ref. 19 for a small friction coefficient. Their model, however, focused primarily on the resulting stress state arising from the inclusion of Coulomb-type friction at the faults, but did not take into account any long-term relaxation processes in the south-polar region caused by this stress state. In addition, also the geometry of the fault system was fixed. Notably, the reference right-lateral slip state appears also in agreement with the interpretation of the tectonic system at the SPT being driven by regional-scale right-lateral shear<sup>29,30</sup>.

**Activity ramp due to normal-stress-controlled ambient flow.** A second major observation corresponds to the ramp of slowly increasing activity between -80° MA and 140° MA, which connects the two peaks. While apparently a secondary characteristic of the activity of the Enceladian plume, it is strongly constrained by Cassini observations (Fig. 1a). Our results demonstrate that this intermediate ramp originates from the normal-stress-controlled aperture mechanism. For the best-fitting model, this mechanism has a single maximum at -120° MA, corresponding to the maximum extensional stress (dash-dotted line in Fig. 3, see also Fig. 1b, c and Supplementary Fig. 1). Two parameters, $E_{\text{eff}}$ and ratio of amplitudes $A_u^0/A_o^0$ control the relative contributions of the two vapor-transport mechanisms, see Fig. 4b and Supplementary Fig. 10. The optimal value of the effective Young modulus $E_{\text{eff}}$ (9.86 × 10⁴ Pa) characterising the secondary normal-stress-controlled mechanism corresponds to 5 orders of magnitude reduction when compared to the value for undamaged ice ($E_{\text{ice}}$ -10 GPa<sup>28</sup>), suggestive of pervasively fractured ice or a highly porous and unconsolidated ice powder environment<sup>31</sup>.

**Hydraulic amplification.** So far, we have not considered the possibility of a significant dynamic water table change<sup>23</sup>. Indeed, our results show that except for a very specific setting, which we refer to as “hydraulic amplification”, the latter is negligible compared to the hydrostatic surface-to-water table distance, making the effect insignificant (see Supplementary Fig. 11). We now consider the hydraulically amplified state in the case of the best-fitting model in order to assess its impact (light blue line in Fig. 4c). In this state, the width of the water channel is such that the imposed normal motion of the channel walls leads to an almost complete closure of the water channel, resulting in the water table rising and then decreasing by hundreds of metres. This would imply that liquid water may reach the surface.

The amplified state of the hydraulic system induces a notable contribution to the plume activity, namely a third (late-period) peak at -320° MA (blue line in Fig. 4c), which is hinted by a slight increase in activity for some observational datasets (Fig. 1a), followed by a less pronounced local activity minimum around the pericenter. The amplitude of this maximum is highly sensitive to the maximal water table elevation. The mid-period maximum is unaffected. The early-period maximum is slightly delayed and shifted from -47° MA to approximately 55° MA. We note that it has been speculated that the lesser early-period activity peak corresponds to a double peak, with two local maxima around 36° MA and 51° MA<sup>27</sup>. This double peak might be in line with the 8° shift predicted by our model in the case where hydraulic amplification affects only a portion of the fissure system. In addition, the activity at the beginning of the tidal cycle (between 0–50° MA) is also suppressed owing to a decreased water table elevation (cf. Fig. 4c).

## Discussion

We showed that the brightness variability curve as observed by Cassini can be explained by combining two complementing mechanisms

<figure><img src="image_8.png" /></figure>

**Fig. 5 | Prediction of vapour mass flux and evolution of Type II particles.** a Predicted vapour mass flux corresponding to the optimal activity fit (shown in Fig. 3) and increasing vapour-to-solid ratio ($r$) between the normal-stress-controlled and slip-controlled surface aperture mechanisms. For $r > 1$ the normal stress-controlled mechanism enhances the vapour mass flux at orbital ranges between -50° MA and 150° MA (a vapour-to-solid ratio of $r = 5$ is displayed as a grey curve and thin grey lines correspond to $r = 2, 3, 4$). A thick orange line highlights the least-squares fit ($r = 2.8$) hinting at a possible explanation of the rescaled UVIS data<sup>7</sup>. The expected influence of the hydraulic amplification, when present, is also displayed (shaded region). b Evolution of Type II particle population on an orbit. Assuming values in range 55–65% at MA = 200°<sup>32</sup> yields various envelopes depending on the ratio measuring the effect of both mechanisms on their emission rates ($c_u^\text{II}/c_u^\text{I}$, changing from 0 to 100%, in colours). The resulting average for Type II population (right of b) varies significantly (between 35% and 65%), which can be compared to the average population observed during flybys E5 and E17<sup>32</sup> denoted by the error bar: 20–40%. Colour-coding of UVIS data and flybys denotes the measurement year (see Fig. 1a, for the corresponding colour bar). Source data, excluding dataset from ref. 7 and ref. 32, are provided as a Source Data file.

controlling the mass flux. The distinctive two-peak pattern is attributed to jets emanating through slip-controlled apertures along the faults. For brevity, we will denote this mechanism as the “primary”. Operating on a finer spatial scale, the ambient flow controlled by the normal stress explains the existence of a more subtle feature, the ramp. We will refer to this mechanism as the “secondary”. These two mechanisms contribute differently to the solid grains flux (activity), and, analogously, they may also yield different vapour fluxes. However, one cannot assume their respective vapour-to-solid fraction to be identical. Namely, transport through more tortuous and narrower pathways envisioned for the secondary mechanism is expected to promote the transmission of vapour compared to solid grains, relative to the primary one that offers a more direct path from the water table. To account for this difference, we introduce the relative vapour-to-solid ratio $r \ge 1$, which imposes that the vapour-to-solid fraction for the secondary mechanism is $r$ times larger than for the primary one (see Methods, subsection Modelling vapour mass flux). By varying the value of $r$, we generated distinct vapour mass flux curves that correspond to the best fit of activity (brightness) curve (see Fig. 5a). When comparing the predicted vapour mass flux with the vapour flux estimate based on the Cassini UVIS data<sup>7</sup>, we obtained the best fit for $r = 2.8$. Consequently, the premise of two distinct surface aperture mechanisms provides a plausible explanation for the apparent differences between the time variability observed in visible-infra-red brightness data and the reduced variability in UV occultation data, which previous models failed to explain<sup>32</sup>.

The UVIS data have been acquired for a limited range of MAs and over a long time period (2005–2017), which makes comparison with our predicted curve activity delicate. However, we notice that the secondary normal-stress-controlled mechanism may result in an enhancement of vapour flux relative to the solid grain flux, thus
explaining the reduced variability observed in the UVIS data between 100° and 250° (Fig. 5). A recent observation by JWST suggests a water vapour flux comparable to the values retrieved from UVIS data⁹. However, it is difficult to make a direct comparison between the two datasets due to differences in observing geometry and the long time span since the UVIS data. Future JWST observations performed with Enceladus at different MAs on its orbit, but acquired during a relatively short period of time will be essential to test our activity model. Based on our prediction, observations should be done in priority for MA values -200° and -310° in order to maximise the variations in water vapour column densities. At the same time, the timing of -320° MA provides the biggest potential for witnessing possible near-subsurface emission of vapour related to the hydraulic amplification.

The two proposed transport mechanisms are also anticipated to contribute differently to the emission of solid icy grains, depending both on their size and varying formation processes. Postberg et al.³³⁻³⁵ have categorised three families of icy grains observed in the E-ring and inside the plume owing to distinct compositions revealed by the Cosmic Dust Analyzer (CDA) onboard Cassini. Type I grains are salt-poor ice grains likely resulting from the condensation of water vapour during the ascent through conduits but also above the surface within the vents, while organic-rich Type II and salt-rich Type III grains are expected to be aerosols directly blasted from the top of the water table. Type II and Type III grains should therefore be more sensitive to the permeability/tortuosity of the conduits above the water table and thus their emission would be reduced compared to Type I grains during periods where the secondary mechanism is dominant, and enhanced during periods of highest primary mechanism activity. A similar fractionation of the grain populations is also expected based on their size as the fraction of nearly pure water ice particles (Type I) appears to decrease with increasing particle size. Again, the secondary transport mechanism would be detrimental to the emission of larger on average Type II and Type III particles.

Among the three solid particle families, organic-rich Type II is the most promising in terms of astrobiology. Based on our model, we predict that the maximal production of Type II should occur between approximately 200° MA and 240° MA when the contribution of the secondary mechanism is minimal, and should be smaller for mean anomalies between 80° and 120° and between 310° and 350° (Fig. 5b). Quantification of this reduction is not straightforward. If the background Type II population observed in the plume during flybys E5 and E17 is used as a guidance for the average population, one should expect large variations in emission rates, corresponding to a concentration of Type II particles transported by the secondary mechanism representing less than 25% of that transported by the primary mechanism. A similar line of reasoning suggests a comparable diurnal variation for Type III ice grains. These “salt-rich” particles amount to 5–10% of the grains analysed in the E-ring, they have been shown to also contain phosphates³⁶ and are generally considered to represent a fresh sampling of the subsurface ocean.

The maximum population of Type III grains (30%) was observed close to the surface during flyby E5, but this value is smaller for E17 and E18. In addition, the value in the core region of the plume, measured at an altitude of several hundreds of km resembles that observed in the E-ring (10%)³⁵. Constraints on the maximum fraction for the Type III population immediately emitted in the plume are thus even more uncertain. Nevertheless, we expect that the largest fraction of Type III grains will also be observed at the MA corresponding to the global mid-period activity peak (around 200° MA, as was coincidentally the case for flybys E5, E17, and E18).

Our model prediction is based on average activity and average fault properties. In reality, the emission rates likely vary depending on the local fault properties. Each jet source may reach a maximal emission rate for a MA slightly different from the average value predicted here. This can be seen for instance by comparing the emission rate

observed by Cassini for flybys E5 and E17, which showed significant variability even if the two flybys occurred at the same MA. Our model shows that there is a clear correlation between the dynamics of water-filled faults and the emission rate, for which at a given time the emission of salt-rich/organic-rich grains is optimal. Careful monitoring of fault dynamics by combining thermal emission mapping, radar sounder, InSar and other geophysical techniques, e.g. ³⁷⁻⁴¹ from a future orbiter platform will permit the identification of the best site and the best timing to sample ejected plume materials at low altitudes. A detailed characterisation of the most active sources from orbital observations will be crucial to identifying the most favourable landing site to collect freshly erupted oceanic materials carrying biosignatures⁴².

## Methods

### Data processing

We use the dataset of Ingersoll and coworkers²⁷ providing the longest published data collection. Ingersoll et al.²⁷ identified several tens of groups based on the days of data acquisition. Due to the temporal variability on longer time scales (mainly 4 and 11 years) and aperiodic contributions, it is challenging to fit the data. Inspired by ref. 27, we have constructed a template from the groups obtained within seven days (groups 1, 2, and 3 from June 2017) due to their global coverage of the orbit and their temporal proximity. For the template, we assumed diurnal periodicity and smoothed the data using the Savitzky-Golay filter⁴³ and its SciPy implementation⁴⁴. Then each of the data groups was rescaled by the factor $A_D$ to the template. The rescaling factor is obtained by the minimisation of misfit S between data $d_i$ (for MA $t_i$) and template T

$$S^2 = \min \left( \sum_{i=1}^{N_D} \left| d_i - \frac{T(t_i)}{A_D} \right|^2 \right), \quad (1)$$

where $N_D$ is the number of measurements in the group D. This procedure allows us to filter out the long-period changes in the plume activity. We denote the rescaled data as $d_i^R = A_D d_i$. The figure graphically representing the data processing flowchart is in the Supplementary Fig. 2.

### Global shell deformation model

For our study, we incorporated modelling data sourced from Pleiner Sládková et al.¹⁹, where we conducted an investigation into the tidal-induced deformation of Enceladus' outer shell, specifically focusing on the role of faults characterised by Coulomb-type friction. We used a finite element model based on the open-source library FEniCS, which was initially introduced by Souček et al. and Běhounková et al.¹⁶,¹⁷,⁴⁵. Within the framework by ref. 19, the frictional faults (based on ref. 5) located in the moon's south-polar region are treated as damaged zones with a width of several hundred metres. These zones traverse the entire shell (with geometry based on ref. 14) and exhibit a specific visco-elasto-plastic rheology designed to replicate the behaviour of a frictional fault with Coulomb-type friction. For the purpose of our study, we exclusively utilise output variables that have been spatially averaged across the shell's thickness and along the faults. The physical and kinematical quantities relevant for this study are the horizontal and normal displacements at the faults ($u_{slip}$ and $u_n$, respectively), normal velocity ($v_n$), and normal traction ($\sigma_n$) at the fault. In Supplementary Fig. 1, we plot these quantities depending on the orbital position (in terms of the MA) and on the effective Coulomb coefficient of friction $\mu_f$, which characterises the frictional contact at the faults and is varied between $\mu_f = 0.01$ to $\mu_f = 0.8$. Previous research by Sládková et al.¹⁹ has demonstrated that these two end-member scenarios, namely the frictionless case and the case devoid of any faults (referred to as the locked state), are accurately represented by these two end-members.
## Vapour-transport model

In order to analyse the evaporation of water vapour from the water table and its subsequent transportation to the surface of Enceladus, we employed a numerical implementation of the one-dimensional thermodynamic vapour-ice grains mixture flow model proposed by Nakajima & Ingersoll<sup>24</sup>. This numerical approach enabled us to characterise the vapour flux analytically, considering two key parameters: the depth of the water table (D), which represents the distance to the surface, and the effective aperture of the fissures (δ). The model<sup>24</sup> encompasses the mass, momentum and energy balance of a vapour column generated by evaporation at the water table, which is expanding to vacuum at the surface. The model also takes into account the mechanical and thermodynamic interaction of the vapour-ice mixture with the channel walls and is capable of predicting the observed vapour flux quantitatively in a certain range of model parameters (for details on the model, see Supplementary Methods, section Vapour flow model).

Our numerical implementation of the model was carefully tested against the numerical results reported in ref. 24 with a very good match (see Supplementary Fig. 12), the slight differences attributed to differences in numerical details. In order to align with the observed mass fluxes of ~200–1000 kg/s<sup>46,47</sup>, we conducted a fitting analysis of the total mass flux emanating from a fissure with a length of 500 km (which serves as an approximate estimation for the total length of all tiger stripes) as a function of both the vapour channel aperture δ, and of the water table depth D within a specific range of these parameters. Fitting first independently the δ and D dependence of F, we propose the following ansatz

$$F \propto \delta^p / D^q. \qquad (2)$$

Considering this function, we performed an independent nonlinear least-squares 2D fit in the range δ = 3–5 cm, D = 300–600 m, localised around the observed values of the surface mass flux. The numerical fit yields exponent values p = 3.74, q = 1.11, (see also Supplementary Fig. 13), where the simulated data points are plotted along with the fit δ'/D' relationship. For all calculations, we used the value s₀ = 0.1 (solid mass fraction at the water table), and we confirmed only a mild variation of the mass flux with this parameter, in agreement with ref. 24.

## Surface aperture model

We additionally examine a parameterisation of the effective aperture δ, pertaining to the vapour channel. Throughout the tidal cycle, the fissures endure both normal and tangential forces. Conventionally, the plume's activity has primarily been interpreted in terms of normal forcing, e.g., refs. 1,4,11. Moreover, Nimmo et al.<sup>4</sup> investigated the possibility of right-lateral-slip origin. Pleiner Sládková et al.<sup>19</sup> and Berne et al.<sup>20</sup> noted a correlation between the timing of the two activity maxima and the occurrence of slip maxima and minima. In this study, we expand this observation by suggesting its physical origin and interpretation. We expect the near-surface section of the fissures to be in a contact-fissure regime. This is due to the fact that the effective hydrostatic compressive pressure in ice (hydrostatic pressure in ice counter-balanced by the water pressure in the fissures) at the depth corresponding to the water table position (i.e., at -1/10 of the local shell thickness), reaches approximately ρᵢgD, i.e., ~50 kPa for a 5 km local shell thickness at the SPT<sup>14</sup>. At the same time, the maximal extensional dynamic stresses reached in our model with frictional faults do not exceed ~25 kPa, see Fig. 1a and Supplementary Fig. 1, indicating that the upper part of the fissure corresponding to the vapour channel is (at least partially) under compression at all times.

We distinguish two different vapour pathways corresponding to a) jet vapour flux in macroscopic voids resulting from slip and geometric incompatibilities between the fissure walls, and b) vapour flux

through ambient medium comprising the gouge and debris within the contacting parts of the fissures, powder snow condensing on channel walls, side cracks, or isolated vents separate from the main fissure system. Both these mechanisms, distinguished in this manuscript as “slip-controlled (primary)” and “normal-stress-controlled (secondary)”, respectively, are visualised in Fig. 2c, d. For parameterising the first mechanism, we employ the analogy from terrestrial rock faults hydraulics. A non-zero aperture for contact-type fissures is created as a result of the slip displacement and geometry incompatibilities (see Fig. 2c). Employing an analytical hydraulic conductivity model for a simplified sawtooth geometry of the fissures by Elsworth et al.<sup>26</sup>, allows one to relate the slip-controlled aperture δ_u to a mutual slip displacement of the fissures u<sub>slip</sub>. Assuming, in the context of our model, fault wavelength larger than the scale of the slip displacement (about 1 m), we obtain a simple relationship<sup>26</sup>

$$\delta_u \propto |u_{\text{slip}} + u_{\text{slip}}^0|, \qquad (3)$$

where now the u<sub>slip</sub> denotes the average mutual horizontal displacement along the faults and we introduced a free parameter u<sub>slip</sub><sup>0</sup> corresponding to slip offset, required since the input data from the 3D shell model have been pre-processed by subtracting the mean values. The slip-controlled aperture (3) is expected to be, to some extent, also modulated by the tidal stresses. Here, we employ the parameterisation by ref. 48, replacing (3) by the following expression

$$\delta_u \propto |u_{\text{slip}} + u_{\text{slip}}^0| \times f(\sigma_n), \qquad (4)$$

with

$$f(\sigma_n) = \tan(\alpha \text{Log}_{10}(-\sigma_c / \bar{\sigma}_n^{\text{tot}})). \qquad (5)$$

Here the parameter σ_c is the uniaxial compressive strength of the ice (using a moderate estimate of 7 MPa<sup>49</sup>) and σ_n<sup>tot</sup> is the vertically averaged total normal stress approximated by

$$\bar{\sigma}_n^{\text{tot}}(t) = \frac{1}{2}(\rho_i - \rho_w) \frac{\rho_i}{\rho_w} gH + \sigma_n(t). \qquad (6)$$

Here σ_n(t) denotes the average of the dynamic normal loading of the fault, see Supplementary Fig. 1, and the first term represents the effective stress within hydrostatically flooded fissures, see ref. 19 (SI). Parameter α is roughly estimated by a constant value α = 2 (corresponding to the choice of a fixed value of joint roughness coefficient approximately JRC<sub>mob</sub> = 2, see ref. 48).

Apart from the above mechanism, we expect that part of the fault zone region is also permeable to vapour flow by a secondary mechanism, not related to strike-slip-induced openings, and responsible to less eruptive ambient vapour flow. Due to incessant mutual grinding of the fissure walls, we could, for instance, expect sections of the fault system to be filled with ice debris gouge behaving as a porous medium. In addition, we might also expect part of the vapour transfer to be facilitated through side cracks or isolated vents. For all these mechanisms combined together, we choose a common simple parameterisation of the corresponding vapour/grains transmissivity in terms of a secondary normal-stress-controlled effective aperture δ_a depending on the normal stress loading σ_n by the following expression

$$\delta_a \propto (1 + \sigma_n / E_{\text{eff}})^+ . \qquad (7)$$

Here E<sub>eff</sub> is an effective Young's modulus of the normal-stress-controlled mechanism medium, and the symbol (+) expresses the positive part of the argument, reflecting that the aperture cannot become negative (preventing self-penetration). Let us note that the simple stress dependence in (7) is motivated by an idealisation of a
vapour channel with a circular cross-section in a medium with Young's modulus $E_{eff}$ in the plane stress approximation$^{50}$.

## Integrated global plume activity

The total vapour/grains mass flux from the SPT is obtained by integrating (2) along the faults:

$$ \mathcal{F}(t) \propto \int_{\text{faults}} \frac{\chi(s)\delta_u^p(s) + (1-\chi(s))\delta_o^p(s)}{(D_0(s) - \Delta h(t,s))^q} ds \quad (8) $$

where $\chi(s) \in (0, 1)$ is the local weight between the two surface aperture mechanisms, and we replaced $D$ in the denominator by $D_0-\Delta h$, taking into account the dynamic water table change $\Delta h$ with respect to a reference (hydrostatic) water table depth $D_0$. Without more detailed knowledge about the spatial distribution of the two mechanisms along the faults, we replace (8) with a weighted average of the two mechanisms. Moreover, for the prediction of the plume activity (slab density), which corresponds to the content of solid ice particles, we will assume that the ice grain-to-vapour ratio is fixed for each of the two mechanisms, but the ratios may differ between them. As a consequence, we parameterise the plume activity $A$ as follows using the scaling factors (amplitudes) for the two mechanisms $A_u^0, A_o^0$, respectively:

$$ A(t) = A_u(t) + A_o(t), \quad (9) $$

with

$$ A_u(t) = A_u^0 \frac{\left( |u_{\text{slip}}(t) + u_{\text{slip}}^0| \times f(\sigma_n(t)) \right)^p}{(D_0 - \Delta h(t))^q}, \quad (10) $$

$$ A_o(t) = A_o^0 \frac{\left( (1 + \sigma_n(t)/E_{\text{eff}})^+ \right)^p}{(D_0 - \Delta h(t))^q}. \quad (11) $$

## Hydraulic water transport model

In order to compute the motion of water in the fissures, i.e., the quantity $\Delta h$ in Equations (10) and (11), we follow a similar strategy as in Kite & Rubin$^{23}$. In contrast to them, we consider kinematic rather than dynamic forcing, i.e., the water table changes in response to the prescribed normal displacement of the ice walls of the channel. For this kinematic forcing, we employ the (spatially averaged) model data from the shell deformation model, in particular, the normal displacement of the fissure walls $u_n$ and its time derivative—the normal wall velocity $v_n$, as shown in Supplementary Fig. 1. Let us note that in doing so, we assume that the shell deformation model provides us with reasonable estimates of these kinematic quantities for the hydraulic part of the fissure. At the same time, we believe it to be inaccurate in the top (vapour transport) section of the fissure, where the fissure walls exhibit maximum compressive stresses and are thus presumably in a contact regime. This is why an independent parameterisation of the surface aperture has been introduced and described in the previous section.

Assuming (as in ref. 23) that the motion of water in the hydraulic channel is characterised by the Darcy-Weisbach Colebrook-White relation$^{51}$ (their Equation (4.10)), we have derived (see Supplementary Methods, subsection Hydraulic model) a nonlinear ODE describing the elevation of the water table $\Delta h$ with respect to the reference hydrostatic level $h_0$:

$$ \dot{\Delta h}(t) = c(t) - v_n \frac{h_0 + \Delta h(t)}{w(t)}, \quad (12) $$

where $w(t) = w_0 + u_n(t)$ is the time-variable width of the water channel, $w_0$ representing the mean value. The parameter $c(t)$ is given for

$v_n(t) \neq 0$ by the solution of

$$ \begin{aligned} & |c(t)|^3 - |c(t) - \frac{v_n(t)}{w(t)}(h_0 + \Delta h(t))|^3 \\ &= -3g\Delta h(t)v_n(t)\left(4\log_{10}\left[\frac{k}{3.7 \times 2w(t)}\right]\right)^2, \end{aligned} \quad (13) $$

where $g$ is the gravity acceleration and $k$ is dimensionless wall roughness parameter. The unique analytical solution is given in SI (see Supplementary Equation (34)). Alternatively, for $v_n(t) = 0$, the parameter $c(t)$ is given by

$$ c(t) = 4\log_{10}\left[\frac{k}{3.7 \times 2w(t)}\right] \sqrt{gw(t)} \times \sqrt{\frac{|\Delta h(t)|}{h_0 + \Delta h(t)}} \operatorname{sgn}(\Delta h(t)). \quad (14) $$

Values of all parameters are provided in Supplementary Table 2. The solution of the water table equation (12) is obtained using the numerical integration routine NIntegrate from Wolfram Mathematica.

## Data fitting

The processed (see Methods, subsection Data processing) dataset was used to identify the free parameters in the expression for the modelled plume activity $A$ given by Equations (9)–(11). As no error estimates are provided for the data by ref. 27, we considered the error of measurement $\sigma_i$ to be 10% of the maximum, i.e., $\sigma_i = 0.1$ for all $i$ for the normalised data. In the inversion routine, we are searching for each considered value of the friction coefficient $\mu_f$ for four parameters $u_{slip}^0, E_{eff}, A_o^0$, and $c = A_u^0/A_o^0$, where the last couple of parameters is equivalent to identifying the two amplitudes $A_u^0$ and $A_o^0$, but admits a more robust inversion.

Time series for $u_{slip}$ and $\sigma_n$ are taken as spatial averages from the numerical simulations$^{19}$, see Supplementary Fig. 1; $\Delta h(t)$ is the result of the hydraulic model and the values of parameters $p = 3.74$ and $q = 1.11$ are based on the vapour-transport model, both described in detail above. The logarithm of the applied likelihood function for $n$ data-points reads

$$ \log p(d^R | t, \sigma, u_{\text{slip}}^0, E_{\text{eff}}, A_o^0, c) = -\frac{1}{l} \sum_{i=1}^{n} \left( \left| \frac{d_i^R - A(t_i, u_{\text{slip}}^0, E_{\text{eff}}, A_o^0, c)}{\sigma_i} \right|^l + \log(\sigma_i^l) \right). $$

A value of $l=1$ has been chosen to suppress the effect of outliers.

To obtain the posterior distribution of models, we employ Bayesian inversion using the MCMC and emcee library for Python$^{52}$ based on the sampling methods proposed by ref. 53; 50 chains and models beyond 50 times the autocorrelation time were used. Besides the likelihood function, we additionally characterise the results using the misfit

$$ \chi = \left( \sum_{i=1}^{n} \left| \frac{d_i^R - A(t_i, u_{\text{slip}}^0, E_{\text{eff}}, A_o^0, c)}{\sigma_i} \right|^l \right)^{\frac{1}{l}}. $$

## Modelling vapour mass flux

We parameterise the vapour mass flux $\mathcal{F}$ in terms of the predicted activity $A$. Accounting for the expected differences between vapour and solid particle transport efficiency for each of the two mechanisms, we choose to express the mass flux of vapour as follows:

$$ \mathcal{F}(t) \propto A_u(t) + rA_o(t), \quad (15) $$

with $A_u(t)$ and $A_o(t)$ given by (10) and (11), respectively. The parameter $r$ reflects the effectiveness of the vapour transport: for $r > 1$, the vapour-
to-solid ratio in the flow realised through the normal-stress-controlled mechanism is *r* times larger than in the slip-controlled one.

## Modelling solid particles populations

In order to predict the distribution of the solid particles into the three families (Types I, II, and III), we define the fraction $f^X(t)$ ($X=I, II, III$) and its variation during one period as:

$$f^X(t) = \frac{c_u^X A_u(t) + c_o^X A_o(t)}{A(t)}. \qquad (16)$$

The coefficients $c_u^X$ and $c_o^X$ denote the fraction of the particle Type X for slip- and normal-stress-controlled mechanisms, respectively. The ratio $c_o^X/c_u^X$ is varied and the resulting average Type X particle population in the plume is computed as follows

$$f_{\text{avg.}}^{X} = \frac{\int_{0}^{T} (c_{u}^{X} A_{u}(t) + c_{o}^{X} A_{o}(t)) dt}{\int_{0}^{T} A(t) dt}. \qquad (17)$$

## Data availability

The datasets generated during and/or analysed during the current study are available from the corresponding author on request. The proprietary data from ref. 27 (Supplementary Fig. 2), ref. 7 (Fig. 5), and ref. 32 (Fig. 5), used in this study, are accessible through these publications. Source data are provided in this paper. Further results, physical parameters' specification, parametric sensitivity study, as well as theoretical formulation of the vapour-transport model and derivation of the hydraulic model are provided in the Supplementary Information.

## Code availability

The Fortran90 routines for vapour flow model, Wolfram Mathematica scripts for the hydraulic model, as well as the Python scripts for statistical data processing and plotting, are available in a public Figshare repository⁵⁴.

## References

1. Hurford, T. A., Helfenstein, P., Hoppa, G. V., Greenberg, R. & Bills, B. G. Eruptions arising from tidally controlled periodic openings of rifts on Enceladus. *Nature* **447**, 292–294 (2007).

2. Hurford, T. A. et al. Geological implications of a physical libration on Enceladus. *Icarus* **203**, 541–552 (2009).

3. Hedman, M. M. et al. An observed correlation between plume activity and tidal stresses on Enceladus. *Nature* **500**, 182–184 (2013).

4. Nimmo, F., Porco, C. & Mitchell, C. Tidally modulated eruptions on Enceladus: Cassini ISS observations and models. *The Astronomical Journal* **148**, 46 (2014).

5. Porco, C. C., DiNino, D. & Nimmo, F. How the Geysers, Tidal Stresses, and Thermal Emission Across the South Polar Terrain of Enceladus Are Related. *The Astronomical Journal* **148**, 45 (2014).

6. Saur, J. et al. Evidence for temporal variability of Enceladus' gas jets: Modeling of Cassini observations. *Geophysical Research Letters* **35**, https://doi.org/10.1029/2008GL035811 (2008).

7. Hansen, C. J. et al. The composition and structure of Enceladus' plume from the complete set of Cassini UVIS occultation observations. *Icarus* **344**, 113461 (2020).

8. Teolis, B. et al. Enceladus plume structure and time variability: comparison of Cassini observations. *Astrobiology* **17**, 926–940 (2017).

9. Villanueva, G.L. et al. JWST molecular mapping and characterization of Enceladus' water plume feeding its torus. *Nature Astronomy*, https://doi.org/10.1038/s41550-023-02009-6 (2023).

10. Denny, K. E., Hedman, M. M., Bockelée-Morvan, D., Filacchione, G. & Capaccioni, F. Constraining time variations in enceladus's water-vapor plume with near-infrared spectra from cassini's visual and infrared mapping spectrometer. *The Planetary Science Journal* **5**, 144 (2024).

11. Běhounková, M. et al. Timing of water plume eruptions on Enceladus explained by interior viscosity structure. *Nature Geosci.* **8**, 601–604 (2015).

12. Beuthe, M., Rivoldini, A. & Trinh, A. Enceladus's and dione's floating ice shells supported by minimum stress isostasy. *Geophysical Research Letters* **43**, 10088–10096 (2016).

13. Čadek, O. et al. Enceladus's internal ocean and ice shell constrained from cassini gravity, shape, and libration data. *Geophysical Research Letters* **46**, 5653–5660 (2016).

14. Čadek, O. et al. Long-term stability of Enceladus' uneven ice shell. *Icarus* **319**, 476–484 (2019).

15. Hemingway, D. J. & Mittal, T. Enceladus's ice shell structure as a window on internal heat production. *Icarus* **332**, 111–131 (2019).

16. Běhounková, M., Souček, O., Hron, J. & Čadek, O. Plume activity and tidal deformation on Enceladus influenced by faults and variable ice shell thickness. *Astrobiology* **17**, 941–954 (2017).

17. Souček, O. et al. Tidal dissipation in Enceladus' uneven, fractured ice shell. *Icarus* **328**, 218–231 (2019).

18. Thomas, P. C. et al. Enceladus's measured physical libration requires a global subsurface ocean. *Icarus* **264**, 37–47 (2016).

19. Pleiner Sládková, K., Souček, O. & Běhounková, M. Enceladus' tiger stripes as frictional faults: Effect on stress and heat production. *Geophysical Research Letters* **48**, 2021–094849 (2021).

20. Berne, A. et al. Jet activity on Enceladus linked to tidally driven strike-slip motion along tiger stripes. *Nature Geoscience* **17**, 385–391 (2024).

21. Schmidt, J., Brilliantov, N., Spahn, F., Kempf, S. Slow dust in Enceladus' plume from condensation and wall collisions in tiger stripe fractures. *Nature* **451**, https://doi.org/10.1038/nature06491 (2008).

22. Ingersoll, A. P. & Pankine, A. A. Subsurface heat transfer on Enceladus: Conditions under which melting occurs. *Icarus* **206**, 594–607 (2010).

23. Kite, E. S. & Rubin, A. M. Sustained eruptions on Enceladus explained by turbulent dissipation in tiger stripes. *Proceedings of the National Academy of Sciences* **113**, 201520507 (2016).

24. Nakajima, M. & Ingersoll, A. P. Controlled boiling on Enceladus. 1. model of the vapor-driven jets. *Icarus* **272**, 309–318 (2016).

25. Spencer, J.R. et al. Plume Origins and Plumbing: From Ocean to Surface. In: Schenk, P.M., Clark, R.N., Howett, C.J.A., Verbiscer, A.J., Waite, J.H. (eds.) Enceladus and the Icy Moons of Saturn, p. 163, https://doi.org/10.2458/azu_uapress_9780816537075-ch008 (2018).

26. Elsworth, D. & Goodman, R. E. Characterization of rock fissure hydraulic conductivity using idealized wall roughness profiles. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts* **23**, 233–243 (1986).

27. Ingersoll, A. P., Ewald, S. P. & Trumbo, S. K. Time variability of the Enceladus plumes: orbital periods, decadal periods, and aperiodic change. *Icarus* **344**, 113345 (2020).

28. Oksanen, P. & Keinonen, J. The mechanism of friction of ice. *Wear* **78**, 315–324 (1982).

29. Yin, A., Zuza, A. V. & Pappalardo, R. T. Mechanics of evenly spaced strike-slip faults and its implications for the formation of tiger-stripe fractures on Saturn's moon Enceladus. *Icarus* **266**, 204–216 (2016).

30. Rossi, C., Cianfarra, P., Salvini, F., Bourgeois, O. & Tobie, G. Tec- tonics of Enceladus' south pole: Block rotation of the tiger stripes. *Journal of Geophysical Research: Planets* **125**, 2020–006471 (2020).

31. Petrovic, J. Review mechanical properties of ice and snow. *Journal of materials science* **38**, 1–6 (2003).
32. Postberg, F. et al. *Plume and Surface Composition of Enceladus*, pp. 129–162. University of Arizona Press, Tucson, http://www.jstor.org/stable/j.ctv65sw2b.16 (2018). Accessed 2023-06-23.

33. Postberg, F. et al. The E-ring in the vicinity of Enceladus: II. probing the moon's interior - the composition of E-ring particles. *Icarus* **193**, 438–454 (2008).

34. Postberg, F. et al. Sodium salts in E-ring ice grains from an ocean below the surface of Enceladus. *Nature* **459**, 1–4 (2009).

35. Postberg, F., Schmidt, J., Hillier, J., Kempf, S. & Srama, R. A salt-water reservoir as the source of a compositionally stratified plume on Enceladus. *Nature* **474**, 620–622 (2011).

36. Postberg, F. et al. Detection of phosphates originating from Enceladus's ocean. *Nature* **618**, 489–493 (2023).

37. Marusiak, A. G. et al. Exploration of icy ocean worlds using geophysical approaches. *The Planetary Science Journal* **2**, 150 (2021).

38. Ermakov, A. I. et al. A recipe for the geophysical exploration of Enceladus. *The Planetary Science Journal* **2**, 157 (2021).

39. Vance, S. et al. Distributed geophysical exploration of Enceladus and other ocean worlds. *Bulletin of the AAS* **53**, https://baas.aas.org/pub/2021n4i127 (2021).

40. Choblet, G. et al. Enceladus as a potential oasis for life: Science goals and investigations for future explorations. *Experimental Astronomy* **54**, 809–847 (2022).

41. Souček, O. et al. Radar attenuation in Enceladus’ ice shell: Obstacles and opportunities for constraining shell thickness, chemistry, and thermal structure. *Journal of Geophysical Research: Planets* **128**, 2022–007626 (2023).

42. MacKenzie, S. M. et al. The Enceladus orbilander mission concept: Balancing return and resources in the search for life. *The Planetary Science Journal* **2**, 77 (2021).

43. Savitzky, A. & Golay, M. J. E. Smoothing and differentiation of data by simplified least squares procedures. *Analytical Chemistry* **36**, 1627–1639 (1964).

44. Virtanen, P. et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods* **17**, 261–272 (2020).

45. Souček, O., Hron, J., Běhounková, M. & Čadek, O. Effect of the tiger stripes on the deformation of Saturn’s moon Enceladus. *Geophysical Research Letters* **43**, 7417–7423 (2016).

46. Dong, Y., Hill, T.W., Teolis, B.D., Magee, B.A., Waite, J.H. The water vapor plumes of Enceladus. *Journal of Geophysical Research: Space Physics* **116**, https://doi.org/10.1029/2011JA016693 (2011).

47. Hansen, C.J. et al. The composition and structure of the Enceladus plume. *Geophysical Research Letters* **38**, https://doi.org/10.1029/2011GL047415 (2011).

48. Olsson, R. & Barton, N. An improved model for hydromechanical coupling during shearing of rock joints. *International Journal of Rock Mechanics and Mining Sciences* **38**, 317–329 (2001).

49. Schulson, E.M., Duval, P. Creep and Fracture of Ice. Cambridge University Press, Cambridge (2009).

50. Maugis, D. Contact, Adhesion and Rupture of Elastic Solids, 1st edn. Springer Series in Solid-State Sciences Nr. 130. Springer, Heidelberg (2000).

51. Marriott, M. *Nalluri And Featherstone’s Civil Engineering Hydraulics: Essential Theory with Worked Examples*, 6th edn. Wiley-Blackwell, Chichester (2016).

52. Foreman-Mackey, D., Hogg, D. W., Lang, D. & Goodman, J. emcee: The MCMC Hammer. *Publications of the Astronomical Society of the Pacific* **125**, 306 (2013).

53. Goodman, J. & Weare, J. Ensemble samplers with affine invariance. *Communications in Applied Mathematics and Computational Science* **5**, 65–80 (2010).

54. Souček, O., Běhounková, M., Lanzendörfer, M., Tobie, G., Choblet, G. Figshare repository for NCOMMS-23-38372B: “Variations in plume activity reveal the dynamics of water-filled faults on Enceladus”. Figshare, https://doi.org/10.6084/m9 figshare. 26356342. (2024).

## Acknowledgements

This research was supported by the Czech Science Foundation through project No. 22-20388S (O.S., M.B.). This study acknowledges the financial support from the French Agence Nationale de Recherche, project ANR COLOSSE, grant No ANR-2020-CE49-0010 (G.C., G.T.).

## Author contributions

Authors' contributions: O.S. conceived the project and led the analysis. All authors (O.S., M.B., M.L., G.T., G.C.) contributed to the data analysis, the results discussion, and the writing of the manuscript. The development of numerical models was carried out by O.S. (Vapour-transport model, Hydraulic model, Aperture model), M.B. (Data and results processing and fitting) and M.L. (Hydraulic model).

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-024-51677-z.

**Correspondence** and requests for materials should be addressed to Ondřej Souček.

**Peer review information** *Nature Communications* thanks Andrew Ingersoll, and the other, anonymous, reviewers for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2024