Simulation of Freezing Cryomagma Reservoirs in Viscoelastic Ice Shells

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₁ to be the reservoir radius and R₂ the viscoelastic shell radius (see Figure 1).

Zoom In Zoom Out Reset image size

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.

Table 1. Parameters Used in This Study

Parameter	Description	Value
cp	Specific heat capacity of ice	2.3 J kg−1 K−1
E	Young modulus	3 · 109 Pa
G	Shear modulus	3 · 109 Pa
H	Reservoir depth	1–10 km
k	Thermal conductivity of ice	2.3 W m−1 K−1
K	Bulk modulus	$\tfrac{E}{3\cdot (1-2\nu )}$
R1	Reservoir radius	10 m–5 km for spherical reservoirs, up to 100 km for lens-shaped reservoirs
R2	Viscoelastic shell radius	See Figures 1 and 4 and Equation (9)
Tcold	Temperature in the ice crust far from any reservoir	Linear gradient from 110 (surface) to 250 K (bottom of conductive layer, 10 km deep)
Tm	Cryomagma melting temperature	273 K for pure water, 268 K for briny cryomagma (Kargel 1991)
T(R2)	Temperature at which ice becomes viscoelastic around a reservoir	See Section 2.3 and Equation (9)
Ve	Erupted volume with nondeformable wall approximation	See Equation (6)
Vev	Erupted volume taking into account the elastic and viscoelastic deformation	Estimated using the iterative scheme described in Section 2.5
ΔPc	Critical overpressure necessary to fracture the reservoir with the nondeformable wall approximation	2(σc + ρs gH), see Equation (4)
η	Viscosity	See Equation (7) (applies to viscoelastic region only, i.e., R1 < r < R2)
κ	Thermal diffusivity of ice	k/(ρs cp )
ν	Poisson ratio	0.325
ρl	Density of liquid cryomagma	1000 kg m−3 for pure water, 1180 kg m−3 for briny cryomagma
ρs	Density of ice	920 kg m−3 for pure water ice, 1130 kg m−3 for briny ice
σc	Tensile strength of ice	Linear gradient from 1.7 · 106 (surface) to 1 · 106 (10 km deep; Litwin et al. 2012)
τ	Viscoelastic deformation time	See Dragoni & Magnanensi (1989), Equation (9)
τM	Maxwell time	η/E
τc	Critical freezing time necessary to fracture the reservoir with the nondeformable wall approximation	See Equation (5)
τcv	Critical freezing time necessary to fracture the reservoir taking into account the ice elastic and viscoelastic deformation	Estimated using the iterative scheme described in Section 2.5
χ	Liquid water compressibility	5 · 10−10 Pa−1

Download table as:  ASCIITypeset image

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 ΔP in the reservoir is

$$\begin{eqnarray}&&{\rm{\Delta }}P=-\displaystyle \frac{1}{\chi }\mathrm{ln}\left(\displaystyle \frac{1-n\left({\rho }_{l}/{\rho }_{s}\right)}{1-n}\right),\end{eqnarray} \tag{ 1 }$$

where χ is the liquid water compressibility, n is the cryomagma frozen fraction, and ρ_l and ρ_s are, respectively, the densities of the cryomagma and associated ice. The frozen fraction n is calculated as a function of time by solving the Stefan problem (see Lesage et al. 2020, Appendix A, for further details),

$$\begin{eqnarray}&&n=1-{\left(1-\displaystyle \frac{S}{{R}_{1}}\right)}^{3}.\end{eqnarray} \tag{ 2 }$$

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

$$\begin{eqnarray}&&S=2\lambda \sqrt{\kappa t},\end{eqnarray} \tag{ 3 }$$

with λ being a constant given by solving the Stefan problem and κ 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.

Zoom In Zoom Out Reset image size

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))

