# Identifying signatures of past and present cryovolcanism on Europa

Elodie Lesage, Samuel M Howell, Marc Neveu, Julia W Miller, Mariam

Naseem, Mohit Melwani Daswani, Justine Villette, Steven D Vance

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

► To cite this version:

Elodie Lesage, Samuel M Howell, Marc Neveu, Julia W Miller, Mariam Naseem, et al. Identifying signatures of past and present cryovolcanism on Europa. Nature Communications, 2025, 16 (1), pp.1886. 10.1038/s41467-025-57070-8 . hal-05230316

HAL Id: hal-05230316

https://hal.science/hal-05230316v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

# Identifying signatures of past and present cryovolcanism on Europa

Received: 30 July 2024

Accepted: 11 February 2025

Published online: 22 February 2025

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

Check for updates

Elodie Lesage Ⓔ<fnref n="1" />✉, Samuel M. Howell Ⓔ<fnref n="1" />, Marc Neveu Ⓔ²,<fnref n="3" />, Julia W. Miller<fnref n="1" />,⁴, Mariam Naseem Ⓔ⁵, Mohit Melwani Daswani Ⓔ<fnref n="1" />, Justine Villette Ⓔ⁶ & Steven D. Vance Ⓔ<fnref n="1" />

Europa, the most visibly active icy moon of Jupiter, is a prime target for the search for life in the outer solar system. Two spacecraft missions, Europa Clipper from the National Aeronautics and Space Administration (NASA) and the Jupiter Icy Moon Explorer (JUICE) from the European Space Agency (ESA), will observe its surface, probe its interior structure, and characterize the space environment starting in 2030. Occasional eruptions of water sourced from Europa's interior may provide a window on the interior conditions and habitability of the moon. Here, we investigate the storage and evolution of briny water in Europa's ice shell and propose a framework to interpret spectral, thermal, radar and gravity data collected by future missions. We show that it is possible to discriminate between water erupting from the deep ocean or from shallow liquid reservoirs using combined measurements of the material's salinity, surface temperature and ice shell thickness.

The eruption of briny solutions on the surfaces of icy ocean worlds has been proposed to explain the origin of vapor plumes, domes, dark deposits, and lobate, flow-like features<fnref n="1" />⁻⁹. Additional surface expressions such as chaos, lenticulae, double ridges, and central pits in impact craters are plausibly explained by liquid brine ("cryomagma") reservoirs in Europa's shallow ice shell¹⁰⁻¹⁵, despite challenges related to their formation and retention¹⁶⁻¹⁹. Subsurface liquid reservoirs may thus store water and non-water material and provide a mechanism for their transport through Europa's ice shell.

Cryomagma freezing is expected to pressurize reservoirs, potentially triggering eruptions²,⁷,²⁰,²¹. Potential erupted flows or plumes, if observed by upcoming missions, could provide new insights into the habitability and evolution of Europa's ice shell by revealing properties of transient and long-lived interior reservoirs of liquid water environments⁶⁻⁸,²¹⁻²⁷. The identification and characterization of erupted material and shallow subsurface liquid reservoirs would therefore be of great interest to understand Europa's habitability and, in the distant future, to access and analyze subsurface water²⁸⁻³³.

However, even if ongoing eruptions are observed directly, the composition of erupted material may not represent the composition of the reservoir at its time of formation, because melting and freezing of trapped brines are expected to affect their chemical composition through salt rejection, entrapment, precipitation, and dissolution²²,²⁵,³⁴,³⁵. Characterization of source reservoirs requires understanding how observed physico-chemical properties of surface features are affected by subsurface reservoir evolution processes, ambient ice temperature, as well as their depth within the ~ 20–40 km ice shell³⁶. In previous studies, Naseem et al.²² characterized the compositional evolution of brines of varying starting compositions relevant to Europa's ocean³⁷ during freezing of a cryoreservoir. Lesage et al.²¹ modeled reservoir rheology, pressurization and fracturing. Finally, Howell (2021)³⁶ presented a framework to model ice shells of time-variable thickness.

In this work, we bridge possible surface expressions of the erup-
ted material and the physical signatures of source reservoirs with their
physicochemical properties. We develop a novel program, CRYOLA-
VASAURUS, publicly available at https://github.com/ElodieLesage/
Cryolavasaurus³⁸. Within this program, we self-consistently simulate

<footnote n="1">¹Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA. ²Department of Astronomy, University of Maryland, College Park, MD, USA.</footnote>

<footnote n="3">³Planetary Environments Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD, USA. ⁴Department of Earth, Environmental, and Planetary Sciences, Brown University, Providence, RI, USA. ⁵Department of Geology, University of Maryland, College Park, MD, USA. ⁶Nantes Université, Univ Angers, Le Mans Université, CNRS, Laboratoire de Planétologie et Géosciences, LPG UMR 6112, 44000 Nantes, France. ✉ e-mail: elodie.lesage@jpl.nasa.gov</footnote>
the interplay of physical and chemical processes driving the freezing and eruption of shallow subsurface cryomagma reservoirs on Europa (see Methods, Section 1). We model an ice shell of time-variable thickness in which is embedded a shallow liquid reservoir of time-variable size and composition (see Methods, Section 2). The thermal, physical and compositional evolution of the ice shell and the reservoir, including convection in the ice shell and phase change at all interfaces, are computed simultaneously. We model heat transfer through a 1-dimensional (1D), spherically symmetric grid using an explicit finite-difference approach forward in time. Due to our 1D approximation, the reservoir is implemented as a layer of liquid water of imposed initial thickness and depth to the top of the layer. We track through time the ice shell temperature and thickness and the reservoir temperature, composition, and frozen fraction. We predict the characteristics of the erupted liquid at the time of each eruption, depending on the initial composition, depth, and thickness of the reservoir. Details on the numerical simulation are given in the Supplementary Model Description.

