# Simulation of Freezing Cryomagma Reservoirs in Viscoelastic Ice Shells

Elodie Lesage, Hélène Massol, Samuel M. Howell, Frédéric Schmidt

► To cite this version:

Elodie Lesage, Hélène Massol, Samuel M. Howell, Frédéric Schmidt. Simulation of Freezing Cryomagma Reservoirs in Viscoelastic Ice Shells. The Planetary Science Journal, 2022, 3, 10.3847/PSJ/ac75bf. insu-03745278

HAL Id: insu-03745278

https://insu.hal.science/insu-03745278

Submitted on 4 Aug 2022

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Simulation of Freezing Cryomagma Reservoirs in Viscoelastic Ice Shells

Elodie Lesage<fnref n="1" /> ID, Hélène Massol<fnref n="2" /> ID, Samuel M. Howell<fnref n="1" /> ID, and Frédéric Schmidt<fnref n="2" />,<fnref n="3" /> ID

Received 2022 March 24; revised 2022 May 23; accepted 2022 June 1; published 2022 July 21

## Abstract

Brine reservoirs in Europa's icy crust, if they exist, could represent the most accessible liquid water bodies in the outer solar system. Any potential subsurface reservoirs are key for the exploration of ocean worlds and the search for habitability and life beyond Earth. Previous studies have demonstrated that freezing cryoreservoirs might trigger eruptions due to the pressurization associated with volume change as liquid water expands to become water ice, but those studies did not take into account the deformation of the reservoir wall. Viscoelastic deformation of the wall can act to accommodate the growing overpressure and prevent eruptions. Here we utilize a semianalytical numerical approach to calculate the pressure-induced deformation of a freezing cryoreservoir on Europa as a function of the far-field ice temperature. We find that shallow reservoirs located in relatively cold ice deform slightly and can be pressurized by freezing, but that deeper reservoirs located in warmer ice deform more easily and in some cases cannot reach the threshold overpressure required for spontaneous eruption. We identify a transition depth around 4–8 km depending on the reservoir radius, assuming a conductive ice shell structure consistent with current best estimates. Further, we find that shallow lens-shaped reservoirs may store a large volume of cryomagma and can be easily pressurized by freezing; they represent the best candidates for the mobilization of cryovolcanic material at Europa's surface.

*Unified Astronomy Thesaurus concepts: Europa (2189); Volcanism (2174); Natural satellite surfaces (2208)*

## 1. Introduction

Jupiter's icy moon Europa possesses one of the youngest surfaces of the solar system, with an average surface age estimated around 40–90 Myr (Zahnle et al. 2003; Bierhaus et al. 2009). The great variety of geological features observed on Europa's surface demonstrates widespread ongoing or geologically recent activity (Greeley et al. 1998; Doggett et al. 2009; Kattenhorn & Prockter 2014; Howell & Pappalardo 2018). Additionally, it is very likely that Europa hosts a global liquid water ocean that is believed to be in contact with its rocky mantle (Khurana et al. 1998; Kuskov & Kronrod 2005) and situated beneath an ~20–50 km thick icy crust (Howell 2021). Based on surface observations and subsurface modeling, several studies have also hypothesized the potential presence of liquid reservoirs in the ice crust (e.g., Fagents 2003; Schmidt et al. 2011; Dombard et al. 2013; Craft et al. 2016; Manga & Michaut 2017). These reservoirs, if they exist, could represent the most accessible liquid water bodies in the outer solar system and are key for the exploration of ocean worlds and the search for life beyond Earth.

Europa is the target of the upcoming National Aeronautics and Space Administration (NASA) Europa Clipper mission (Howell & Pappalardo 2020) and one target of the upcoming European Space Agency (ESA) JUICE mission (Grasset et al. 2013). In anticipation of the arrival of these spacecraft, it is important to address the question of whether or not liquid water could reach Europa's surface and to constrain the sources of any potential observed flows.

Some geological features observed in images collected by the Galileo spacecraft, such as domes (Quick et al. 2017) and smooth plains (Fagents 2003), have been interpreted as flows of water-based material. However, the origin and plausibility of ocean-surface water transport and/or in situ melting remain a topic of heavy debate. Previous studies have demonstrated that eruption of liquid water rising directly from the ocean through cracks is unlikely because of the unrealistically high pressure this mechanism requires: oceanic pressure of the order of $10^5–10^6$ Pa, which may be implausible (Manga & Wang 2007). Even with significant overpressure from global ice shell thickness changes, deep viscoelastic relaxation and significant gravity are likely to prevent any throughgoing crack formation on Europa (Rudolph et al. 2022).

Hydrologically isolated reservoirs of briny liquid water, canonically referred to as cryomagma, are one alternative eruption source (Fagents 2003; Neveu et al. 2015; Lesage et al. 2020). Although the formation and stability of these reservoirs are challenged by the difficulty of melting more than a few percent of ice (Vilella et al. 2020) and a potentially short percolation timescale to the ocean (Kalousová et al. 2014), the eruption of fluid from these reservoirs could be initiated by the freezing of only a few percent of their volume (Neveu et al. 2015; Lesage et al. 2020). This work chiefly addresses whether isolated bodies of liquid water contained within the ice shell are capable of spontaneous eruption due to overpressure induced by freezing, while considering the viscoelastic deformation of their surroundings. Thus, we remain agnostic to the origin of any specific isolated liquid water body and instead focus on their eruptability.

A few authors have studied the pressurization of such cryoreservoirs. Fagents (2003) first proposed that the freezing of a few percent of cryomagma could be sufficient to generate a pressure great enough to overcome the negative buoyancy of

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹NASA Jet Propulsion Laboratory, California Institute of Technology, CA, USA; elodie.lesage@jpl.nasa.gov</footnote>

<footnote n="2">²Université Paris-Saclay, CNRS, GEOPS, 91405, Orsay, France</footnote>

<footnote n="3">³Institut Universitaire de France (IUF), France</footnote>
liquid water in the ice crust. Neveu et al. (2015) also predicted that the freezing of a small fraction of cryomagma in a reservoir can generate pressures up to several MPa, which is enough to fracture the ice shell and trigger an eruption. Finally, Lesage et al. (2020) proposed a model to simulate the freezing of a briny cryoreservoir and predicted the critical freezing time required to trigger an eruption and the erupted volume of cryolava. Nevertheless, these studies do not take into account the ice viscoelastic rheology, which can allow the reservoirs to accommodate the freezing overpressure. This mechanism may increase the critical freezing time or even prevent the eruption in some cases.

On Earth, the deformation of pressurized magma chambers is well described compared to their planetary counterparts and has been related to surface deflection of the ground before and after an eruption (Folch et al. 2000; Segall 2016) or the formation of volcanic activity-related structures such as calderas (Jellinek & DePaolo 2003). McTigue (1987) first calculated the deformation of a spherical cavity in an elastic medium due to the increase of its internal pressure. Nevertheless, as demonstrated by Hodge (1974), a shell of partially molten material develops around magma reservoirs because of the thermal transfer between magma and the surrounding rock, locally decreasing the parent rock viscosity. Thus, to more realistically model the evolution of pressurized magma reservoirs, Dragoni & Magnanensi (1989) proposed a new analytical model to calculate the deformation and stress fields around pressurized spherical reservoirs surrounded by a viscoelastic shell. However, the model of Dragoni & Magnanensi (1989) was developed in the context of terrestrial magma chambers with a predictable inner pressure source function, which is consistent with sustained magma chamber replenishment from mantle sources.