$$\begin{eqnarray}&&{\rm{\Delta }}{P}_{c}=2({\sigma }_{c}+{\rho }_{s}{gH}),\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}&&{\tau }_{c}=\left(\frac{{S}_{c}}{2\lambda \sqrt{\kappa }}\right),\end{eqnarray} \tag{ 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,

$$\begin{eqnarray}&&{V}_{e}={V}_{i}{n}_{c}\displaystyle \frac{{\rho }_{l}}{{\rho }_{s}},\end{eqnarray} \tag{ 6 }$$

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

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 η of the viscoelastic ice as a function of its temperature T, we use the Newtonian law of Hillier & Squyres (1991):

$$\begin{eqnarray}&&\eta ={10}^{14}\exp \left\{25.2\left(273/T-1\right)\right\}.\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&\bar{\eta }={\left(\displaystyle \frac{1}{n}\sum _{i=1}^{n}\displaystyle \frac{1}{\eta ({T}_{i})}\right)}^{-1},\end{eqnarray} \tag{ 8 }$$

where T₁ 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 ≪ 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, τ_M = η/E, with E being the ice Young modulus and η the temperature-dependent viscosity; and (ii) the critical freezing time τ_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 τ_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 τ_M and τ_c for a 500 m radius reservoir. One can see that in the first 6 km of ice extending from the surface, τ_M ≫ τ_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, τ_M ≪ τ_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.

Zoom In Zoom Out Reset image size

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₂ − R₁. This can be accomplished by finding the location r = R₂, where ice is cold enough to behave as a purely elastic material, i.e., where τ_M becomes greater than τ_c . To determine R₂, we first find the temperature T(R₂) at which τ_M = τ_c using Equation (7) to represent the ice viscosity. The radius R₂ can then be determined by finding the location where the temperature around the reservoir equals T(R₂). An expression for R₂ 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₂),

$$\begin{eqnarray}&&\begin{array}{c}{R}_{2}={R}_{1}-2\sqrt{{\kappa }{{\tau }}_{c}}{{\rm{e}}{\rm{r}}{\rm{f}}}^{-1}\\ \,\times \,\left(\displaystyle \frac{T({R}_{2})-{T}_{{\rm{c}}{\rm{o}}{\rm{l}}{\rm{d}}}(H)}{{T}_{m}-{T}_{{\rm{c}}{\rm{o}}{\rm{l}}{\rm{d}}}(H)}\left(1+{\rm{e}}{\rm{r}}{\rm{f}}({\lambda })\right)-1\right),\end{array}\end{eqnarray} \tag{ 9 }$$

where κ is the ice thermal diffusivity, T_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 · R₁. 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.

Zoom In Zoom Out Reset image size

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 (σ_rr (r, t)) and two tangential (σ_{θ θ} (r, t) and σ_{ϕ ϕ} (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), σ_rr (r, t), σ_{θ θ} (r, t), and σ_{ϕ ϕ} (r, t) are the following:

$$\begin{eqnarray}&&u(r,t)=\displaystyle \frac{{A}_{1}(r)}{{D}_{0}}{f}_{1}(t)+\displaystyle \frac{{A}_{2}(r)}{{D}_{0}}{f}_{2}(t),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\sigma }_{{rr}}(r,t)=\displaystyle \frac{{B}_{1}(r)}{{D}_{0}}{f}_{1}(t)+\displaystyle \frac{{B}_{2}(r)}{{D}_{0}}{f}_{2}(t),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\sigma }_{\theta \theta }(r,t)={\sigma }_{\phi \phi }(r,t)=\displaystyle \frac{{C}_{1}(r)}{{D}_{0}}{f}_{1}(t)+\displaystyle \frac{{C}_{2}(r)}{{D}_{0}}{f}_{2}(t).\end{eqnarray} \tag{ 12 }$$

The coefficients A₁, A₂, B₁, B₂, C₁, C₂, and D₀ and the two functions f₁(t) and f₂(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 μ in Dragoni & Magnanensi (1989). We use Equation (7) to calculate these coefficients.

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

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}P(t) & = & {\rm{\Delta }}{P}_{c}\left\{\displaystyle \frac{t}{{t}_{1}}\left[1-H(t-{t}_{1})\right]\right.+\,H(t-{t}_{1})-H(t-{t}_{2})\\ & & +\left.\displaystyle \frac{{t}_{3}-{t}_{1}}{{t}_{3}-{t}_{2}}\left[H(t-{t}_{2})-H(t-{t}_{3})\right]\right\},\end{array}\end{eqnarray} \tag{ 13 }$$