# Results

## Reservoir signatures and effect of salinity

Subsurface cryoreservoirs may leave several observable hints of their current or recent presence within the ice shell, including surface composition changes and thermal anomalies (Fig. 1a and b).

Continuous reservoir freezing induces overpressure, triggering repeated eruptions until no liquid remains or until pressure is relieved by viscous relaxation in the icy shell. Reservoirs 5 km thick (the thickest tested) remain active for up to 60 kyrs and trigger up to a few hundred eruptions, depending on their volume and depth. For a narrative explanation of the results, we chose to use a nominal reservoir that is 1 km thick and 1 km below the surface. For this reservoir, we obtain a total erupted volume of $4.5 \cdot 10^{-2} \text{ km}^3$ at the end of its 15 kyr active life (Fig. 1a), which represents an eruption of about 8% of the initial volume. The erupted volume is approximately the difference between the densities of ice and liquid water multiplied by the initial volume of the reservoir, which is expected for the assumptions we made: (i) the reservoir is disconnected from the ocean with no water refilling, (ii) the frozen ice layer has no porosity, keeping the water from infiltrating the ice shell, and (iii) the cryomagma is maintained at a pressure high enough that it does not volatilize before it reaches the surface.

For the nominal 1 km thick reservoir located 1 km below the surface, the surface thermal anomaly (Fig. 1b) resulting from the conduction and diffusion of the heat and latent heat from the reservoir peaks at +1.5 K after 20 kyr have elapsed since reservoir emplacement. The propagation of this heat in the ice shell also induces a local thinning below the reservoir (Fig. 1c). The thinning is insignificant for the shallow, nominal reservoir (<3 meters), but increases with reservoir

<figure><img src="image_5.png" /><figcaption>**Fig. 1** | Observable signatures and freezing characteristics for the 1 km thick reservoir located 1 km deep below Europa's surface. a Cumulative erupted volume, b surface temperature anomaly, c local ice shell thickness, d cryomagma</figcaption></figure>

freezing temperature, e freezing rate, and (f) volumic frozen fraction. Time (x-axis) scales differ between panels. The emplacement mechanism of the source reservoir is beyond the scope of this study.
<figure><img src="image_6.png" /><figcaption>**Fig. 2 | Evolution of the composition of a freezing 1 km thick, 1 km deep reservoir through time.** Here we show reservoirs injected from three different plausible ocean compositions corresponding to three different primordial Europa compositions²²,³⁷: a Cl chondrites, b CM chondrites, and c a mix of chondrites</figcaption></figure>

scaled to distance to the Sun with the addition of cometary material. Salts present in very low concentration are not represented on the plots (see Methods for more information on input salt molalities). 20 eruptions occur between each vertical gray line.

depth. In particular, thinning is strongest ~500 kyr after the onset of emplacement. Freezing is delayed due to the timescale of heat propagation by a factor of 30 relative to the freezing time of the reservoir. The thinning of the ice shell is also a long-lived signature, retaining 50% of its maximum value after >1 Myr.

In our simulations, freezing occurs slowly enough (Fig. 1d–f) that, combined with efficient mixing in the reservoir (see Methods, Section 3), it likely takes place at chemical equilibrium, meaning that the salts rejected in the liquid phase continue to participate in chemistry. This assumption is consistent with experimental data for terrestrial seawater and previous Europa literature³⁹. This remains true for all reservoir sizes (0.2–5 km) and depths (1–15 km) simulated.

Notably, despite the antifreeze properties of dissolved salts that decrease the solidus temperature, cryomagma composition has only a minor effect on freezing time, freezing rate, eruption interval and total erupted volume. The nominal reservoir freezes at 11.4 kyr if it is composed of pure water and <15 kyr if it is briny, with freezing rates of <0.6 cm/yr in all cases. As a result, the induced surface thermal anomaly, and ice shell and lithosphere thinning do not inform on source composition, but rather on reservoir physical properties. The effect of salinity is minor because only ≈ 1/15th of the energy exchanged with the reservoir arises from the transfer of thermal energy, the bulk arising from the latent heat released by liquid water upon freezing. These results confirm previous predictions from fixed eutectic freezing models²⁷,²¹.

## Erupted composition

Simulated eruption intervals and reservoir compositions are shown in Fig. 2 for three initial liquid compositions. Current-day oceanic compositions are predicted by rocky mantle devolatilization and water-rock interaction simulations of candidate Europa interiors; they depend on the assumed composition of material accreted during Europa's formation³⁷. We consider two early formation scenarios in Cl (Fig. 2a, "Cl") or CM chondrite-rich environments (Fig. 2b, "CM") and one formation scenario in mixed chondrites and water ice-depleted environment with subsequent addition of cometary material (Fig. 2c, "MC-Scale+comet"). Total dissolved solid concentrations for these formation scenarios are respectively 4.8, 10.5, and 45 g/kg, dominated by dissolved silica in the MC-Scale+Comet scenario. More information on the input compositions is given in Methods.