Sustained sources of pressure, like those invoked for silicate magmatism on Earth, may not represent potential water volcanism on Europa. There, reservoirs might be isolated from the internal ocean, with their inner pressure increasing with time due to freezing. For this case, the pressure evolution and reservoir deformation are intrinsically linked and cannot be decoupled, and it is necessary to adapt the deformation model of Dragoni & Magnanensi (1989) to a different pressure source that takes into account the reservoir deformation. Gregg et al. (2013) proposed a two-dimensional finite element model of magma chamber deformation adapted to increasing inner pressure. To address the first-order controls on the eruption of deformable wall magma chambers, we simplify the approach of Gregg et al. (2013) to a pseudoanalytical, one-dimensional, spherically symmetric finite element model.

Lesage et al. (2020) speculated that, to first order, the ability to erupt was intuitively dependent on the Maxwell time of the viscoelastic ice and the critical freezing time in the case of a nondeformable reservoir wall. Here we improve upon our first cryomagma freezing model, taking into account the viscoelastic deformation of ice. To do so, we use the formulation proposed by Dragoni & Magnanensi (1989) and adapt it to Europa's ice shell and putative cryoreservoirs. We thus refine our predictions regarding the feasibility for cryomagma reservoirs on Europa to erupt and calculate the overpressure needed to overcome the reservoir deformation, as well as its associated critical freezing time and erupted volume. We predict the conditions for which reservoir deformation inhibits the eruption by accommodating the freezing overpressure. We then compare

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic model of a cryomagma reservoir surrounded by viscoelastic and elastic ice. The cryomagma reservoir has a radius $R_1$, and the viscoelastic shell has a radius $R_2$. Here $K_1$ and $K_2$ are the bulk moduli of, respectively, the viscoelastic and elastic ice, and $G_1$ and $G_2$ are their respective shear moduli. Here $r$ is the radial coordinate, and $\eta$ is the viscosity of the viscoelastic ice (where $R_1 < r < R_2$). In the reservoir, the overpressure $\Delta P$ comes from cryomagma freezing. Adapted from Dragoni & Magnanensi (1989).</figcaption></figure>

these results to the first-order estimates of Lesage et al. (2020). Finally, we extend the reservoir freezing model to lens-shaped reservoirs with more realistic, horizontally elongated geometries.

## 2. Model

### 2.1. Principle

We begin by assuming that a liquid water reservoir is emplaced in Europa's ice shell. Whether emplaced by melting or by fracturing and ascent of oceanic water, it acts as a heat source that warms and decreases the viscosity of the surrounding ice, promoting viscoelastic behavior (further details are given in Section 2.3). For this reason, we use the approach of Dragoni & Magnanensi (1989), which is summarized in Figure 1. We define $R_1$ to be the reservoir radius and $R_2$ the viscoelastic shell radius (see Figure 1).

Each of the physical parameters used below are summarized in Table 1. The Cryoreservoir Freezing-Induced Eruption Simulation (CryoFRIES) code developed and used in this study is available at https://github.com/ElodieLesage/CryoFRIES.git.

### 2.2. Critical Freezing Time of a Nondeformable Reservoir

As described in Lesage et al. (2020), we make the assumption that cryomagma reservoirs are pressurized by freezing as the liquid phase transitions to the less dense solid phase. In this model, cryomagma gradually freezes from the reservoir wall as heat is conducted into the surrounding ice. We made, in the first model, the approximation of a nondeformable reservoir wall that does not accommodate the overpressure. We call this the “nondeformable wall approximation” throughout the paper. Under this assumption, the overpressure $\Delta P$ in the reservoir is

$$ \Delta P = -\frac{1}{\chi} \ln \left( \frac{1 - n(\rho_l / \rho_s)}{1 - n} \right), \quad (1) $$

where $\chi$ is the liquid water compressibility, $n$ is the cryomagma frozen fraction, and $\rho_l$ and $\rho_s$ are, respectively, the densities of the cryomagma and associated ice. The frozen fraction $n$ is
<table><caption>**Table 1**

Parameters Used in This Study</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>c<sub>p</sub></td><td>Specific heat capacity of ice</td><td>2.3 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>E</td><td>Young modulus</td><td>3 &middot; 10<sup>9</sup> Pa</td></tr><tr><td>G</td><td>Shear modulus</td><td>3 &middot; 10<sup>9</sup> Pa</td></tr><tr><td>H</td><td>Reservoir depth</td><td>1–10 km</td></tr><tr><td>k</td><td>Thermal conductivity of ice</td><td>2.3 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>K</td><td>Bulk modulus</td><td>E/(3(1-2ν))</td></tr><tr><td>R<sub>1</sub></td><td>Reservoir radius</td><td>10 m–5 km for spherical reservoirs, up to 100 km for lens-shaped reservoirs</td></tr><tr><td>R<sub>2</sub></td><td>Viscoelastic shell radius</td><td>See Figures 1 and 4 and Equation (9)</td></tr><tr><td>T<sub>cold</sub></td><td>Temperature in the ice crust far from any reservoir</td><td>Linear gradient from 110 (surface) to 250 K (bottom of conductive layer, 10 km deep)</td></tr><tr><td>T<sub>m</sub></td><td>Cryomagma melting temperature</td><td>273 K for pure water, 268 K for briny cryomagma (Kargel 1991)</td></tr><tr><td>T(R<sub>2</sub>)</td><td>Temperature at which ice becomes viscoelastic around a reservoir</td><td>See Section 2.3 and Equation (9)</td></tr><tr><td>V<sub>e</sub></td><td>Erupted volume with nondeformable wall approximation</td><td>See Equation (6)</td></tr><tr><td>V<sub>ev</sub></td><td>Erupted volume taking into account the elastic and viscoelastic deformation</td><td>Estimated using the iterative scheme described in Section 2.5</td></tr><tr><td>ΔP<sub>c</sub></td><td>Critical overpressure necessary to fracture the reservoir with the nondeformable wall approximation</td><td>2(σ<sub>c</sub> + ρ<sub>s</sub>gH), see Equation (4)</td></tr><tr><td>η</td><td>Viscosity</td><td>See Equation (7) (applies to viscoelastic region only, i.e., R<sub>1</sub> &lt; r &lt; R<sub>2</sub>)</td></tr><tr><td>κ</td><td>Thermal diffusivity of ice</td><td>k/(ρ<sub>s</sub>c<sub>p</sub>)</td></tr><tr><td>ν</td><td>Poisson ratio</td><td>0.325</td></tr><tr><td>ρ<sub>l</sub></td><td>Density of liquid cryomagma</td><td>1000 kg m<sup>-3</sup> for pure water, 1180 kg m<sup>-3</sup> for briny cryomagma</td></tr><tr><td>ρ<sub>s</sub></td><td>Density of ice</td><td>920 kg m<sup>-3</sup> for pure water ice, 1130 kg m<sup>-3</sup> for briny ice</td></tr><tr><td>σ<sub>c</sub></td><td>Tensile strength of ice</td><td>Linear gradient from 1.7 &middot; 10<sup>6</sup> (surface) to 1 &middot; 10<sup>6</sup> (10 km deep; Litwin et al. 2012)</td></tr><tr><td>τ</td><td>Viscoelastic deformation time</td><td>See Dragoni &amp; Magnanensi (1989), Equation (9)</td></tr><tr><td>τ<sub>M</sub></td><td>Maxwell time</td><td>η/E</td></tr><tr><td>τ<sub>c</sub></td><td>Critical freezing time necessary to fracture the reservoir with the nondeformable wall approximation</td><td>See Equation (5)</td></tr><tr><td>τ<sub>cv</sub></td><td>Critical freezing time necessary to fracture the reservoir taking into account the ice elastic and viscoelastic deformation</td><td>Estimated using the iterative scheme described in Section 2.5</td></tr><tr><td>χ</td><td>Liquid water compressibility</td><td>5 &middot; 10<sup>-10</sup> Pa<sup>-1</sup></td></tr></tbody></table>