where t = 0 is the beginning of the pressure increase, t = t₁ is the beginning of the pressure plateau, t = t₂ is the beginning of the pressure decrease, and t = t₃ 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 = τ_c when the eruption begins. We set t₁ = t₂ = t₃ = τ_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).

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 derive an analytical form of the unique realistic critical freezing time to reach the critical overpressure, ΔP_c , and trigger an eruption. To solve this problem and find the critical freezing time, hereafter called τ_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, τ_c , of a reservoir and the associated overpressure, ΔP_c , by making the nondeformable wall assumption using the model from Lesage et al. (2020).
• 2.  
  We calculate the corresponding deformation u(R₁, τ_c ) from Equations (10)–(12) at time t = τ_c and the induced pressure drop, δ P, using the definition of the water compressibility, $\delta P=-\tfrac{1}{\chi }\tfrac{{V}_{f}}{{V}_{i}}$, with χ 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 ΔP_c + δ P and simulate new reservoir deformation after this longer freezing time, τ_c + δ τ. This deformation induces a new pressure drop of δ P₂.
• 4.  
  We verify whether δ P − δ P₂ < , where is an imposed convergence criterion; we use = ΔP_c · 10⁻⁵. If the model converges, the algorithm has found the realistic critical freezing time τ_cv = τ_c + δ τ to reach the critical overpressure, ΔP_c . If not, the algorithm iterates beginning with step 3, replacing δ P with δ P₂ and δ τ with δ τ₂, and then verifies the convergence, and so on.

We iterate until the program converges to a realistic critical freezing time, τ_cv, that takes into account the feedback between pressurization and deformation. As an example, Figure 6 shows the pressure drop, δ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In some cases, the deformation is great enough to completely accommodate the freezing overpressure at each iteration, and the algorithm cannot find a realistic τ_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 τ_c > τ_M. Hence, perhaps intuitively, the Maxwell time of the ice can be used to discriminate reservoirs that can or cannot trigger eruptions.

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¹² 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 · 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 τ_c and Maxwell time τ_M to determine whether an eruption is possible or not, the efficacy of which is demonstrated in Section 2.5.

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₂ equal to half of the thickness of Europa's viscoelastic ice shell (R₂ = 2 km here), and we only model the region spanning R₁ < r < R₂.

The deformation and stress fields around the reservoirs at t = τ_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 = τ_c , which is the critical freezing time under the nondeformable wall approximation calculated at the first iteration of our algorithm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As shown in Figure 7(a), reservoir A (elastic case) is lightly deformed due to freezing; after the critical freezing time τ_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 = τ_c . As σ_{θ θ} reaches the critical value Δ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 = τ_c but only after a longer critical freezing time τ_cv. In the particular example of reservoir A, we obtained τ_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 τ_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, σ_{θ θ} , never reaches the critical value, Δ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 iterations following model divergence to save time and conserve sanity.

The results presented in this section are calculated with the iterative process introduced in Section 2.5 to calculate realistic critical freezing times τ_cv that take into account the feedback between pressure and deformation. Figure 9 shows the critical freezing times τ_c and how they compare with the realistic critical freezing times τ_cv calculated with the nondeformable wall approximation. We show τ_c and τ_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 τ_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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 τ_M > τ_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 (τ_M > τ_c ) and viscoelastic ice (τ_M < τ_c ). These lines are copied in panels (c) and (d). We propose that the comparison between τ_M and τ_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 τ_M < τ_c , the reservoir deformation is always sufficient to accommodate the freezing overpressure and prevent an eruption.

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 · 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¹³ m³ for an elongated reservoir with F = 500, compared with 10⁹ m³ 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.

Zoom In Zoom Out Reset image size

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 · 10⁹ m³. 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 · 10⁹ m³ 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.