Importantly, despite the total salinity increasing during freezing, we find that relative proportions of the different salts in the erupted material are very similar to those in the initial reservoir, suggesting that determination of compositions of erupted material would provide meaningful insight into the composition of parent subsurface liquids. We also observe that for all three tested brine compositions, efficient salt rejection during freezing results in a gradual increase in solution salinity by an order of magnitude throughout the lifetime of the reservoir. This piece of information represents a first step in discriminating between eruptions sourced from the deep ocean or from evolved, local reservoirs. However, we do not address here how the composition and structure of the material are modified on the surface during solid-vapor phase changes or irradiation⁴⁰⁻⁴², which are discussed in the following.

## Influence of reservoir size and depth

Reservoir size and depth have a determining role in surface signatures (Fig. 3). The surface temperature anomaly above the reservoirs, assumed here to be of the same spatial extent as the reservoir itself, mainly depends on the reservoir depth. We obtain a maximum anomaly value of ~1.5 K for the shallowest reservoirs simulated (Fig. 3a). Variations in ice shell thickness are more pronounced with reservoir thickness and depth (Fig. 3b). This result is intuitive, as larger reservoirs release more energy stored in the liquid phase of water, while deeper reservoirs more rapidly and efficiently transport heat to the ice-ocean interface. Variations in ice shell thickness range from <1 m for reservoirs <1 km thick and 1 km deep, to 100 m for reservoirs >5 km thick and 5 km deep. Finally, the total amount of fluid that erupted from a reservoir depends on its initial volume and the density difference between the liquid and solid phases, but also on the temperature of the ice shell surrounding it. As the ice shell temperature increases, viscous relaxation becomes more efficient and prevents the reservoir internal pressure from building up, preventing the pressurization of reservoirs from achieving the critical stress required for eruption. Figure 3c shows the erupted volumes assuming a spherical reservoir. In our simulations, the deepest reservoirs (>10 km) never trigger eruptions for reservoir thicknesses >200 m.

Furthermore, we show that deeper reservoirs, located in warmer ice, freeze more slowly because the thermal gradients driving
<figure><img src="image_7.png" /><figcaption>**Fig. 3** | Effect of reservoir thickness and depth on various signatures. a Surface temperature anomaly, b ice shell thickness variation, and (c) total erupted volume. Blank cells in (c) indicate reservoirs that do not trigger eruptions.</figcaption></figure>

heat conduction are smaller, lengthening the interval between eruptions. As a result, the eruption frequency decreases with the depth of the reservoir. The local onset and/or enhancement of solid-state convection develops faster below deeper reservoirs, and convection cells take less time to reach the base of the ice shell. As a result of the exponential decrease in the viscosity of the ice with temperature, the active reservoir lifetime decreases with increasing reservoir depth. This might seem counter-intuitive as deeper, warmer reservoirs freeze at slower rates, but this finding is explained by the fact that eruptions cannot be triggered in ice that dissipates accumulating stress viscously at or faster than the freezing timescale²¹.

## Discussion

We find that the composition of the erupted material reflects proportionally the composition at the reservoir formation time. Thus, a comparison between the surface composition at a candidate cryovolcanic site and planned independent measurements of the salinity of Europa's ocean³,⁴,³,⁴⁴ will contribute evidence as to whether the parent reservoir liquid is sourced from the ocean or is instead due to melting of ice shell material. Eutectic mixtures melted in situ would refreeze near the eutectic, and the material stored in the reservoir and erupted at the surface would thus be relatively enriched in lower-eutectic-temperature salts, such as chlorides, compared to the oceanic compositions tested. For example, chloride enrichment has been reported
from Earth-based observations at higher concentrations in dark terrains of the trailing hemisphere⁴⁵.

Our findings also show that cryomagma salinity increases by an order of magnitude between reservoir emplacement and total freezing. Thus, regions identified as candidate cryovolcanic sites may show gradational spatial variations in surface composition associated with the eruptive sequence. These changes in surface properties may be detectable by the Europa Clipper mission's suite of optical remote sensing instruments comprising the Europa Imaging System (EIS)⁴⁶ and the Mapping Imaging Spectrometer for Europa (MISE)³⁵. The MISE camera is dedicated to the detection and mapping of non-ice materials (salts, acids, organics) and has a spatial resolution of up to 25 m/px⁴⁷.

We note that inference of reservoir composition from observations of erupted brines is limited by the following caveat: several processes can affect the composition and structure of ice grains formed from the erupted material. Detection of frozen brines and measurement of their composition is complicated by processes occurring during solidification, such as the vitrification of grains⁴¹. Enrichment in certain species and depletion in others also occur during phase changes, as demonstrated by laboratory work. Enrichment in NaCl preferentially to Na₂CO₃ has been observed under conditions relevant to Enceladus and Ceres⁴². Enrichment in Na₂SO₄ and MgCl₂ was also observed under Europa conditions⁴⁸. Finally, sublimation and irradiation of frozen brines also affect their composition on longer time scales, as shown by⁴⁰ for frozen chloride brines under Europa conditions. In general, a disconnect is often observed between source compositions and frozen brines, comprising an important and ongoing area of laboratory work and the development of compositional models⁴¹,⁴²,⁴⁹,⁵⁰. Potential opportunities to directly sample freshly erupted plume material³³,³⁵,⁴³ may be key in resolving the uncertainties introduced by exposure of cryomagma to space.

Although this study did not model the detection of any predicted signature presented here, these may be detectable by several instruments aboard the Europa Clipper spacecraft. Short-lived thermal anomalies as found here are plausibly detectable by the Europa Thermal Emission Imaging system (E-THEMIS), which has an expected temperature precision of 1.2 K³³,⁴³,⁵¹,⁵². Europa Clipper's radar instrument REASON may also be able to detect subsurface liquid brine or salt-rich layers⁵³, as well as more subtle changes in the depths at which thermally diagnostic interfaces occur (e.g. pore curing or liquid brines) using sparse echoes⁵⁴. Finally, density changes associated with the presence of a large-scale liquid reservoir and any additional local thinning of the ice shell might result in a gravity anomaly detectable by gravity and radio science⁴⁴.