calculated as a function of time by solving the Stefan problem (see Lesage et al. 2020, Appendix A, for further details),

$$n = 1 - \left( 1 - \frac{S}{R_1} \right)^3. \qquad (2)$$

Here $R_1$ is the reservoir radius, and $S$ is the cryomagma frozen layer thickness as a function of time,

$$S = 2\lambda\sqrt{\kappa t}, \qquad (3)$$

with $\lambda$ being a constant given by solving the Stefan problem and $\kappa$ the ice thermal diffusivity. Combining Equations (1)–(3), we can calculate the evolution of the overpressure in a freezing reservoir with time, noting that we are still under the nondeformable wall approximation. Figure 2(a) shows the pressure evolution in a 500 m radius reservoir located 2 km beneath the ice shell surface.

We solve the Stefan problem for the mean temperature around a reservoir (calculated by averaging the temperature along the reservoir wall), rather than a discrete approach that allows the wall temperature to vary spatially. This makes the calculation much faster and is a valid approximation, as the overpressure in the reservoir is calculated using the mean thickness of the frozen cryomagma layer, regardless of its spatial distribution. Given that the heat transfer equation is

linear as a function of the temperature, this approximation is not expected to impact the freezing rate.

The critical overpressure necessary to break a reservoir wall is (Lesage et al. 2020, Equation (12))

$$\Delta P_c = 2(\sigma_c + \rho_s g H), \qquad (4)$$

where $\sigma_c$ is the ice tensile strength (Litwin et al. 2012), and $H$ is the reservoir depth. By combining Equations (1)–(3) and setting $\Delta P = \Delta P_c$, we can calculate the critical freezing time $\tau_c$ necessary to trigger an eruption with the nondeformable wall approximation

$$\tau_c = \left( \frac{S_c}{2\lambda\sqrt{\kappa}} \right), \qquad (5)$$

where $\lambda$ is a constant obtained by solving the Stefan problem (see Lesage et al. 2020 for details) and $\kappa$ is the ice thermal diffusivity. Then, during an eruptive event, the cryolava erupted volume $V_e$ corresponds to the cryomagma volume increase in the reservoir due to phase change prior to the eruption. This volume is thus proportional to the solid and liquid densities,

$$V_e = V_t n_c \frac{\rho_l}{\rho_s}, \qquad (6)$$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Pressure in the reservoir due to cryomagma freezing as modeled by Lesage et al. (2020) under the nondeformable wall approximation. Here $\Delta P_c$ is the critical overpressure necessary to trigger an eruption, and $\tau_c$ is the critical freezing time necessary to reach $\Delta P_c$. (b) Linear pressure source function used in this model.</figcaption></figure>

where $V_i$ is the volume of liquid cryomagma before freezing, and $n_c$ is the critical frozen fraction of cryomagma.

## 2.3. Ice Rheology

In this study, we assume a linear thermal gradient in the outermost 10 km of Europa's icy crust, which corresponds to the portion of the ice shell where the dominant mechanism of global heat transfer is thermal conduction. This imposed conductive layer ranges from 250 K at 10 km depth to 110 K at the surface. This corresponds to the current best estimate for the thickness (10.4 km) and basal temperature (~252 K) spanned by the conductive layer of the ice shell (Howell 2021). This initial ice shell temperature profile does not take into account the presence of a warm reservoir and is called $T_{cold}$ hereafter. More locally, the ice temperature around a freezing reservoir is calculated by solving the Stefan problem (see Lesage et al. 2020 for details).

Because of the thermal gradients both in the ice crust and around reservoirs, we expect some parts of the ice shell to behave as a viscoelastic, rather than fully elastic, material. To calculate the viscosity $\eta$ of the viscoelastic ice as a function of its temperature $T$, we use the Newtonian law of Hillier & Squyres (1991):

$$ \eta = 10^{14} \exp \{25.2(273/T - 1)\}. \quad (7) $$

At high stresses and strain rates, like those associated with tectonic deformation, the ice will exhibit non-Newtonian behavior that can locally increase the viscosity. At lower stresses and strain rates (such as those associated with convection), Newtonian diffusion creep dominates (Goldsby & Kohlstedt 2001; Howell & Pappalardo 2018). Thus, we capture best the first-order temperature dependence of viscosity, and while incorporating non-Newtonian behavior was beyond the scope of this study's iterative approach, incorporating a more complex rheology may act to mitigate reservoir deformation by increasing the icy viscosity.

As reservoirs may extend vertically across several hundred or even thousand meters, we must calculate the mean viscosity along the wall of a given reservoir. To take into account the maximum deformation possible of each reservoir and thus not overestimate their eruption likelihood, we use a harmonic mean

for the viscosity, preferentially weighting the lower viscosity values. The viscosity harmonic mean is defined as

$$ \bar{\eta} = \left( \frac{1}{n} \sum_{i=1}^{n} \frac{1}{\eta(T_i)} \right)^{-1}, \quad (8) $$

where $T_1$ and $T_n$ are the temperatures, respectively, at the top and bottom of the reservoir. We choose a small temperature step, i.e., $T_{i+1} - T_i \ll T_n$.

To understand the ice rheology around a reservoir at a given depth, it is necessary to compare two timescales: (i) the Maxwell relaxation time of the ice, $\tau_M = \eta/E$, with $E$ being the ice Young modulus and $\eta$ the temperature-dependent viscosity; and (ii) the critical freezing time $\tau_c$ of the reservoir, which determines the timescale over which the stresses generated by the freezing overpressure are applied to the reservoir wall. A first-order approximation of $\tau_c$ is calculated by making the assumption of a nondeformable reservoir wall (see Lesage et al. 2020, Equation (17)).

Figure 3 shows a comparison between $\tau_M$ and $\tau_c$ for a 500 m radius reservoir. One can see that in the first 6 km of ice extending from the surface, $\tau_M \gg \tau_c$ (light blue region of Figure 3), which means the viscous component of the strain in response to the freezing stress is negligible in this region. Above 6 km, the ice behaves mostly as an elastic material during the freezing, and for more simplicity, we refer to this region as the elastic ice hereafter. Deeper than 6 km, $\tau_M \ll \tau_c$, which means the ice behaves as a viscoelastic material during freezing. We call it viscoelastic ice in the rest of this paper. One should note that the elastic-viscoelastic transition depth is a function of the critical freezing time and thus depends on the reservoir radius. It is also inversely proportional to the reservoir depth and varies between approximately 7 and 5 km for reservoirs of radius ranging from 100 m to 1 km, respectively.

On a more local scale, the calculation of the Maxwell time around a reservoir allows derivation of the thickness of the viscoelastic shell surrounding a reservoir located in elastic ice, $R_2 - R_1$. This can be accomplished by finding the location $r = R_2$, where ice is cold enough to behave as a purely elastic material, i.e., where $\tau_M$ becomes greater than $\tau_c$. To determine $R_2$, we first find the temperature $T(R_2)$ at which $\tau_M = \tau_c$ using Equation (7) to represent the ice viscosity. The radius $R_2$ can
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Black line: Maxwell time $\tau_M$ of the ice, far from any reservoir. Blue line: critical freezing time $\tau_c$ of a 500 m radius reservoir. When $\tau_M \gg \tau_c$ (light blue region), the ice can be approximated as an elastic medium, and the viscous component is negligible. When $\tau_M \lesssim \tau_c$ (white region), the ice behaves as a viscoelastic material. The transition takes place at $H \simeq 6$ km for a 500 m radius reservoir at a temperature $T \simeq 194$ K.</figcaption></figure>

then be determined by finding the location where the temperature around the reservoir equals $T(R_2)$. An expression for $R_2$ can be obtained by manipulating the solution to the Stefan problem for the temperature around the reservoir (see Lesage et al. 2020, Equation (31) and Figure 12) and substituting $T = T(R_2)$,

$$R_2 = R_1 - 2\sqrt{\kappa\tau_c}\text{erf}^{-1}\left(\frac{T(R_2) - T_{\text{cold}}(H)}{T_m - T_{\text{cold}}(H)}(1 + \text{erf}(\lambda)) - 1\right), \quad (9)$$

where $\kappa$ is the ice thermal diffusivity, $T_{\text{cold}}(H)$ is the far-field temperature, and $T_m$ is the cryomagma melting temperature.

As an example, the thickness of the viscoelastic shell around a 500 m radius reservoir is shown in Figure 4 as a function of the reservoir depth. One can see that the viscoelastic shell is relatively thin, with a maximum thickness of $1.3 \cdot R_1$. Again, the elastic–viscoelastic transition appears in Figure 4 at around 6 km; deeper than 6 km, the thickness of the viscoelastic shell diverges, as there is no elastic region any more. Thus, for reservoirs located deeper than the transition depth, we adapt the model by considering a fully viscoelastic material.

## 2.4. Deformable Wall

To predict the response of a cryomagma reservoir to the increase of its inner pressure, we adapt the formulation of Dragoni & Magnanensi (1989). That study calculated an analytical solution for the one radial ($\sigma_{rr}(r, t)$) and two tangential ($\sigma_{\theta\theta}(r, t)$ and $\sigma_{\phi\phi}(r, t)$) stress fields around a pressurized magma reservoir, as well as the wall displacement ($u(r, t)$). Each of these quantities is a function of the radial position, $r$, and time, $t$. We follow the approach of Dragoni & Magnanensi (1989), whereby we pose these equations for the purely elastic case and then use the correspondence principle (Biot 1954) to generalize the equations to a viscoelastic material. The final expressions we used for $u(r, t)$, $\sigma_{rr}(r, t)$, $\sigma_{\theta\theta}(r, t)$, and $\sigma_{\phi\phi}(r, t)$ are the following:

$$u(r, t) = \frac{A_1(r)}{D_0}f_1(t) + \frac{A_2(r)}{D_0}f_2(t), \quad (10)$$

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Dimensionless thickness of the viscoelastic shell around a 500 m radius reservoir at the end of the freezing ($t = \tau_c$ calculated with the nondeformable wall approximation) as a function of its depth in the ice crust.</figcaption></figure>

$$\sigma_{rr}(r, t) = \frac{B_1(r)}{D_0}f_1(t) + \frac{B_2(r)}{D_0}f_2(t), \quad (11)$$

$$\sigma_{\theta\theta}(r, t) = \sigma_{\phi\phi}(r, t) = \frac{C_1(r)}{D_0}f_1(t) + \frac{C_2(r)}{D_0}f_2(t). \quad (12)$$

The coefficients $A_1, A_2, B_1, B_2, C_1, C_2$, and $D_0$ and the two functions $f_1(t)$ and $f_2(t)$ used in Equations (10)–(12) are provided in Appendix 2 of Dragoni & Magnanensi (1989). One should note that here we used the more common planetary notation for shear modulus, $G$, referred to as $\mu$ in Dragoni & Magnanensi (1989). We use Equation (7) to calculate these coefficients.

Equations (10)–(12) are based on a trapezoidal overpressure function, written as

$$\Delta P(t) = \Delta P_c \left\{ \frac{t}{t_1} [1 - H(t - t_1)] + H(t - t_1) - H(t - t_2) + \frac{t_3 - t_1}{t_3 - t_2} [H(t - t_2) - H(t - t_3)] \right\}, \quad (13)$$

where $t=0$ is the beginning of the pressure increase, $t=t_1$ is the beginning of the pressure plateau, $t=t_2$ is the beginning of the pressure decrease, and $t=t_3$ is the time at which the excess pressure has been fully released. In our case, we assume that the pressure increases linearly during freezing and decreases abruptly at $t = \tau_c$ when the eruption begins. We set $t_1 = t_2 = t_3 = \tau_c$ to obtain a linear pressure source function. By doing so, we adapt the trapezoidal pressure function of Dragoni & Magnanensi (1989) to the simpler pressure function shown in Figure 2(b).

## 2.5. Pressure Source Function: Iterative Numerical Model

As described above, we improve upon Lesage et al. (2020) by investigating and quantifying the pressure changes induced by (i) the cryomagma freezing and (ii) the reservoir deformation. Because these behaviors are coupled, it is not possible to
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Numerical procedure used to calculate the critical pressure and freezing time in the case of a deformable wall.</figcaption></figure>

derive an analytical form of the unique realistic critical freezing time to reach the critical overpressure, $\Delta P_c$, and trigger an eruption. To solve this problem and find the critical freezing time, hereafter called $\tau_{cv}$ (where the index $v$ stands for viscoelastic), we use the following iterative process, illustrated in Figure 5.

1. We calculate the critical freezing time, $\tau_c$, of a reservoir and the associated overpressure, $\Delta P_c$, by making the nondeformable wall assumption using the model from Lesage et al. (2020).

2. We calculate the corresponding deformation $u(R_1, \tau_c)$ from Equations (10)–(12) at time $t = \tau_c$ and the induced pressure drop, $\delta P$, using the definition of the water compressibility, $\delta P = -\frac{1}{\chi} \frac{V_f}{V_i}$, with $\chi$ being the water compressibility, $V_i$ the reservoir volume before deformation, and $V_f$ its volume after deformation.

3. We calculate the critical freezing time necessary to reach an overpressure $\Delta P_c + \delta P$ and simulate new reservoir deformation after this longer freezing time, $\tau_c + \delta\tau$. This deformation induces a new pressure drop of $\delta P_2$.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Pressure drop in the reservoir at each iterative step of the algorithm and the associated critical freezing time. The algorithm converges to a realistic critical freezing time $\tau_{cv}$ in 17 iterations. The results were obtained with a shallow reservoir (see reservoir A in Table 2).</figcaption></figure>

4. We verify whether $\delta P - \delta P_2 < \epsilon$, where $\epsilon$ is an imposed convergence criterion; we use $\epsilon = \Delta P_c \cdot 10^{-5}$. If the model converges, the algorithm has found the realistic critical freezing time $\tau_{cv} = \tau_c + \delta\tau$ to reach the critical overpressure, $\Delta P_c$. If not, the algorithm iterates beginning with step 3, replacing $\delta P$ with $\delta P_2$ and $\delta\tau$ with $\delta\tau_2$, and then verifies the convergence, and so on.

We iterate until the program converges to a realistic critical freezing time, $\tau_{cv}$, that takes into account the feedback between pressurization and deformation. As an example, Figure 6 shows the pressure drop, $\delta P$, caused by the reservoir deformation at each iteration for a 500 m radius reservoir located 2 km deep in the ice crust. We can see that the iterative process converges relatively quickly to a realistic pressure function in this case; 17 iterations were sufficient to reach the convergence criterion.

In some cases, the deformation is great enough to completely accommodate the freezing overpressure at each iteration, and the algorithm cannot find a realistic $\tau_{cv}$. In this case, the model diverges. This result indicates that the modeled reservoir cannot be pressurized by freezing and thus cannot trigger eruptions. We demonstrate in the Results section (see Section 3.2) that reservoirs for which the algorithm diverges largely correspond to cases where $\tau_c > \tau_M$. Hence, perhaps intuitively, the Maxwell time of the ice can be used to discriminate reservoirs that can or cannot trigger eruptions.

## 2.6. Application to Lens-shaped Reservoirs