In conclusion, the combined measurements of thermal, spectral, radar and gravity anomalies may enable the identification of current or recent (0 - ~10² kyr old) reservoir emplacement and freezing. This identification could constrain the solution space comprising the timing of emplacement, reservoir depth, and initial reservoir volume. Our simulations and results demonstrate how the upcoming Europa Clipper and JUICE missions will investigate problems central to the evolution and habitability of the Galilean moons through multidisciplinary measurements and their synergy³³,⁴³,⁵¹,⁵⁵⁻⁵⁸.

## Methods

### Model principle

CRYOLAVASAURUS simulates the thermal, mechanical, compositional, and temporal evolution of the ice shell and cryomagma reservoir by solving for the conservation of enthalpy. It uses finite differences in a one-dimensional (1D) spherically symmetric shell, propagated explicitly forward in time. The surface temperature, the thickness of the ice shell and its thermal structure, and the reservoir material phases and compositions are determined self-consistently by the simulation. Extensive documentation on the CRYOLAVASAURUS

framework and model can be found in Supplementary Model Description.

The thermal solver described in Supplementary Information, Section 1, includes:

* the exchange of energy due to temperature differences through conduction, diffusion, and convection;

* the exchange of latent heat during melting and freezing at the ice-reservoir and ice-ocean interfaces;

* the release or capture of energy due to changes in the specific enthalpy of formation of brines and hydrated salts;

* the addition of radiogenic heating from Europa's rocky interior at the seafloor;

* the rheology-dependent viscous dissipation of gravitational potential energy within the ice shell due to diurnal eccentricity tides;

* the conservative balance of energy associated with insolation, Europa's black body radiation, ice sublimation, and the escape of heat from Europa's icy surface.

Within this framework, we account for the following ways in which conservative energy transfer can affect the system:

* **Mechanical energy:** Pressure and volume change, allowing useful work to be done through deformation (e.g. fracturing) and mass transport (which fundamentally affects heat transport through advection).

* **Rate of mechanical change:** The mechanical properties of the material change (e.g. compressibility, elastic moduli, Lamé parameters, sound speed, rheology), affecting the future partitioning of energy storage and dissipation.

* **Thermal energy:** The temperature changes according to the thermal properties of the material and its proximity to phase boundaries in both the state of matter and the composition.

* **Composition:** Mineralogy may change, both in abundance and in the creation and destruction of some minerals, according to the specific enthalpies of formation of the possible materials, their physical and thermal properties, and changes in stability of materials as pressure and temperature evolve.

* **Phase:** Materials can change their state of matter according to the pressure and temperature evolution of their surroundings with respect to their phase space.

* **Rates of thermal and compositional change:** The thermal properties of the material change, controlling how much the next Joule of energy will alter the temperature, phase, and composition.

## Boundary and Initial conditions

In this section, we summarize the conditions imposed at key interfaces within the simulation.

**Ice shell properties.** The temperature at the ice-ocean interface is set as the melting temperature of the oceanic material. At the surface (ice-space interface), the temperature is self-consistently calculated at each time step using a heat flux balance. This balance accounts for radiogenic heat flux from Europa's interior, solar insolation heat flux, sublimation heat flux, and radiative heat flux. Details on this calculation are given in the Supplementary Model Description (Section 1). The thickness of the ice shell is not imposed, but is determined self-consistently by the simulation. Models are initiated with a conductive temperature profile near the equilibrium thickness and then allowed to evolve thermally for 5 Myr prior to reservoir emplacement so that a steady-state ice shell structure is ensured.

**Reservoir emplacement.** Whether it occurs rapidly by injection into the ocean or more slowly by in situ melting, reservoir emplacement must happen much faster than the timescale for the reservoir to freeze¹⁶,¹⁷,⁵⁹. Here, we consider that the reservoir is emplaced
<table><caption>**Table 1 | Molalities (mol element per kg of H₂O) of the three input compositions used in this study**</caption><thead><tr><th></th><th>CI</th><th>CM</th><th>MC-Scale+Comet</th></tr></thead><tbody><tr><td>Ca</td><td>9.39 × 10<sup>-3</sup></td><td>1.86 × 10<sup>-3</sup></td><td>3.02 × 10<sup>-6</sup></td></tr><tr><td>Mg</td><td>5.16 × 10<sup>-3</sup></td><td>2.12 × 10<sup>-2</sup></td><td>0</td></tr><tr><td>Na</td><td>3.05 × 10<sup>-1</sup></td><td>9.19 × 10<sup>-2</sup></td><td>5.59 × 10<sup>-1</sup></td></tr><tr><td>K</td><td>2.63 × 10<sup>-3</sup></td><td>5.84 × 10<sup>-3</sup></td><td>1.62 × 10<sup>-2</sup></td></tr><tr><td>Cl</td><td>1.52 × 10<sup>-2</sup></td><td>2.00 × 10<sup>-2</sup></td><td>4.07 × 10<sup>-2</sup></td></tr><tr><td>S</td><td>1.46 × 10<sup>-2</sup></td><td>4.47 × 10<sup>-2</sup></td><td>0</td></tr><tr><td>C</td><td>2.63 × 10<sup>-2</sup></td><td>4.24 × 10<sup>-2</sup></td><td>1.51 × 10<sup>-2</sup></td></tr><tr><td>Si</td><td>0</td><td>2.86 × 10<sup>-5</sup></td><td>4.89 × 10<sup>-1</sup></td></tr></tbody></table>