The model described above is adapted to spherical reservoirs. Nevertheless, as is the case for terrestrial magma chambers (Sigurdsson et al. 1999), cryomagma reservoirs on Europa may have an elongated, lens, or sill shape (Michaut & Manga 2014; Craft et al. 2016; Manga & Michaut 2017). As Europa's conductive upper icy shell is assumed to be 10 km thick in our model (Howell 2021), the assumption of spherical reservoirs limits their radius to 5 km and thus their volume to $<10^{12}$ m³. Moreover, for reservoirs with a significant vertical extent, the wall viscosity can be several orders of magnitude lower at the reservoir bottom than at the top, so large spherical reservoirs deform mostly from their bottom. This behavior may be very different than that of elongated, shallow reservoirs embedded in high-viscosity ice. For these reasons, we adapted our model to simulate the freezing of elongated lens-shaped reservoirs, which

are modeled as ellipsoids with $a$ and $b$ as two equal horizontal semiaxes and $c$ as the vertical semiaxis. For simplicity, we denote the geometry of ellipsoidal reservoirs by calculating a semiaxis aspect ratio, $F$, where $a=b=F \cdot c$.

To robustly adapt our model to ellipsoidal reservoirs, we primarily modify the temperature field around the reservoirs to match their new vertical extent. The Stefan problem is solved in Cartesian coordinates in our model, which remains a realistic approximation because the freezing layer is thin compared to the radius of the cavity.

The deformation model of Dragoni & Magnanensi (1989) is formulated in spherical geometry and thus cannot be applied to ellipsoid reservoirs. We demonstrate in Section 3.1 that the deformation of reservoirs located entirely within elastic ice is negligible, and deformation of reservoirs stored in viscoelastic ice prevents their eruption. Hence, we chose to calculate the critical freezing time and erupted volume of ellipsoid reservoirs without taking into account their viscoelastic deformation, and we use the comparison between their critical freezing time $\tau_c$ and Maxwell time $\tau_M$ to determine whether an eruption is possible or not, the efficacy of which is demonstrated in Section 2.5.

## 3. Results

We first show in Section 3.1 the evolution of the deformation and stress fields around a 500 m radius reservoir. The initial iteration considers the linear pressure source function obtained with the nondeformable wall approximation. Each successive iteration considers coupled changes in the stress and displacement fields. We find that the algorithm never converges for the viscoelastic reservoirs located in warm, low-viscosity ice and justifies the need to study the nondeformable case as a comparison between the elastic and viscoelastic reservoir environments. In Section 3.2, we use the iterative model to calculate the realistic pressure source function, taking into account the feedback between pressure and deformation. We then generalize these results to reservoirs of varying radius and depth and predict the critical freezing time and erupted volume for reservoirs able to trigger eruptions. These results are more realistic than the one predicted by Lesage et al. (2020) and, more importantly, allow us to identify which reservoirs can or cannot trigger eruptions.

### 3.1. Deformation and Stress Fields around a Freezing Reservoir

The results presented in this section are obtained for a spherical 500 m radius reservoir filled with pure liquid water (see Section 3.2 for results with varying radius, depth, and composition). The elastic-viscoelastic transition for these reservoirs is located approximately 6 km deep within the ice shell for a 10 km thick conductive layer. We model these two different cases as follows and summarized in Table 2.

1. The elastic case is modeled with a 2 km deep reservoir located in elastic ice and surrounded by a thin viscoelastic shell, the thickness of which is calculated to be around 10 m (reservoir A).

2. The viscoelastic case is modeled as an 8 km deep reservoir located in fully viscoelastic ice (reservoir B). In this case, as there is no elastic ice far from the reservoir, we set $R_2$ equal to half of the thickness of Europa's viscoelastic ice shell ($R_2 = 2$ km here), and we only model the region spanning $R_1 < r < R_2$.
Reservoir A: shallow elastic ice

<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a) Deformation and (b) radial and (c) tangential stress fields around a spherical 500 m radius reservoir located 2 km deep within the ice shell (reservoir A; see Table 2). Reservoir A is surrounded by a 10 m viscoelastic ice shell, itself embedded in elastic ice. Results are plotted at time $\tau_c$, which is the critical freezing time required to reach an eruption under the nondeformable wall approximation calculated at the first iteration of our algorithm. As $\sigma_{\theta\theta}$ reaches the critical value $\Delta P/2$ (see Equation (4)), an eruption is possible in the case of a nondeformable wall, but in reality, a longer critical freezing time $\tau_{cv}$ will be necessary to reach the critical overpressure once we take into account the wall deformation.</figcaption></figure>

<table><caption>**Table 2**
Parameters Used to Model Spherical Reservoirs A (Shallow) and B (Deep)</caption><thead><tr><td>Reservoir</td><td>Radius R<sub>1</sub></td><td>Depth</td><td>Viscoelastic Shell Radius R<sub>2</sub></td><td>Temperature T<sub>cold</sub></td><td>Critical Freezing Time τ<sub>c</sub></td></tr></thead><tbody><tr><td>A</td><td>500 m</td><td>2 km</td><td>1.025 · R<sub>1</sub></td><td>142 K</td><td>3 yr</td></tr><tr><td>B</td><td>500 m</td><td>8 km</td><td>4 · R<sub>1</sub></td><td>226 K</td><td>110 yr</td></tr></tbody></table>

The deformation and stress fields around the reservoirs at $t = \tau_c$ in the elastic and viscoelastic cases are given in Figures 7 and 8, respectively. These results are obtained by applying the linear pressure source function obtained with the nondeformable wall approximation (see Figure 2). We show the solutions at $t = \tau_c$, which is the critical freezing time under the nondeformable wall approximation calculated at the first iteration of our algorithm.

As shown in Figure 7(a), reservoir A (elastic case) is lightly deformed due to freezing; after the critical freezing time $\tau_c$, the reservoir radius increases by 30 cm. Figures 7(b) and (c) show that the critical radial and tangential stresses are reached at $t = \tau_c$. As $\sigma_{\theta\theta}$ reaches the critical value $\Delta P_c/2$, the eruption may be triggered. However, these results do not take into account the pressure drop caused by the wall deformation, so the critical stress required to trigger the eruption will not actually be reached by $t = \tau_c$ but only after a longer critical

freezing time $\tau_{cv}$. In the particular example of reservoir A, we obtained $\tau_{cv} = 12$ yr, constituting a fourfold increase in critical freezing time and thus a significant increase in the intereruption timescale.

Figure 8(a) shows that reservoir B (viscoelastic case) is much more deformed due to the pressure increase; the reservoir radius increases by 200 m over the critical freezing time $\tau_c$. In this case, the viscous rheology of the ice dominates, and tangential stresses are in an extensional state around the reservoir (Figure 8(c)), which is typical of the viscoelastic rheology (e.g., Dragoni & Magnanensi 1989). The tangential stress, $\sigma_{\theta\theta}$, never reaches the critical value, $\Delta P_c/2$, necessary to fracture the reservoir. This indicates that the overpressure is accommodated by the reservoir deformation, so reservoir B is unable to trigger an eruption. As the eruption cannot take place, the overpressure continues to increase as long as some liquid remains in the reservoir. We interrupt the algorithm after a few
Reservoir B: deep viscoelastic ice

<figure><img src="image_11.png" /><figcaption>**Figure 8.** (a) Deformation and (b) radial and (c) tangential stress fields around a spherical 500 m radius reservoir located 8 km deep in the ice shell (reservoir B; see Table 2). Reservoir B is surrounded by viscoelastic ice only. Results are plotted at time $t = \tau_c$, which is the critical freezing time required to reach an eruption under the nondeformable wall approximation calculated at the first iteration of our algorithm. As $\sigma_{\theta\theta}$ never reaches the critical value $\Delta P/2$ (see Equation (4)), eruptions cannot be triggered ($\tau_{cv}$ tends toward $\infty$).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Critical freezing time before eruption of a 500 m radius reservoir as a function of its depth in the nondeformable wall approximation (gray line) and with the wall deformation (black line). Two cryomagma compositions are modeled: (a) pure liquid water and (b) a briny cryomagma ($81 \text{ wt\% H}_2\text{O} + 16 \text{ wt\% MgSO}_4 + 3 \text{ wt\% Na}_2\text{SO}_4$).</figcaption></figure>

iterations following model divergence to save time and conserve sanity.

## 3.2. Critical Freezing Time and Erupted Volume

The results presented in this section are calculated with the iterative process introduced in Section 2.5 to calculate realistic

critical freezing times $\tau_{cv}$ that take into account the feedback between pressure and deformation. Figure 9 shows the critical freezing times $\tau_c$ and how they compare with the realistic critical freezing times $\tau_{cv}$ calculated with the nondeformable wall approximation. We show $\tau_c$ and $\tau_{cv}$ as functions of the reservoir depth and fixed $R = 500$ m for this example. Although
the critical freezing time is increased because of the wall
deformation, it still remains in the same order of magnitude as
$\tau_c$ for a given reservoir depth and radius. The same effect is
observed for the (first eruption event) erupted volume $V_e$ in our
calculations that is of the same order of magnitude as $V_{ev}$. This
demonstrates that the nondeformable wall approximation is
acceptable for reservoirs located in elastic ice. This first
important result may simplify future studies and further
modeling of cryomagmatic reservoirs.

Figures 10 and 11 summarize our results in terms of
critical freezing timescale and erupted volume for a large
range of reservoir radii (from 10 m to a few kilometers) and
depth (from 1 to 10 km below the surface, assuming the
conductive layer is 10 km thick). Figure 10 shows the critical
freezing time necessary to trigger an eruption for reservoirs
filled with pure liquid water (panels (a), (c), and (e)) or briny
cryomagma sourced from the ocean (81 wt% H₂O + 16 wt%
MgSO₄ + 3 wt% Na₂SO₄; Kargel 1991; panels (b), (d), and
(f)). This ocean composition was obtained with models of
brine evolution in CI chondrite bodies and represents an
upper bound for oceanic salt concentration. More recent
models take into account the possible early differentiation of
Europa's interior and expect a lower salt concentration of
maximum a few wt% (Melwani Daswani et al. 2021). To
model briny cryomagma, we vary the density of the liquid
and solid phases, respectively, equal to 920 and 1000 kg m⁻³
for pure water and 1130 and 1180 kg m⁻³ for briny
cryomagma (Kargel 1991). We also set the freezing
temperature to 268 K for the briny solution. Figure 11
shows the initial eruption volume for these two cryomagma
compositions.

For these two figures, panels (a) and (b) are obtained with the nondeformable wall approximation as in Lesage et al. (2020). The dark blue areas of these panels simply indicate the reservoirs that are too large to be stored in the 10 km conductive layer. Panels (c) and (d) are obtained using the iterative algorithm described in Section 2.5. The dark blue area of these panels indicates the reservoirs for which the algorithm diverges and thus that cannot be pressurized by freezing according to the algorithm. Finally, panels (e) and (f) of both Figures 10 and 11 are calculated using the iterative model, but reservoirs that cannot trigger eruptions are identified using the Maxwell time criterion; for each reservoir, if $\tau_M > \tau_c$, then the ice behaves as a viscoelastic material, and the reservoir is considered as unable to erupt.

Here we propose to demonstrate the similarity of the results
obtained with the deformation algorithm and the Maxwell time
criterion when identifying the reservoirs that can or cannot
trigger eruptions. In panels (e) and (f) of both Figures 10 and
11, we use white dashed lines to distinguish between reservoirs
that are stored in the elastic ice ($\tau_M > \tau_c$) and viscoelastic ice
($\tau_M < \tau_c$). These lines are copied in panels (c) and (d). We
propose that the comparison between $\tau_M$ and $\tau_c$ demonstrates
that a simple viscoelastic analysis is a convenient and efficient
method to evaluate the effects of viscoelastic deformation
without employing the iterative deformation model. We
demonstrate here that for a reservoir of a given radius and
depth, if $\tau_M < \tau_c$, the reservoir deformation is always sufficient
to accommodate the freezing overpressure and prevent an
eruption.

3.3. *Effect of Reservoir Geometry*

We finally model the freezing of lens-shaped reservoirs by
adapting our model to this new geometry as described in
Section 2.6. An ellipsoidal geometry is more realistic to model
putative sills or elongated reservoirs, which could better match
the extent of the geological features observed on Europa, such
as chaos, lenticulae, smooth plains, or double ridges (Michaut
& Manga 2014; Craft et al. 2016). For these reservoirs, we
cannot use the deformation model from Dragoni & Magnanensi
(1989), which only applies to spherical cavities. Nevertheless,
as demonstrated in the previous sections, reservoir deformation
can be neglected in the elastic part of the ice, and we thus
calculate the critical freezing time and erupted volume of lens-
shaped reservoirs using the nondeformable wall approximation.
This time, the reservoirs are modeled as ellipsoids of horizontal
semiaxes *a* and *b* and vertical semiaxis *c*, with an aspect ratio
$F$ ($a=b=F \cdot c$).

Figure 12 shows the results obtained with an aspect ratio $F = 500$. We consider reservoirs with horizontal extents up to several tens of kilometers. Compared to the spherical approximation, larger volumes can be stored in the elastic part of the ice with elongated reservoirs. We can see in Figure 12 that the erupted volume of cryolava goes up to $10^{13} \text{ m}^3$ for an elongated reservoir with $F = 500$, compared with $10^9 \text{ m}^3$ for a spherical one (see Figures 11(e) and (f)). This result demonstrates that elongated reservoirs or sills are necessary to explain the emplacement of large cryovolcanic features by freezing.

For example, Lesage et al. (2021) measured the volume of four relatively small (a few kilometers wide) smooth plains observed in Galileo images and found associated volumes of erupted cryolava up to $3 \cdot 10^9 \text{ m}^3$. This erupted volume corresponds to the limit of the plausible spherical reservoir volumes but can easily be explained with an elongated reservoir. The upper limit of $3 \cdot 10^9 \text{ m}^3$ of cryolava can be explained by the initial eruption from a 3 km deep or shallower spherical reservoir (Figures 11(e) and (f)) or a 6 km deep or shallower lens-shaped reservoir (Figure 12). The features studied in Lesage et al. (2021) are only a few kilometers wide, and larger potentially cryovolcanic features such as the lobate flows on the west side of Thrace Macula (47°S, 172°W; e.g., Fagents 2003) would thus imply shallower and larger lens-shaped reservoirs.

4. Discussion

If indeed freezing-induced cryolava eruptions have occurred
on Europa, our results provide insights into the potential source
depths of reservoir formation and eruption mechanisms.
Mechanisms leading to the formation of reservoirs in the
shallow part of the ice shell (≈5 km, varies with the reservoir
size and shape) could be compatible with freezing-induced
eruptions; deeper reservoirs would accommodate the over-
pressure by viscoelastic deformation of the surrounding ice.
Here we discuss those different hypotheses in view of our
results, which are summarized in Figure 13.

Several modeling studies have predicted that partial melting
occurs below Europa's conductive upper ice shell, in or on top
of the convective layer, where the ice is warmer and the
dissipation of tidal energy more efficient (Sotin et al. 2002;
Tobie et al. 2003; Han & Showman 2010; Vilella et al. 2020).
We demonstrate here that reservoirs formed by deep partial
Pure water:

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

(c) Freezing time with deformation

$t_{cv}$ (years)

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

(e) Freezing time with
deformation and Maxwell time filter

$t_{cv}$ (years)

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

81 wt% H₂O + 16 wt% MgSO₄ + 3wt% Na₂SO₄:

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