Compositional data from<sup>22,37</sup>.

instantaneously within the ice shell once the ice shell temperature reaches equilibrium to avoid the practical implications of simultaneously resolving numerous questions regarding reservoir formation and emplacement, and to maintain the focus of this study on reservoir evolution and detection. Because of our 1D approximation for the model grid, the liquid reservoir does not technically own a horizontal extent, but is rather described by a global layer of water of imposed initial thickness. However, we use a spherical description to better simulate energy transfer in and around the reservoir. We interpolate heat loss at the surface and base of the reservoir around the shell perimeter, conservatively scaled by the far-field ice temperature. This is a good approximation with only small uncertainties introduced through *de minimis* heat transport radially along the thin, nearly isothermal boundary normal to the horizontal plane.

**Initial cryomagma composition.** We use three plausible compositions for Europa's ocean, as well as a reference model with pure water. Melwani Daswani et al.<sup>37</sup> modeled the water-rock interaction for a young Europa to infer plausible present-day compositions of the ocean. The resulting ocean composition depends on the environment during Europa's accretion, e.g. the composition of chondrites and the bulk amount of available water. As detailed in Melwani Daswani et al.<sup>37</sup>, the oceanic composition “CI” refers to a formation scenario in an environment rich in CI chondrites, and the composition “CM” in an environment rich in CM chondrites. The “MC-Scale+comet” composition was obtained using Monte Carlo simulations where the probability of accreting different types of carbonaceous chondrites (CI, CM, CV, CK, CR, and CO) was weighted by their formation distance from Jupiters. Cometary material was then added to compensate for insufficient volatiles produced by this accretion scenario<sup>60</sup>, using comet 67P as a reference. These fundamentally different formation models explain the order-of-magnitude differences in elemental concentrations such as Si and Ca between scenarios. Table 1 summarizes the total element molalities for the three compositions used in the present study.

**Ice shell and reservoir thermal evolution**

The thicknesses of the reservoir and ice shell evolve over time to account for phase-related enthalpy exchange at all interfaces. That is, the top and bottom reservoir boundaries and the ice-ocean interface move over time as a result of freezing and melting.

At the surface, the temperature evolves self-consistently by balancing incoming diurnally averaged equatorial insolation, $\dot{q}_{sol}$, and geologic heat flow through the surface, $\dot{q}_{geo}$, with heat losses due to blackbody radiation, $\dot{q}_{bb}$, and sublimation to maintain a rarefied atmosphere, $\dot{q}_{sub}$:

$$ \dot{q}_{sol} + \dot{q}_{geo} = \dot{q}_{bb} + \dot{q}_{sub}. \qquad (1) $$

See the Supplementary Model Description (Section 2) for further information.