(d) Freezing time with deformation

$t_{cv}$ (years)

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

(f) Freezing time with
deformation and Maxwell time filter

$t_{cv}$ (years)

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

**Figure 10.** Critical freezing time necessary to trigger an eruption from a reservoir filled with ((a), (c), and (e)) pure liquid water and ((b), (d), and (f)) a mixture of liquid water and salts. The following situations are represented. (a) and (b) The reservoir wall does not deform during freezing (nondeformable wall approximation as in Lesage et al. 2020). The dark blue area represents reservoirs that are too large to be stored in a 10 km thick conductive ice crust. (c) and (d) The reservoir wall deformation is taken into account using the algorithm described in Section 2.5. The algorithm diverges for the reservoirs in the dark blue area; these reservoirs accommodate the overpressure by deformation and cannot trigger eruptions. (e) and (f) We take into account the wall deformation as in panels (c) and (d), but we add a filter using the Maxwell time criterion: the reservoirs for which $\tau_M < \tau_c$ cannot trigger eruptions (dark blue area). The white dashed line corresponds to the limit between the elastic and viscoelastic reservoirs (Maxwell time criterion) in panels (e) and (f). This line is reported in panels (c) and (d) to show the similarity with the reservoirs for which the algorithm converges (eruption possible) or diverges (overpressure accommodated by deformation).
Pure water:

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

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

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

81 wt% H₂O + 16 wt% MgSO₄ + 3wt% Na₂SO₄:

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

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

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

**Figure 11.** Volume of the first eruption of a reservoir filled with ((a), (c), and (e)) pure liquid water or ((b), (d), and (f)) a mixture of liquid water and salts. See caption of Figure 10 for the three different cases represented.
<figure><img src="image_25.png" /><figcaption>**Figure 12.** Critical freezing time and erupted volume for an ellipsoidal reservoir filled with pure water. The aspect ratio of the reservoir is 500 ($a = b = 500c$).</figcaption></figure>

melting, if they remain at depth, would be very unlikely to erupt because of the low ice viscosity involved and the large critical freezing time required, even when assuming that tidal heating ceases prior to freezing.<fnref n="4" />

Shallower reservoir formation mechanisms have been proposed in the literature, and eruptions from these reservoirs may be consistent with our results. Partial melting may be possible around tidally activated strike-slip faults, as explored by Hammond (2020); cyclic slip around fractures can produce frictional heat in the near surface, potentially enough to locally reach the water melting temperature. Simulations conducted by Kalousová et al. (2016) demonstrated that no more than 10% of the melt volume fraction could be stored in the shallow subsurface because of the efficient downward water percolation through the ice shell. These studies do not take into account antifreezing components like sulfates and chlorides that may be present within Europa's ice. Nevertheless, to allow the freezing of such a partially melted reservoir, heat input through tidal dissipation and advection (recharge from continuing activity) must decline or cease; our model thus only considers areas that are no longer geologically active.

Recent radar observations revealed the presence of a liquid sill beneath a double ridge in Greenland (Culberg et al. 2022). This double ridge may have been formed by freezing-induced fracturation and eruption of the reservoir as demonstrated by Culberg et al. (2022). In Lesage et al. (2021), we studied four smooth plains that are indeed located close to or cross-cutting double ridges. In the view of these new results, implying that elongated reservoirs are more likely to produce large deposits, potential relationships between ridges and long, shallow reservoirs may be of particular interest.

Craft et al. (2016) proposed another formation mechanism of shallow reservoirs that relates to the intrusion of sills in the ice shell. If these sills reach Europa's shallow lithosphere, they

could represent good candidates for freezing-induced eruptions. In previous studies, shallow reservoirs have been considered in relation to the formation of several geological structures, including pits, domes and microchaos (Michaut & Manga 2014), and double ridges (Johnston & Montési 2014; Culberg et al. 2022). One challenge still persists with these formation mechanisms, that is the refilling of the sills. A recent study demonstrated that fractures are unable to propagate over sufficient distance to link Europa's ocean and near subsurface (Walker et al. 2021), which makes the shallow sill formation by ocean injection questionable.

Finally, impact-induced melting has also recently been considered on Europa. Steinbrügge et al. (2020) demonstrated that impact craters can cause local melting and brine migration, which could have formed a cryomagma reservoir and triggered an eruption by freezing in Manannán crater, among others. This mechanism is also compatible with our results, as the impact-induced reservoir would be located in shallow elastic ice (around 4 km deep; Steinbrügge et al. 2020). This melting and freezing mechanism has also been invoked on Ceres to explain the lobate deposits observed at the bottom of Occator crater (Hesse & Castillo-Rogez 2019; Raymond et al. 2020).

In Lesage et al. (2021), we investigated the possibility of cyclic cryolava eruptions in the case where reservoir freezing continues after the first eruption. In that case, the erupted cryolava volume can theoretically reach about one-tenth of the reservoir volume. Nevertheless, because the new reservoir radius and the temperature of the ice around the reservoir vary, it would be valuable to recalculate the critical freezing time, $τ_c$, before each eruption and compare it with the Maxwell time, $τ_M$, in order to know if the next eruption can be triggered. This is particularly true for reservoirs close to the elastic-viscoelastic transition depth corresponding to their radius and shape. We expect that, as the far-field ice temperature surrounding a reservoir increases with depth, viscoelastic relaxation becomes important. The intereruption timescale also increases with temperature (and depth), as observed as an increase in the critical freezing times for a given reservoir volume. The total number of eruptions and total erupted

<footnote n="4">⁴ Nevertheless, as pointed out by Tobie et al. (2003), this first conclusion might be modulated in the case of a highly fractured crust that could produce enhanced tidal heating along preexisting faults and migration of the melt generated deeper.</footnote>
<figure><img src="image_26.png" /></figure>

**Figure 13.** Review of the mechanisms proposed in the literature to explain the emplacement of liquid reservoirs invoked to correspond to a variety of geological features at Europa's surface. Tidally activated strike-slip faults may produce local partial melting and draining, as modeled by Kalousová et al. (2016) and Hammond (2020). Meteorite impacts may generate melted brine that, after percolating, can form shallow lenses per Steinbrügge et al. (2020). Double ridges may form by freezing of shallow liquid water intrusions as modeled by Johnston & Montési (2014) and recently observed in Greenland (Culberg et al. 2022). A wide range of features, such as chaos (Schmidt et al. 2011), microchaos, lenticulae, pits and domes (Manga & Michaut 2017), and smooth plains (Lesage et al. 2021), may be morphologically explained by the presence of water sills or elongated reservoirs in the shallow subsurface. Craft et al. (2016) demonstrated that horizontal sills can propagate in the ice shell (this study is applied to the formation of double ridges but is not represented here for simplicity of the figure), but this mechanism is challenging. Connection between the sills and the subsurface ocean is unlikely, as demonstrated by Walker et al. (2021). Finally, deep in situ partial melting on top of convective ice plumes has been modeled in several studies, e.g., Vilella et al. (2020) and references therein. We indicate with a dashed line the elastic-viscoelastic transition determined by comparing the ice Maxwell time, $τ_M$, and the reservoir critical freezing time, $τ_c$, assuming that all reservoirs represented have an identical critical freezing time. The figure is not to scale.

volume thus decrease with depth because more of the ice shell freezes between eruptive episodes. Finally, as the volume of melt required to relieve the overpressure decreases with decreasing reservoir volume, we expect the volumes of individual eruptions to decrease with depth.

**5. Conclusion**

In the context of the freezing of a cryomagma reservoir of a given depth, volume, and geometry, we demonstrate that we can separate Europa's conductive upper icy shell into two distinct layers:

1. The uppermost part of the conductive layer, where the Maxwell time of the ice, $τ_M$, is greater than the reservoir critical freezing time, $τ_c$, and thus behaves as an elastic material at the critical freezing timescale. This part of the ice shell deforms slightly but not rapidly enough to accommodate the freezing overpressure. It is possible to trigger an eruption by freezing a reservoir stored in this part of the icy shell.

2. The lower part of the conductive layer, where $τ_M < τ_c$. This region behaves as a viscous material at the critical freezing timescale. It deforms easily and fully accommodates the

freezing overpressure. It is not possible to trigger an eruption from a reservoir stored in this part of the ice shell or anywhere at greater depths.

Our results indicate that the approximation of a nondeformable reservoir wall is realistic to model reservoirs for which $τ_M > τ_c$. These results may simplify further modeling of cryovolcanic eruptions by freezing.

Finally, we demonstrate that spherical reservoirs of more than a few kilometers in radius are unable to trigger eruptions, which limits their erupted volume to $10^{10} \text{ m}^3$ in that case, for the assumed ice shell structure. Elongated ellipsoidal reservoirs may explain larger erupted volumes, up to $10^{14} \text{ m}^3$ for reservoirs with a realistic aspect ratio of 500 between the horizontal and vertical semiaxes, and potentially more for greater aspect ratios. Depending on the reservoir geometry, the maximum depth of active reservoirs is limited to $\lesssim 75\%$ of the conductive ice shell thickness.

The upcoming NASA Europa Clipper (Howell & Pappalardo 2020) and ESA JUICE (Grasset et al. 2013) missions will allow for the correlation of potential surface features and plumes, seen in both visible and spectroscopic imagery, and subsurface structures, observed through radar, to better
constrain the evolution of potential reservoirs. This study represents a step in constraining the source depth and geometry of subsurface reservoirs using our understanding of the ice shell thermal and mechanical state. If erupting reservoirs are discovered and found to be relegated to the shallow depths predicted in this study, they may represent the most accessible candidate target for the in situ exploration of liquid water environments beyond Earth.

We acknowledge support from the Institut National des Sciences de l'Univers (INSU), France, the Centre National de la Recherche Scientifique (CNRS), France, and the Centre National d'Etudes Spatiales (CNES), France, through the Programme National de Planétologie, France. 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 partially supported by Scientific Exploration Subsurface Access Mechanism for Europa (SESAME) grant No. 80NSSC19K0614. The code developed and used in this study is available at https://github.com/ElodieLesage/CryoFRIES.

## ORCID iDs

Elodie Lesage https://orcid.org/0000-0001-5425-9665

Hélène Massol https://orcid.org/0000-0002-0599-3844

Samuel M. Howell https://orcid.org/0000-0002-5126-3228

Frédéric Schmidt https://orcid.org/0000-0002-2857-6621

## References

Bierhaus, E. B., Zahnle, K., Chapman, C. R., et al. 2009, in Europa, ed. R. T. Pappalardo et al. (Tucson, AZ: Univ. Arizona Press), 161

Biot, M. A. 1954, JAP, 25, 1385

Craft, K. L., Patterson, G. W., Lowell, R. P., & Germanovich, L. 2016, Icar, 274, 297

Culberg, R., Schroeder, D. M., & Steinbrügge, G. 2022, NatCo, 13, 2007

Doggett, T., Greeley, R., Figueredo, P., & Tanaka, K. 2009, in Europa, ed. R. T. Pappalardo et al. (Tucson, AZ: Univ. Arizona Press), 137

Dombard, A. J., Patterson, G. W., Lederer, A. P., & Prockter, L. M. 2013, Icar, 223, 74

Dragoni, M., & Magnanensi, C. 1989, PEPI, 56, 316

Fagents, S. A. 2003, JGR, 108, 5139

Folch, A., Fernández, J., Rundle, J. B., & Martí, J. 2000, GeoJI, 140, 37

Goldsby, D., & Kohlstedt, D. L. 2001, JGR, 106, 11017

Grasset, O., Dougherty, M., Coustenis, A., et al. 2013, P&SS, 78, 1

Greeley, R., Sullivan, R., Klemaszewski, J., et al. 1998, Icar, 135, 4

Gregg, P., de Silva, S., & Grosfils, E. 2013, E&PSL, 384, 100

Hammond, N. P. 2020, JGRE, 125, e06170

Han, L., & Showman, A. P. 2010, Icar, 207, 834

Hesse, M. A., & Castillo-Rogez, J. C. 2019, GeoRL, 46, 1213

Hillier, J., & Squyres, S. W. 1991, JGR, 96, 15665

Hodge, D. S. 1974, Natur, 251, 297

Howell, S. M. 2021, PSJ, 2, 129

Howell, S. M., & Pappalardo, R. T. 2018, GeoRL, 45, 4701

Howell, S. M., & Pappalardo, R. T. 2020, NatCo, 11, 1311

Jellinek, A. M., & DePaolo, D. J. 2003, BVol, 65, 363

Johnston, S. A., & Montési, L. G. 2014, Icar, 237, 190

Kalousová, K., Souček, O., Tobie, G., Choblet, G., & Čadek, O. 2014, JGRE, 119, 532

Kalousová, K., Souček, O., Tobie, G., Choblet, G., & Čadek, O. 2016, JGRE, 121, 2444

Kargel, J. S. 1991, Icar, 94, 368

Kattenhorn, S. A., & Prockter, L. M. 2014, NatGe, 7, 762

Khurana, K. K., Kivelson, M. G., Stevenson, D. J., et al. 1998, Natur, 395, 777

Kuskov, O., & Kronrod, V. 2005, Icar, 177, 550

Lesage, E., Massol, H., & Schmidt, F. 2020, Icar, 335, 113369

Lesage, E., Schmidt, F., Andrieu, F., & Massol, H. 2021, Icar, 361, 114373

Litwin, K. L., Zygielbaum, B. R., Polito, P. J., Sklar, L. S., & Collins, G. C. 2012, JGRE, 117, E08013

Manga, M., & Michaut, C. 2017, Icar, 286, 261

Manga, M., & Wang, C.-Y. 2007, GeoRL, 34, L07202

McTigue, D. F. 1987, JGR, 92, 12931

Melwani Daswani, M., Vance, S. D., Mayne, M. J., & Glein, C. R. 2021, GeoRL, 48, e94143

Michaut, C., & Manga, M. 2014, JGRE, 119, 550

Neveu, M., Desch, S., Shock, E., & Glein, C. 2015, Icar, 246, 48

Quick, L. C., Glaze, L. S., & Baloga, S. M. 2017, Icar, 284, 477

Raymond, C. A., Ermakov, A. I., Castillo-Rogez, J. C., et al. 2020, NatAs, 4, 741

Rudolph, M. L., Manga, M., Walker, M., & Rhoden, A. R. 2022, GeoRL, 49, e94421

Schmidt, B. E., Blankenship, D. D., Patterson, G. W., & Schenk, P. M. 2011, Natur, 479, 502

Segall, P. 2016, JGRA, 121, 8501

Sigurdsson, H., Houghton, B., Rymer, H., Stix, J., & McNutt, S. 1999, Encyclopedia of Volcanoes (New York: Academic)

Sotin, C., Head, J. W., & Tobie, G. 2002, GeoRL, 29, 1233

Steinbrügge, G., Voigt, J. R., Wolfenbarger, N. S., et al. 2020, GeoRL, 47, e90797

Tobie, G., Choblet, G., & Sotin, C. 2003, JGRE, 108, 5124

Vilella, K., Choblet, G., Tsao, W.-E., & Deschamps, F. 2020, JGRE, 125, e06248

Walker, C. C., Bassis, J. N., & Schmidt, B. E. 2021, PSJ, 2, 135

Zahnle, K., Schenk, P., Levison, H., & Dones, L. 2003, Icar, 163, 263