We assume that the cryomagma is convecting during freezing, because the estimated Rayleigh number in the reservoir exceeds values that typically indicate convection ($Ra > 10^6 - 10^8$). For physical parameters on the order of those used in this study ($\rho = 1000 \text{ kg m}^{-1}$ the liquid density, $g = 1.35 \text{ m s}^{-2}$ Europa's gravity, $\alpha = 10^{-3} \text{ K}^{-1}$ the water thermal expansivity, $cp = 4.2 \text{ kJ kg}^{-1}$ the liquid water specific heat capacity, $\eta = 10^{-3} \text{ Pa s}$ the liquid water viscosity, $k = 0.5 \text{ Wm}^{-1} \text{K}^{-1}$ the liquid water thermal conductivity; see<sup>61</sup>), and with a realistic maximum temperature gradient of 1 K across the reservoir and a characteristic reservoir vertical extent of 2 km, we get $Ra \approx 10^{19}$. Vigorous convection implies that heat is efficiently transported in the liquid phase, and thus we consider it isothermal.

Finally, we model phase changes at the reservoir interfaces and the ice-ocean interface. We do so by tracking the energy flux across these boundaries (see Section 4 below). In the reservoir, we use the freezing paths from Naseem et al.<sup>22</sup> (Section 4) to estimate the amount of cryomagma that freezes given the energy exchanged at the interfaces at each time step. We update the temperature solution of the boundary elements consistently with the compositional evolution of reservoir and icy material, and freezing, using the above-referenced simulation to partition changes in energy. We track the interface locations and update all physical parameters. The process we use to calculate the effective thermal conductivity at interface nodes is described in the Supplementary Model Description, Section 1.

**Cryomagma chemical evolution during freezing**

Naseem et al.<sup>22</sup> modeled the freezing of inclusions of ocean material in the ice shell using the oceanic compositions predicted by Melwani Daswani et al.<sup>37</sup> as inputs. In Naseem et al.<sup>22</sup>, compositions are treated as a function of decreasing temperature<sup>22</sup> and obtained using PHREEQC v3<sup>62</sup> and the thermodynamic data converted by Toner et al.<sup>63</sup> from previous work<sup>34,64</sup>. The model in Naseem et al.<sup>22</sup> calculates the chemical sequence of forming species with decreasing temperatures, while ignoring the time evolution of the system temperatures. The link to the data spreadsheets used in this study, as well as information on their labeling, are given in the Supplementary Model Description, Section 3.

The freezing rates we obtain are under 1 cm/yr, indicating a slow freezing scenario in which salts are efficiently rejected from the formed ice and concentrate in the liquid phase, where the salinity increases over time. Because convection is efficient (Section 3), precipitated salts are assumed to remain suspended in the liquid phase and still participate in chemical reactions (equilibrium crystallization). This end-member case was simulated in<sup>22</sup> and we use the corresponding data as input for our model.

We do not account for the effect of phase change-induced pressure changes on compositions, which is expected to be negligible<sup>22</sup>. Liquid density is an output of PHREEQC and is used to calculate the overpressure in the reservoir. Density change due to compression is calculated and used to determine if an eruption is triggered at each time step (Section 5).

To bridge the PHREEQC outputs of composition as a function of decreasing temperature with the computed reservoir thermal evolution, we track energy removed from the reservoir in both models. Our simulation uses time steps, and the PHREEQC data use temperature steps, but we can reconcile these by expressing both models in terms of energy steps. The initial total amount of energy stored in the reservoir can be calculated using the liquid temperature $T_{liq}$:

$$ E_{init} = \rho_{liq} cp T_{liq}. \qquad (2) $$

We then track how much energy is removed from the reservoir at each time step of the simulation and subtract it from $E_{init}$. In parallel, we calculate how much energy is left in the reservoir for each PHREEQC
data point using the changes in temperature and amount of liquid and solid materials in the PHREEQC output files. Using interpolations between the PHREEQC data points, we retrieve the composition associated with the energy state of the reservoir at each time step. More details about this process are available in the Supplementary Model Description, Section 4.

## Reservoir pressurization and eruption

We further develop the framework of Lesage et al.²¹ to track reservoir pressurization and eruption. Reservoir pressurization arises from the increase in water volume after freezing, with volume expansion opposed by the incompressibility of the melt. When the stress on the reservoir wall exceeds a threshold value, determined as function of the ice tensile strength, a fracture propagates to the surface and an eruption is triggered.

We self-consistently evaluate the accumulation of elastic stresses due to pressurization of the compressible medium, the viscoelastic relation of those stresses with time according to the rheology, and the ability of the system to trigger a new eruption mechanically. Details on the model implementation of these processes are included in the Supplementary Model Description, Section 5.

## Data availability

The data generated in this study have been deposited in the public GitHub repository https://github.com/ElodieLesage/Cryolavasaurus/tree/main/Onput and permanently archived at https://doi.org/10.5281/zenodo.14768339. The compositional data from Naseem et al.²², as well as the input files used to run the CRYOLAVASAURUS simulation, are available at https://github.com/ElodieLesage/Cryolavasaurus/tree/main/Onput and permanently archived at https://doi.org/10.5281/zenodo.14768339. These files contain the input parameters used to obtain all the results presented in this study. All input and output data is fully available at the links above and downloadable without restrictions.

## Code availability

The compositional data showed in Fig. 2 was obtained using the PHREEQC routine available at https://github.com/MarcNeveu/frezchem and following the methods in Naseem et al.²². The CRYOLAVASAURUS software used to obtain the results of this study is open-source and fully available at https://github.com/ElodieLesage/Cryolavasaurus, and archived at https://doi.org/10.5281/zenodo.14768339.

## References

1. Fagents, S. A. et al. Cryomagmatic mechanisms for the formation of Rhadamanthys linea, triple band margins, and other low-albedo features on Europa. *Icarus* **144**, 54–88 (2000).

2. Fagents, S. A. Considerations for effusive cryovolcanism on Europa: the post-Galileo perspective. *J. Geophys. Res.* **108**, 5139 (2003).

3. Porco, C. C. et al. Cassini observes the active south pole of Enceladus. *Science* **311**, 1393–1401 (2006).

4. Sparks, W. B. et al. Probing for evidence of plumes on Europa with HST/STIS. *Astrophys. J.* **829**, 121 (2016).

5. Quick, L. C., Glaze, L. S. & Baloga, S. M. Cryovolcanic emplacement of domes on Europa. *Icarus* **284**, 477–488 (2017).

6. Quick, L. C. & Hedman, M. M. Characterizing deposits emplaced by cryovolcanic plumes on Europa. *Icarus* **343**, 113667 (2020).

7. Lesage, E., Massol, H. & Schmidt, F. Cryomagma ascent on Europa. *Icarus* **335**, 113369 (2020).

8. Lesage, E., Schmidt, F., Andrieu, F. & Massol, H. Constraints on effusive cryovolcanic eruptions on Europa using topography obtained from Galileo images. *Icarus* **361**, 114373 (2021).

9. Quick, L. C. et al. Cryolava dome growth resulting from active eruptions on Jupiter's moon Europa. *Icarus* **387**, 115185 (2022).

10. Collins, G. & Nimmo, F. Chaotic Terrain on Europa. https://europa.nasa.gov/resources/120/chaotic-terrain-on-europa/ (2009).

11. Schmidt, B., Blankenship, D., Patterson, G. & Schenk, P. Active formation of 'chaos terrain' over shallow subsurface water on Europa. *Nature* **479**, 502–505 (2011).

12. Michaut, C. & Manga, M. Domes, pits, and small chaos on Europa produced by water sills. *J. Geophys. Res. Planets* **119**, 550–573 (2014).

13. Steinbrügge, G. et al. Brine migration and impact-induced cryovolcanism on Europa. *Geophys. Res. Lett.* **47**, e2020GL090797 (2020).

14. Culberg, R., Schroeder, D. M. & Steinbrügge, G. Double ridge formation over shallow water sills on Jupiter's moon Europa. *Nat. Commun.* **13**, 1–10 (2022).

15. Chivers, C., Buffo, J. & Schmidt, B. Stable brine layers beneath Europa's chaos. *Planetary Sci. J.* **4**, 159 (2023).

16. Craft, K. L., Patterson, G. W., Lowell, R. P. & Germanovich, L. Fracturing and flow: Investigations on the formation of shallow water sills on Europa. *Icarus* **274**, 297–313 (2016).

17. Vilella, K., Choblet, G., Tsao, W.-E. & Deschamps, F. Tidally heated convection and the occurrence of melting in icy satellites: application to Europa. *J. Geophys. Res. Planets* **125**, e2019JE006248 (2020).

18. Carnahan, E., Vance, S. D., Cox, R. & Hesse, M. A. Surface-to-ocean exchange by the sinking of impact generated melt chambers on Europa. *Geophys. Res. Lett* **49**, e2022GL100287 (2022).

19. Hesse, M. A., Jordan, J. S., Vance, S. D. & Oza, A. V. Downward oxidant transport through Europa's ice shell by density-driven brine percolation. *Geophys. Res. Lett.* **49**, e2021GL095416 (2022).

20. Neveu, M., Desch, S., Shock, E. & Glein, C. Prerequisites for explosive cryovolcanism on dwarf planet-class Kuiper belt objects. *Icarus* **246**, 48–64 (2015).

21. Lesage, E., Massol, H., Howell, S. M. & Schmidt, F. Simulation of freezing cryomagma reservoirs in viscoelastic ice shells. *Planetary Sci. J.* **3**, 170 (2022).

22. Naseem, M. et al. Salt distribution from freezing intrusions in ice sheets on ocean worlds. *Planetary Sci. J.* **4**, 181 (2023).

23. Quick, L. C., Barnouin, O. S., Prockter, L. M. & Patterson, G. W. Constraints on the detection of cryovolcanic plumes on Europa. *Planetary Space Sci.* **86**, 1–9 (2013).

24. Hand, K., Sotin, C., Hayes, A. & Coustenis, A. On the habitability and future exploration of ocean worlds. *Space Sci. Rev.* **216**, 1–24 (2020).

25. Chivers, C., Buffo, J. & Schmidt, B. Thermal and chemical evolution of small, shallow water bodies in Europa's ice shell. *J. Geophys. Res. Planets* **126**, e2020JE006692 (2021).

26. Vance, S. D., Journaux, B., Hesse, M. & Steinbrügge, G. The salty secrets of icy ocean worlds. *J. Geophys. Res. Planets* **126**, e2020JE006736 (2021).

27. Wolfenbarger, N. S., Fox-Powell, M. G., Buffo, J., Soderlund, K. M. & Blankenship, D. D. Brine volume fraction as a habitability metric for Europa's ice shell. *Geophys. Res. Lett.* **49**, e2022GL100586 (2022).

28. Howell, S. M. & Leonard, E. J. In Ocean worlds: Interior processes and physical environments. In *Handbook of Space Resources* (eds. Badescu, V., Zacny, K., Bar-Cohen, Y.) 873–906 (Springer, 2023).

29. Howell, S. et al. Ocean Worlds exploration and the search for life. *Bull. Am. Astron. Soc.* **53**, 191 (2021).

30. Hockman, B. J. et al. PRIME: Probe using Radioisotopes for Icy Moons Exploration-a comprehensive cryobot architecture for accessing Europa's ocean. *Astrobiol. Sci. Conf.* 502–04 (2022).

31. Lawrence, J. D. et al. Subsurface science and search for life in ocean worlds. *Planetary Sci. J.* **4**, 22 (2023).

32. Winkenstern, J. & Saur, J. Detectability of local water reservoirs in Europa's surface layer under consideration of coupled induction. *J. Geophys. Res. Planets* **128**, e2023JE007992 (2023).
33. Pappalardo, R. T. et al. Science overview of the Europa Clipper mission. *Space Sci. Rev.* **220**, 1–58 (2024).

34. Marion, G. M., Kargel, J. S., Catling, D. C. & Jakubowski, S. D. Effects of pressure on aqueous chemical equilibria at subzero temperatures with applications to Europa. *Geochimica. Cosmochimica. Acta* **69**, 259–274 (2005).

35. Becker, T., Zolotov, M., Gudipati, M. et al. Exploring the composition of Europa with the upcoming Europa Clipper mission. *Space Sci. Rev.* **220**, 49 (2024).

36. Howell, S. M. The likely thickness of Europa's icy shell. *Planetary Sci. J.* **2**, 129 (2021).

37. Melwani Daswani, M., Vance, S. D., Mayne, M. J. & Glein, C. R. A metamorphic origin for Europa's ocean. *Geophys. Res. Lett.* **48**, e2021GL094143 (2021).

38. Lesage, E., Howell, S. M. & Miller, J. W. CRYOLAVASURUS. https://github.com/ElodieLesage/Cryolavasaurus/releases/tag/v1 (2025).

39. Zolotov, M. Y. & Shock, E. L. Composition and stability of salts on the surface of Europa and their oceanic origin. *J. Geophys. Res. Planets* **106**, 32815–32827 (2001).

40. Thomas, E. C., Hodyss, R., Vu, T. H., Johnson, P. V. & Choukroun, M. Composition and evolution of frozen chloride brines under the surface conditions of Europa. *ACS Earth Space Chem.* **1**, 14–23 (2017).

41. Vu, T. H., Choukroun, M., Hodyss, R. & Johnson, P. V. Probing Europa's subsurface ocean composition from surface salt minerals using in-situ techniques. *Icarus* **349**, 1137–1146 (2020).

42. Naseem, M. & Neveu, M. *Changes in Salt Composition of Ocean Materials During Exposure to Vacuum On Icy Ocean Worlds*. https://agu.confex.com/agu/abscicon24/meetingapp.cgi/Paper/1497131 (2024).

43. Roberts, J. H. et al. Exploring the interior of Europa with the Europa clipper. *Space Sci. Rev.* **219**, 46 (2023).

44. Mazarico, E. et al. The Europa Clipper gravity and radio science investigation. *Space Sci. Rev.* **219**, 30 (2023).

45. Ligier, N., Poulet, F., Carter, J., Brunetto, R. & Gourgeot, F. VLT/SINFONI observations of Europa: new insights into the surface composition. *Astron. J.* **151**, 163 (2016).

46. Turtle, E. et al. The Europa imaging system (EIS) investigation. *Space Sci. Rev.* **220**, 1–68 (2024).

47. Blaney, D. L. et al. The mapping imaging spectrometer for Europa (MISE). *Space Sci. Rev.* **220**, 80 (2024).

48. Vu, T. H., Hodyss, R., Choukroun, M. & Johnson, P. V. Chemistry of frozen sodium–magnesium–sulfate–chloride brines: Implications for surface expression of Europa’s ocean composition. *Astrophys. J. Lett.* **816**, L26 (2016).

49. Johnson, P. V., Hodyss, R., Vu, T. H. & Choukroun, M. Insights into Europa's ocean composition derived from its surface expression. *Icarus* **321**, 857–865 (2019).

50. Howell, S. M. et al. *Enceladus: An Enthalpic Jet Propulsion Laboratory*. AGU23 https://ui.adsabs.harvard.edu/abs/2023AGUFM.P43E3321H/abstract (2023).

51. Daubar, I. et al. Planned geological investigations of the Europa clipper mission. *Space Sci. Rev.* **220**, 18 (2024).

52. Christensen, P. R. et al. The Europa thermal emission imaging system (E-THEMIS) investigation for the Europa clipper mission. *Space Sci. Rev.* **220**, 1–65 (2024).

53. Blankenship, D. D. et al. Radar for Europa assessment and sounding: ocean to near-surface (REASON). *Space Sci. Rev.* **220**, 51 (2024).

54. Schroeder, D. et al. Constraining the thickness of the conductive portion Europa's ice shell using sparse radar echoes. *Geophys. Res. Lett.* **51**, e2024GL110635 (2024).

55. Howell, S. M. & Pappalardo, R. T. NASA's Europa clipper—a mission to a potentially habitable ocean world. *Nat. Commun.* **11**, 1311 (2020).

56. Vance, S. D. et al. Investigating Europa's habitability with the Europa clipper. *Space Sci. Rev.* **219**, 81 (2023).

57. Grasset, O. et al. JUPITER ICy moons explorer (JUICE): an ESA mission to orbit ganymede and to characterise the Jupiter system. *Planetary Space Sci.* **78**, 1–21 (2013).

58. Howell, S. M. et al. Jupiter's ocean worlds: dynamic ices and the search for life. *Ice Solar Syst.* **2024**, 283–314 (2024).

59. Kalousová, K., Souček, O., Tobie, G., Choblet, G. & Čadek, O. Ice melting and downward transport of meltwater by two-phase flow in Europa's ice shell. *J. Geophys. Res. Planets* **119**, 532–549 (2014).

60. Mousis, O. et al. Early stages of Galilean moon formation in a water-depleted environment. *Astrophys. J. Lett.* **944**, L37 (2023).

61. Wagner, W. & Kretzschmar, H.-J. IAPWS Industrial Formulation 1997 For The Thermodynamic Properties of Water and Steam. https://iapws.org/relguide/IF97-Rev.html (2008).

62. Parkhurst, D. L. & Appelo, C. Description of input and examples for PHREEQC version 3–a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. *US Geol. Surv. Tech. Methods* **6**, 497 (2013).

63. Toner, J. D. & Sletten, R. S. The formation of Ca-Cl-rich groundwaters in the dry valleys of Antarctica: field measurements and modeling of reactive transport. *Geochimica. Cosmochimica. Acta* **110**, 84–105 (2013).

64. Marion, G., Catling, D., Crowley, J. & Kargel, J. Modeling hot spring chemistries with applications to martian silica formation. *Icarus* **212**, 629–642 (2011).

## Acknowledgements

Portions of this research were carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration (NASA). This work was primarily supported by NASA's Solar System Workings program grants #80NM0018F0612 (E.L., S.M.H. and J.W.M.) and #80NSSC20K0139 (E.L., S.M.H., M.Ne., M.Na., M.M.D., S.D.V.). Parts of this work were also supported by NASA under award #80GSFC24M0006 (M.Ne., M.Na., S.D.V.) and through the Precursor Science Investigations for Europa program grant #80NM0018F0612 (S.D.V).

## Author contributions

E.L., S.M.H. and J.W.M. developed the CRYOLAVASURUS numerical model and ran the simulations. M.Ne. and M.Na. developed the PHREEQC routine and provided the composition data. E.L., S.M.H., M.Ne., M.Na., J.W.M., M.M.D., J.V., and S.D.V. discussed the findings of the study and contributed to the final manuscript.

## 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-025-57070-8.

**Correspondence** and requests for materials should be addressed to Elodie Lesage.

**Peer review information** *Nature Communications* thanks Sarah Fagents 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 4.0 International License, which permits use, sharing, adaptation, 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 changes were made. 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/4.0/.

© The Author(s) 2025