# Supplement: ``A mushy source for the geysers of Enceladus''

Colin R. Meyer<EMAIL_ADDRESS>Thayer School of Engineering,
Dartmouth College, Hanover, NH 03755  Jacob J. Buffo  Thayer School of
Engineering, Dartmouth College, Hanover, NH 03755  Francis Nimmo  Department
of Earth & Planetary Sciences, University of California, Santa Cruz, CA 95064
USA  Andrew J. Wells  Atmospheric, Oceanic & Planetary Physics, Department of
Physics, University of Oxford, OX1 3PU, UK  Samuel Boury  Courant Institute of
Mathematical Sciences, New York University, New York, NY 10012  Tara C.
Tomlinson  Thayer School of Engineering, Dartmouth College, Hanover, NH 03755
Jamie R. G. Parkinson  Atmospheric, Oceanic & Planetary Physics, Department of
Physics, University of Oxford, OX1 3PU, UK  Geoffrey M. Vasil  School of
Mathematics, University of Edinburgh, EH9 3FD, UK

##  1 Fracture depths

It is unclear how easy it is to produce a fracture spanning the entire ice
shell of Enceladus. Due to the limited tidal stresses ($\sim 100$ kPa), when
compared to the overburden pressure ($\sim 600$ kPa at the base of a 6 km ice
shell), it is hard for near-surface fractures to extend downwards more than
1–2 km (Smith-Konter and Pappalardo, 2008; Olgin et al., 2011; Walker et al.,
2021). This is less than the estimated minimum shell thickness of 6 km at the
south pole (Hemingway et al., 2018). Fractures extending upwards from the base
of the shell are less affected by overburden pressure because they would be
water-filled (Crawford and Stevenson, 1988; Shibley and Laughlin, 2021).
However, ice at the base of the shell is more ductile and thus harder to
fracture than ice at the surface (Schulson and Duval, 2009).

The current paradigm for the source of geyser material is that fractures
extend down to the subsurface ocean (Hurford et al., 2007). To get an idea for
the depth of the tiger stripe fractures, we balance the extensional stress
required to open a fracture with the overburden ice pressure acting to close
it (i.e. the zero-stress Nye (1955) criterion for closely spaced glacier
crevasses). Estimates for the extensional stress are on the order of $\tau\sim
1$ bar (Nimmo et al., 2007, 2011). Gravity on Enceladus is $g=0.113$ m s-2
(Hansen et al., 2006, 2008) and the salty ice density is about
$\rho_{\textrm{ice}}=940$ kg m-3 (Timco and Frederking, 1996). The fracture
close-off depth $h_{\textrm{top}}$ is given as

$h_{\textrm{top}}\sim\frac{\tau}{\rho_{\textrm{ice}}g}\sim 1~{}\mathrm{km},$
(1)

which is an estimate for the likely depth of surface fractures and agrees with
more sophisticated calculations (Smith, 1976; Broberg, 1999; Smith-Konter and
Pappalardo, 2008; Olgin et al., 2011). The fact that the surface fractures are
only likely to propagate to $\sim 1~{}\mathrm{km}$ depth is important because
$\sim 1~{}\mathrm{km}$ is shallower than the minimum predicted shell thickness
of $\sim 6-30~{}\mathrm{km}$ (Hemingway et al., 2018), which casts doubt on
the idea that the fractures penetrate from the surface down to the ocean.
Moreover, it would be a considerable feat for a fracture to extend up from the
ocean (i.e. a basal crevasse) to within $\sim 1~{}\mathrm{km}$ of the surface,
based on the extensional stress just mentioned. A water-filled basal crevasse
does lead to a larger close-off height $h_{\textrm{bottom}}$, given by

$h_{\textrm{bottom}}\sim\frac{\tau}{\left(\rho_{\textrm{water}}-\rho_{\textrm{ice}}\right)g}\sim
10\mbox{ km},$ (2)

where we use a value of seawater density, $\rho_{\textrm{water}}=1030$ kg m-3.
This is sufficient to fracture the shell if (i) the fracture emanates up from
the bottom and (ii) if the shell is on the lower end of thickness range (i.e.
$\sim 6-30~{}\mathrm{km}$ Hemingway et al., 2018). However, the lower shell is
likely ductile and less prone to fracturing (Nimmo, 2004; Schulson and Duval,
2009). On the other hand, there are additional sources of stress available to
open cracks which likely exceed the tidal stresses. These include the effects
of shell thickening (Rudolph et al., 2022) and flexural stresses from plume
loading (Hemingway et al., 2020). Thus, opening of through-going fractures
appears possible but not certain. At the same time, the dike mechanism
discussed in the main text could allow for the initiation of downward-
propagating through-going fractures that could be maintained by turbulent
water flow and shear heating. We discuss the dike mechanism further in section
5.

##  2 Analytical temperature solution

Here we demonstrate that the leading order thermal structure can be understood
from steady state conduction, with reactive flow causing subsequent
modifications to the porosity structure, composition and shape of the
inclusion. We consider a fracture of depth $h$ in a shell of thickness $H$
with a shear heat flux of $\mu\rho_{i}ghu$. So long as the latent heat sink of
partial melting and heat transport by interstitial flow are relatively small,
then the steady state temperature $T$ around the fracture is given by the
solution to Laplace's equation, i.e.

$\nabla^{2}T=0,$ (3)

where we treat the thermal conductivity as temperature independent to
facilitate analytical progress. Equation (3) is subject to the boundary
conditions

$\displaystyle-k\frac{\partial T}{\partial x}$ $\displaystyle=$ $\displaystyle
0~{}~{}~{}\mbox{on}~{}~{}~{}x=0~{}(z<h)~{}~{}~{}\mbox{and}~{}~{}~{}x=L,$ (4)
$\displaystyle-k\frac{\partial T}{\partial x}$ $\displaystyle=$
$\displaystyle\mu\rho_{i}ghu~{}~{}~{}\mbox{on}~{}~{}~{}x=0~{}(z\geq h),$ (5)
$\displaystyle-k\frac{\partial T}{\partial z}$ $\displaystyle=$ $\displaystyle
r_{1}+r_{2}(T-T_{s})~{}~{}~{}\mbox{on}~{}~{}~{}z=H,$ (6) $\displaystyle T$
$\displaystyle=$ $\displaystyle T_{e}~{}~{}~{}\mbox{on}~{}~{}~{}z=0.$ (7)

The surface condition is a linearized radiative condition, with respect to the
background temperature $T_{s}$ of space near Enceladus (see table 1), rather
than the full black-body radiation condition as used by Abramov and Spencer
(2009). This choice for the surface boundary condition allows for an
analytical solution and is the condition used in SOFTBALL. To derive the
linearized condition, we start from the surface energy balance where
conductive heat fluxes from the ice balance net black-body radiation as

$-k\frac{\partial T}{\partial z}=\epsilon\sigma\left(T^{4}-T_{s}^{4}\right),$
(8)

where $\epsilon$ is the emissivity and $\sigma$ is the Stefan-Boltzmann
constant (values in table 1). Expanding the temperature as a small deviation
from the space temperature, i.e. $T=T_{s}+\tilde{T}$, where $\tilde{T}\ll
T_{s}$, then we can use the binomial expansion and discard terms of second
order and higher, to find that

$-k\frac{\partial T}{\partial z}=4\epsilon\sigma T_{s}^{3}(T-T_{s}),$ (9)

which is in the form of the boundary condition above, with $r_{1}=0$ and
$r_{2}=4\epsilon\sigma T_{s}^{3}$.

We nondimensionalized the problem, by taking

$T=T_{e}+\Delta T\theta,~{}~{}~{}x=H\tilde{x},~{}~{}~{}z=H\tilde{z},$

which yields

$\tilde{\nabla}^{2}\theta=0,$ (10)

subject to the boundary conditions

$\displaystyle-\frac{\partial\theta}{\partial\tilde{x}}$ $\displaystyle=$
$\displaystyle
0~{}~{}~{}\mbox{on}~{}~{}~{}\tilde{x}=0~{}(\tilde{z}<\delta)~{}~{}~{}\mbox{and}~{}~{}~{}\tilde{x}=\lambda,$
(11) $\displaystyle-\frac{\partial\theta}{\partial\tilde{x}}$ $\displaystyle=$
$\displaystyle
F~{}~{}~{}\mbox{on}~{}~{}~{}\tilde{x}=0~{}(\tilde{z}\geq\delta),$ (12)
$\displaystyle-\frac{\partial\theta}{\partial\tilde{z}}$ $\displaystyle=$
$\displaystyle
G+b\left(\theta-\theta_{s}\right)~{}~{}~{}\mbox{on}~{}~{}~{}\tilde{z}=1,$ (13)
$\displaystyle\theta$ $\displaystyle=$ $\displaystyle
0~{}~{}~{}\mbox{on}~{}~{}~{}\tilde{z}=0.$ (14)

Here we define the parameters as

$\delta=\frac{h}{H},~{}~{}~{}\lambda=\frac{L}{H},~{}~{}~{}F=\frac{\mu\rho_{i}ghHu}{k_{o}\Delta
T},~{}~{}~{}G=\frac{Hr_{1}}{k_{o}\Delta
T},~{}~{}~{}b=\frac{Hr_{2}}{k_{o}},~{}~{}~{}\theta_{s}=\frac{T_{s}-T_{e}}{\Delta
T}.$

Due to convergence issues in the SOFTBALL simulations, we could not use large
values for $b$, and so the value quoted in table 1 is artificially lower than
what we would calculate from the parameters. The main effect of this
difference is on the surface temperature profile and not the size, shape, or
liquid volume of the mushy zone. Asymptotically, our analytical solutions
shows that large values of $b$ pin the surface temperature to the space
temperature, except very near the crack, where there is a narrow boundary
layer region of rapid temperature change. The boundary layer becomes thinner
as $b$ increases. Additionally, since $G=0$ here, we disregard it in
subsequent calculations.

To simplify the problem in the vertical direction, we write the total solution
$\theta$ as a sum of a particular solution $\theta_{I}$ to the inhomogeneous
problem and the homogenous solution $\theta_{H}$, i.e.

$\theta=\theta_{I}(\tilde{z})+\theta_{H}(\tilde{x},\tilde{z}).$

Solving first for the inhomogeneous part, we have that

$\frac{\partial^{2}\theta_{I}}{\partial\tilde{x}^{2}}+\frac{\partial^{2}\theta_{I}}{\partial\tilde{z}^{2}}=0\longrightarrow\frac{\partial^{2}\theta_{I}}{\partial\tilde{z}^{2}}=0,$

so

$\theta_{I}=A\tilde{z}+B.$

Using the boundary conditions, we find that

$B=0\ \ \ \mbox{ and }\ \ \ A=\frac{b\theta_{s}}{1+b}.$

For large values of $b$, the slope $A$ is approximately $\theta_{s}$, which is
consistent with the asymptotic description above, that the surface temperature
mostly stays close to the space temperature.

Now we can examine the homogeneous problem using separation of variables, i.e.

$\theta_{H}(\tilde{x},\tilde{z})=f(\tilde{x})g(\tilde{z}),~{}~{}~{}\mbox{which
implies}~{}~{}~{}f^{\prime\prime}g+fg^{\prime\prime}=0.$

Dividing through and choosing the separation constant $\alpha^{2}$, we find
the two ODEs

$f^{\prime\prime}=\alpha^{2}f\ \ \ \mbox{ and }\ \ \
g^{\prime\prime}=-\alpha^{2}g.$

Solving the vertical equation first, we have that

$g=C\sin(\alpha\tilde{z})+D\cos(\alpha\tilde{z}).$

With $g=0$ at $\tilde{z}=0$, we know that $D=0$. The homogeneous surface
boundary condition is

$-g^{\prime}=bg\ \ \ \mbox{ {at} }\ \ \ \tilde{z}=1,$

which results in the transcendental eigenvalue condition

$-\alpha\cos(\alpha)=b\sin(\alpha).$ (15)

This yields a full set of distinct eigenvalues $\alpha_{n}$ for $n=1,2,3,...$,
which can be computed numerically. In the horizontal direction, we have that

$f=De^{\alpha\tilde{x}}+Ee^{-\alpha\tilde{x}}.$

Using the no heat flux boundary condition at $\tilde{x}=\lambda$, we have that

$-f^{\prime}=0\ \ \ \mbox{ {at} }\ \ \ \tilde{x}=\lambda\ \ \ \mbox{ {implies}
}\ \ \ E=De^{2\alpha\lambda}.$

Thus, the full solution with unknown coefficients is given by

$\theta=\frac{b\theta_{s}}{1+b}\tilde{z}+\sum_{n=1}^{\infty}{c_{n}\left[e^{\alpha_{n}\tilde{x}}+e^{2\alpha_{n}\lambda}e^{-\alpha_{n}\tilde{x}}\right]\sin(\alpha_{n}\tilde{z})}.$
(16)

We find the values for $c_{n}$ by using orthogonality and integrating over the
heat flux condition at $\tilde{x}=0$, which gives

$\int_{\delta}^{1}{F\sin(\alpha_{m}\tilde{z})~{}d\tilde{z}}=\sum_{n=1}^{\infty}{c_{n}\alpha_{n}\left[1-e^{2\alpha_{n}\lambda}\right]\int_{0}^{1}{\sin(\alpha_{n}\tilde{z})\sin(\alpha_{m}\tilde{z})~{}d\tilde{z}}}.$
(17)

The integral on right is zero for $\alpha_{n}\neq\alpha_{m}$, using the
transcendental eigenvalue condition from equation (15) implying orthogonality
of the eigenfunctions $\sin(\alpha_{n}\hat{z})$. When $\alpha_{n}=\alpha_{m}$,
we have that

$\int_{0}^{1}{\sin(\alpha_{m}\tilde{z})\sin(\alpha_{m}\tilde{z})~{}d\tilde{z}}=\frac{1}{2}-\frac{\sin(2\alpha_{m})}{4\alpha_{m}}.$

We evaluate the integral on the left as

$\int_{\delta}^{1}{F\sin(\alpha_{m}\tilde{z})~{}d\tilde{z}}=-\frac{F}{\alpha_{m}}\left[\cos(\alpha_{m})-\cos(\alpha_{m}\delta)\right].$

Putting it all together gives,

$c_{n}=-\frac{4F}{\left[2\alpha_{n}^{2}-\alpha_{n}\sin(2\alpha_{n})\right]}\frac{\left[\cos(\alpha_{n})-\cos(\alpha_{n}\delta)\right]}{1-e^{2\alpha_{n}\lambda}}.$
(18)

We show the analytical solution in figure 1 and find that a mushy zone will
form along the crack face in regions with $T>T_{e}$. These results are similar
to what we find in the main text as well as Hammond (2020), Kalousová et al.
(2016), and Stevenson (1996). This analytical solution also allows us to
predict the size of the mushy zone, by looking for the contour where
$T=T_{e}$.

Figure 1: Analytical solution for the temperature around a heated crack with
$F=260$.

Figure 2: Liquid volume in the mushy zone for $F=500$: (a) evolution to a
steady state and average effective heat flux as well as (b) rate of liquid
volume production, which can support the entire flux observed by the Cassini
spacecraft from a single tiger stripe fracture.

##  3 Mushy zone physics

To analyze the effect of shear heating on a fracture in an icy shell, we use
SOFTBALL (SOlidification, Flow and Thermodynamics in Binary ALLoys), a mushy
layer code with adaptive mesh refinement. The github repository for this open-
source code is: https://github.com/jrgparkinson/mushy-layer (accessed 8 April
2023). This state-of-the-art code solves the equations of mushy layer theory,
as described below, in two-dimensional geometries (Wells et al., 2019;
Parkinson et al., 2020). The code uses a finite volume discretization and
adaptive time steps. In our simulations we use a spatial scale of $N_{x}=512$
grid points in the horizontal direction and $N_{z}=256$ in the vertical
direction. This gives a resolution of $\sim 40$ meters in both the horizontal
and vertical directions. The time steps adjust and are on the order of $\sim
30$ years. The simulations take on the order of a day to run on the linux
system (2-16 cores) at Dartmouth College.

Conservation of mass in a mushy layer dictates that

$\frac{\partial\left(\phi\rho\right)}{\partial
t}=M_{f}~{}~{}~{}\mbox{and}~{}~{}~{}\frac{\partial\left[\left(1-\phi\right)\rho\right]}{\partial
t}+\boldsymbol{\nabla}\cdot\boldsymbol{q}=-M_{f},$ (19)

where the solid and liquid densities $\rho$ are approximated as equal (using
an extended form of the Boussinesq approximation; Parkinson et al., 2020), the
solid fraction is $\phi$, the fluid flux relative to the solid is
$\boldsymbol{q}$, and the rate of mass phase change per volume is $M_{f}$. The
sum of equations (19) yields $\boldsymbol{\nabla}\cdot\boldsymbol{q}=0$, i.e.
a solenoidal fluid flow field. Here we assume that the porous solid matrix is
stationary. Fluid flux is a function of pressure gradients and gravity, which
can be described by Darcy's law, as

$\boldsymbol{q}=-\frac{\Pi}{\eta}\left(\nabla
p+\rho_{b}g\boldsymbol{z}\right),$ (20)

where $\eta$ is the liquid viscosity, $\rho_{b}$ is the density of liquid
brine (which is a function of the brine concentration and temperature),
$\boldsymbol{z}$ is the vertical unit vector, and $\Pi$ is the solid-fraction-
dependent permeability (Tait and Jaupart, 1992; Katz and Worster, 2008), given
by

$\Pi=\left[\frac{12}{d^{2}}+\frac{1}{K_{0}\phi^{3}}\right]^{-1}.$ (21)

Following Katz and Worster (2008), SOFTBALL uses a permeability regularization
for pure liquid regions (i.e. as $\phi\rightarrow 1$) such that there is a
background permeability and Darcy's law is solved everywhere, rather solving
Navier-Stokes in the pure fluid region (see Parkinson et al., 2020; Buffo et
al., 2021). The ratio of the background permeability $d^{2}/12$ to the
permeability prefactor $K_{0}$ is the reluctance $\mathtt{R}=d^{2}/(12K_{0})$.

Conservation of energy in the mushy layer is given by

$\overline{\rho c_{p}}\frac{\partial T}{\partial
t}+\rho_{b}c_{b}\boldsymbol{q}\cdot\boldsymbol{\nabla}T=\boldsymbol{\nabla}\cdot\left(\overline{k}\boldsymbol{\nabla}T\right)+\mathscr{L}M_{f},$
(22)

where $\overline{\rho c_{p}}$ and $\overline{k}$ are the solid-fraction-
weighted functions of ice and brine properties for heat capacity and thermal
conductivity, respectively. The latent heat is $\mathscr{L}$. We write a
conservation law for the concentration of brine, i.e.

$\left(1-\phi\right)\frac{\partial C}{\partial
t}+\boldsymbol{q}\cdot\boldsymbol{\nabla}C=\boldsymbol{\nabla}\cdot\left(\overline{D}\boldsymbol{\nabla}C\right)+\frac{\left(C-C_{s}\right)}{\rho}M_{f},$
(23)

where $\overline{D}$ is the phase-weighted salt diffusivity arising from
transport through the interstitial liquid.

The evolution of temperature and concentration are coupled in the mushy
region, as shown in figure 1 in the main text, where the ice and liquid phases
are at an equilibrium temperature given by the liquidus. This imposes a
constraint on the evolution equations (22) and (23) in mushy regions, where
the concentration and temperature are tied together by the liquidus
relationship $T=T_{L}(C)$. In idealized form, the liquidus relationship is
given as

$T_{L}=T_{m}-mC,$ (24)

where $T_{m}$ is the pure substance melting temperature and $m$ is the
liquidus slope. This constraint no longer applies when $T<T_{e}$, where
$M_{f}=0$, $\boldsymbol{q}=\boldsymbol{0}$, $\phi=1$, (22) becomes the heat
equation for thermal conduction and salt conservation (23) implies constant
bulk salinity.

##  4 Liquid volume and rates of production

In the main text, we show steady state temperature and porosity fields. We
treat the rate of shear heating as constant, given that the thermal diffusion
timescale is much longer than the tidal cycle (Nimmo and Gaidos, 2002; Nimmo
et al., 2007). Here, in figure 2(left), we show the approach to steady state
for the liquid volume in the mushy zone. Initially, there is no mushy zone and
therefore no liquid volume. Then a mushy zone forms and the liquid volume
increases rapidly. Due to the temperature field reaching a steady state and
the porosity-weakening of the fault friction saturates causing the heat flux
$F(1-\overline{\phi})$ to saturate in figure 2(left). In all simulations, the
local value of $\phi$ is used to calculate the shear heating rate. Here we
show the average value over fracture as an aggregate measure of the effect of
reducing the shear heating flux. From the curve of liquid volume with time, we
take the derivative to find the rate of liquid production, as shown in figure
2(right). This shows the same three regimes: no initial liquid volume, growth
of a mushy zone, and approach to steady state.

\begin{overpic}[width=138.76157pt]{pltU_tempcontour180.pdf}
\put(33.0,75.0){$F=180$}\end{overpic}\begin{overpic}[width=138.76157pt]{pltU_tempcontour240.pdf}
\put(33.0,75.0){$F=240$}\end{overpic}\begin{overpic}[width=138.76157pt]{pltU_tempcontour420.pdf}
\put(33.0,75.0){$F=420$}\end{overpic}\begin{overpic}[width=138.76157pt]{pltU_porosity180.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{pltU_porosity240.pdf}
\end{overpic}\begin{overpic}[width=138.76157pt]{pltU_porosity420.pdf}
\end{overpic}

Figure 3: Temperature and porosity plots for three different values of the
nondimensional heat flux: (left) $F=180$, (middle) $F=240$, and (right)
$F=420$. Note that porosity is shown with different axis than the temperature
plots. Similar to the main text, we show vertical velocity arrows with
logarithmic magnitude with the largest arrows representing about 1 mm/year.

The observed geyser ejection rate for all four tiger stripes is $\sim$0.2 m3
s-1 (Hansen et al., 2006, 2008). Following a step change in heating and for
the range of heat fluxes we used in these simulations, we find that the
maximum rate of liquid volume production is greater than $0.2$ m3 s-1
($F=500$), which is sufficient to explain the geyser flux from all of the
tiger stripe fractures. Figure 2 shows that the melt rate is equal to the
observed geyser ejection rate initially and stays at a comparable magnitude
for $\sim$200 thousand years. At later times, the melt rate reduces and melt
volume saturates because a region of large porosity develops and weakens the
frictional heating on the crack. However, venting of liquid from the crack
will limit porosity growth, and enable sustained heating and melting rates.

Figure 4: Flow velocity in the mushy zone for $F=500$: (left) vertical
velocity, showing buoyant fluid rising along the crack and (right) horizontal
velocity, demonstrating the buoyant plume is stable at $\mathtt{Ra_{C}}=200$
(Boury et al., 2021) with periodic structures in the outer region of the mushy
zone.

The fact that peak melt rates exceeds the observed geyser flux is an exciting
demonstration of the viability of the shear heating model for the genesis of
plumes on Enceladus and other icy satellites. While the simple configuration
with steady heating doesn't account for the inherent complexities of the
transient evolution, our impulse-response calculations likely elucidate the
order of magnitude and timescale of heating. In figure 3, we show the steady
state temperature and porosity for three increasing nondimensional flux
values, $F=180$, $F=240$, and $F=420$. These figures mirror figures 2(right)
and 3(right) in the main text, where we show the temperature and porosity
structure for $F=500$. For $F=260$, the maximum rate of liquid volume
production is on the order of $0.05$ m3 s-1, which is less than rate required
to indefinitely produce plumes with the observed ejection rate, but the mushy
zone liquid reservoir could sustain geysers for tens of thousands of years.
This suggests that there is a parallel to the Enceladus heat signal, which may
be an indication of oscillations in thermal activity and hence plume activity
(Spencer and Nimmo, 2013; Nimmo et al., 2018). In other words, the liquid
volume production rate by shear heating may allow for a large reservoir to
develop periodically, leading to geysers that transiently exhausts the supply.
Lastly, in our simulations, we reduce the shear heating rate as liquid is
generated in the vicinity of the fracture. If we were to remove the liquid as
geyser material, however, the shear heating could stay high, leading to
additional liquid volume production. We currently do not remove the liquid
from our simulation domain and leave these simulations to future work.

Within the mushy zone, there is convection of the pore fluid and the fluid
near the fracture rises buoyantly, as shown in the main text and in figures 4
and 5. The intensity of convection is controlled by the thermal and
compositional mushy Rayleigh numbers,

$\displaystyle\mathtt{Ra_{T}}=\frac{\rho_{w}g\alpha\Delta
TK_{0}h}{\kappa_{o}\eta},$ (25)
$\displaystyle\mathtt{Ra_{C}}=\frac{\rho_{w}g\beta\Delta
CK_{0}h}{\kappa_{o}\eta},$ (26)

respectively. Boury et al. (2021) show that the rising plume along the
fracture is unstable, with the development of the instability here resulting
in rapid variation in horizontal velocity leading to numerical convergence
issues in SOFTBALL (figure 5). To find solutions that converge in a reasonable
amount of time for our grid resolution, we follow Boury et al. (2021) and
restrict the Rayleigh numbers to values below the instability threshold,
$\mathtt{Ra_{C}}\approx 205$ for our system and use $\mathtt{Ra_{C}}=200$ for
the results we present. There is a large uncertainty in the permeability
prefactor $K_{0}$ and the Rayleigh numbers that we present use a value of
$K_{0}$ that is on the lower end for sea ice simulations (Polashenski et al.,
2017; Buffo et al., 2018; Boury et al., 2021). The low permeability limits
convection and allows the simulations to converge numerically. In figure 5, we
compare three increasing Rayleigh numbers
$\mathtt{Ra_{C}}=200,~{}400,~{}\mbox{and}~{}800$, where the latter two
simulations do not converge. The higher Rayleigh numbers show that the flow-
and-phase-change instability leads to fine-scale structures that ultimately
approach the grid resolution, affecting the pattern of the horizontal velocity
field. The porosity field, however, is similar across the three Rayleigh
numbers, suggesting that the results for total volume and peak melt rate would
approximately hold as the Rayleigh number increased.

\begin{overpic}[width=138.76157pt]{pltU_horizvel500_Ra200.pdf}\put(33.0,76.0){$\mathtt{Ra_{C}}=200$}\end{overpic}\begin{overpic}[width=138.76157pt]{pltU_horizvel500_Ra400.pdf}\put(33.0,76.0){$\mathtt{Ra_{C}}=400$}\put(22.0,11.0){not
converged}\end{overpic}\begin{overpic}[width=138.76157pt]{pltU_horizvel500_Ra800.pdf}\put(33.0,76.0){$\mathtt{Ra_{C}}=800$}\put(22.0,11.0){not
converged}\end{overpic}

Figure 5: Increasing instability with stronger convection for $F=500$: (left)
$\mathtt{Ra_{C}}=200$ showing a stable interior plume and an onset of
instability in the region around the liquid inclusion. (middle)
$\mathtt{Ra_{C}}=400$ showing an onset of instability in the interior plume.
(right) $\mathtt{Ra_{C}}=800$ showing an unstable interior plume and a stable
outer region that becomes unstable at the edge of the mushy zone. The three
porosity plots, however, are similar.

Near the fracture, shear heating melts the ice and dissolves the salt.
Therefore, the pore fluid in this region has a higher salinity (figure 6) than
the surrounding shell. Here we plot the bulk salinity $C$, which is given as

$C=\phi C_{\ell}+(1-\phi)C_{s},$ (27)

where $\phi$ is the porosity, $C_{\ell}$ is the liquid salinity, and $C_{s}$
is the solid salinity. Starting from a step change in heating, the evolution
of the bulk salinity is shown in figure 7. First, shear heating warms the
salty ice above the eutectic temperature generating a permeable zone. Since
the brine salinity and temperature are connected through the liquidus
constraint within the mushy zone, colder and saltier brine along the top then
sinks and the resulting convective desalination reduces the bulk salinity when
compared to the background value. Next, the sinking brine cannot escape and
collects near the base of the porous inclusion along the fracture, generating
a salt-rich zone with porosity near unity. Convection continues within the
mushy zone. The convective desalination leaves behind regions of high bulk
salinity in the outer parts of the porous inclusion, that are trapped due to
the reduction of permeability in the neighbouring desalinated ice which now
has low porosity. Due to the partial melting within the mushy zone, the
chemical composition of the geyser ejecta sourced from the shear heating model
that we propose here will be different from geyser ejecta sourced from the
ocean. In our simulations, the salinity in the fluid around the crack is
higher than the background shell salinity. The shell likely froze from the
underlying ocean and will have a salinity structure that varies with depth
(Buffo et al., 2021).

Figure 6: Salinity for the mushy zone with $F=500$. (left) Bulk salinity $C$,
showing that the fluid within the interstices is saltier than the surrounding
shell, which froze from the underlying ocean. (right) Liquid salinity
$C_{\ell}$ showing the increase in salinity with depth in the mushy zone.
High-salinity channels form at an oblique angle, due to the competition
between melting and gravity drainage (Buffo et al., 2023).

##  5 Dike propagation

In our steady state simulations, a region of salty brine develops along the
fracture with a porosity approaching a value of $\phi\approx 1$. Here we do
not allow the liquid to escape as geysers, therefore the liquid volume can
build up within the mushy zone. The question then arises as to whether the
brine pocket can hydraulically fracture through the ice shell and into the
underlying ocean, thereby generating a dike that opens a surface-to-base
connection. We analyze this possibility through a scaling analysis, following
Lister (1990) and Kalousová et al. (2016). For $F=500$, we take the half width
of brine region to be $a\sim 200$ m, the breadth into the page as $w\sim 500$
km, and the height as $\ell\sim 4$ km, we have two options for the balance of
gravity. Either, viscous stresses resist propagation or the fracture toughness
$K_{I}\sim 100$ kPa m1/2 resists opening. Starting with the case of viscous
stress, we have that

$\phi\Delta\rho ga\ell w\sim\mu\frac{V}{a^{2}}a\ell w,$ (28)

where $\phi\sim 0.25$ is a representative average porosity in the brine
pocket; $\Delta\rho\sim 90$ kg m-3 is the density difference between the brine
and the surrounding ice; $\mu\sim 10^{14}$ Pa$\cdot$s is the brine-saturated
temperate ice viscosity; and $V$ is the rate of dike propagation. Canceling
terms, we find that the rate of propagation scales as

$V\sim\frac{\phi\Delta\rho ga^{2}}{\mu}.$ (29)

The time $t_{\textrm{dike}}$ for a dike to propagate a distance $\ell$, then
scales as

$t_{\textrm{dike}}\sim\frac{\ell}{V}\sim\frac{\mu\ell}{\phi\Delta\rho
ga^{2}}\sim 110~{}\mbox{kyr}.$ (30)

This timescale is on the same order of magnitude as the time to steady state
in our shear heating calculations, indicating that the brine pocket could form
and then propagate as a dike within that time frame, if the porosity grows to
order 1 without significant venting of the liquid from the crack.

Figure 7: Bulk salinity evolution with $F=500$, showing the drainage and brine
concentration along the fracture: (left) 541 kyr. (middle) 865 kyr. (right)
13,000 kyr.

If the balance is instead between gravity and the fracture toughness, we have
that

$\phi\Delta\rho ga\ell w\sim\frac{K_{I}}{\ell^{1/2}}aw,$ (31)

which determines whether the fracture will propagate or not. Simplifying and
plugging in numbers, we have that

$\phi\Delta\rho g\ell^{3/2}\sim 700~{}\mbox{kPa
m${}^{1/2}$}~{}~{}>~{}~{}K_{I}\sim 100~{}\mbox{kPa m${}^{1/2}$}.$ (32)

Thus, we find that $\phi\Delta\rho g\ell^{3/2}$ is greater than the ice
fracture toughness $K_{I}$ and much larger than values measured in brine-
saturated sea ice, e.g. $\sim 50$ kPa m1/2 (Schulson and Duval, 2009),
demonstrating that the ice will not provide significant resistance to
fracture.

In both cases, we see that dike propagation is possible if the porosity
reaches order 1 over an extended region, before the liquid is vented. The time
scale for diking is on the same order of magnitude as the for our simulations
to reach a steady state and the ice fracture toughness is not large enough to
prevent fracture. These scaling calculations spawn a few ideas. An internal
dike is interesting because there could be periodic brine pocket generation,
leading to opening a connection between surface and ocean, then Kite and Rubin
(2016)-style eruptions. Or the dike could drain all of the fluid internally,
leading to a stronger fault with larger shear heating, and a periodic mushy
zone source for the geysers. However, by not allowing for interstitial liquid
to escape, our estimates for the liquid brine volume are artificially high,
meaning that a region with porosity approaching 1 may never form within the
mushy zone. In this case, there would be no prospect for a dike and the
material generated would only emanate out as geysers. We leave additional
analysis exploring these ideas to future work.

parameters |  |
---|---|---
specific heat of ice | $c_{i}$ | 2000 J kg-1 K-1
specific heat of ocean | $c_{o}$ | 4000 J kg-1 K-1
ice density | $\rho_{i}$ | 940 kg m-3
seawater density | $\rho_{w}$ | 1030 kg m-3
density coefficient for temperature | $\beta_{T}$ | 2.1$\times 10^{-4}$ kg m-3 K-1
density coefficient for salt | $\beta_{C}$ | 7.7$\times 10^{-4}$ kg m-3 ppt-1
thermal conductivity of ice | $k_{i}$ | 2.0 W m-1 K-1
thermal conductivity of ocean | $k_{o}$ | 0.6 W m-1 K-1
latent heat of fusion | $\mathscr{L}$ | 335000 J kg-1
water viscosity | $\eta$ | 2$\times 10^{-3}$ Pa s
salt diffusivity | $\kappa_{s}$ | 2$\times 10^{-9}$ m2 s-1
thermal diffusivity | $\kappa_{o}$ | 1.5$\times 10^{-7}$ m2 s-1
Stefan-Boltzmann constant | $\sigma$ |  $5.67\times 10^{-8}$ W m-2 K-4
shell salt composition | $C_{i}$ | 35 ppt
eutectic composition | $C_{e}$ | 230 ppt
liquidus slope | $m$ | -0.0913 K ppt-1
eutectic temperature | $T_{e}$ | 252 K
surface temperature | $T_{s}$ | 75 K
gravity on Enceladus | $g$ | 0.113 m s-2
fracture depth | $h$ | 5 km
height scale | $H$ | 10 km
permeability prefactor | $K_{0}$ | 3.2$\times 10^{-13}$ m2
Hele-Shaw cell spacing | $d$ | 5$\times 10^{-5}$ m
time scale | $t_{c}$ | 2.2$\times 10^{7}$ years
slip velocity | $u$ | 3.4$\times$10-6 m s-1
temperature scale | $\Delta T$ | 18 K
emissivity | $\epsilon$ | 1
coefficient of friction | $\mu$ | 0.3
Stefan number | $\mathtt{St}$ | 4.7
Prandtl number | $\mathtt{Pr}$ | 12
Lewis number | $\mathtt{Le}$ | 73
mushy compositional Rayleigh number | $\mathtt{Ra_{C}}$ | 200
mushy thermal Rayleigh number | $\mathtt{Ra_{T}}$ | 2
reluctance | $\mathtt{R}$ | 9.6
shear heating | $F$ | 500
linearized radiation, nondimensional $r_{1}$ | $G$ | 0
linearized radiation, nondimensional $r_{2}$ | $b$ | 100
linearized radiation, nondimensional $T_{s}$ | $\theta_{s}$ | -10
Table 1: Parameters used in the SOFTBALL numerical simulations, grouped by
constants, variables, nondimensional parameters. Definitions follow Parkinson
et al. (2020). Values from Nimmo and Gaidos (2002), Nimmo et al. (2007),
Spencer and Nimmo (2013) as well as Buffo et al. (2021).

## References

  * Abramov and Spencer (2009) O. Abramov and J. R. Spencer.  Endogenic heat from Enceladus' south polar fractures: New observations, and models of conductive surface heating.  _Icarus_ , 199(1):189 – 196, 2009.  doi: 10.1016/j.icarus.2008.07.016.
  * Boury et al. (2021) S. Boury, C. R. Meyer, G. M. Vasil, and A. J. Wells.  Convection in a mushy layer along a vertical heated wall.  _J. Fluid Mech._ , 926:A33, 2021.  doi: 10.1017/jfm.2021.742.
  * Broberg (1999) K. B. Broberg.  _Cracks and fracture_.  Academic Press, 1999.
  * Buffo et al. (2018) J. J. Buffo, B. E. Schmidt, and C. Huber.  Multiphase reactive transport and platelet ice accretion in the sea ice of McMurdo Sound, Antarctica.  _J. Geophys. Res._ , 123(1):324–345, 2018.  doi: 10.1002/2017JC013345.
  * Buffo et al. (2021) J. J. Buffo, C. R. Meyer, and J. R. G. Parkinson.  Dynamics of a solidifying icy satellite shell.  _J. Geophys. Res._ , 126(5):e2020JE006741, 2021\.  doi: 10.1029/2020JE006741.
  * Buffo et al. (2023) J. J. Buffo, C. R. Meyer, C. J. Chivers, C. C. Walker, C. Huber, and B. E. Schmidt.  Geometry of freezing impacts ice composition: Implications for icy satellites.  _J. Geophys. Res._ , 128(3):e2022JE007389, 2023\.  doi: 10.1029/2022JE007389.
  * Crawford and Stevenson (1988) G. D. Crawford and D. J. Stevenson.  Gas-driven water volcanism in the resurfacing of Europa.  _Icarus_ , 73:66–79, 1988.  doi: 10.1016/0019-1035(88)90085-1.
  * Hammond (2020) N. P. Hammond.  Estimating the magnitude of cyclic slip on strike-slip faults on Europa.  _J. Geophys. Res._ , 125(7):e2019JE006170, 2020\.  doi: 10.1029/2019JE006170.
  * Hansen et al. (2006) C. J. Hansen, L. Esposito, A. I. F Stewart, J. Colwell, A. Hendrix, W. Pryor, D. Shemansky, and R. West.  Enceladus' water vapor plume.  _Science_ , 311(5766):1422–1425, 2006.  doi: 10.1126/science.1121254.
  * Hansen et al. (2008) C. J. Hansen, L. W. Esposito, A. I. F. Stewart, B. Meinke, B. Wallis, J. E. Colwell, A. R. Hendrix, K. Larsen, W. Pryor, and F. Tian.  Water vapour jets inside the plume of gas leaving enceladus.  _Nature_ , 456(7221):477–479, 2008.  doi: 10.1038/nature07542.
  * Hemingway et al. (2018) D. Hemingway, L. Iess, R. Tajeddine, and G. Tobie.  The interior of Enceladus.  In P. M. Schenk, R. N. Clark, C. J. A. Howett, A. J. Verbiscer, and J. H. Waite, editors, _Enceladus and the Icy Moons of Saturn_ , pages 57–78. University of Arizona Press, 2018.  ISBN 9780816537075.  doi: 10.2307/j.ctv65sw2b.12.
  * Hemingway et al. (2020) D. J. Hemingway, M. L. Rudolph, and M. Manga.  Cascading parallel fractures on Enceladus.  _Nat. Astron._ , 4:234–239, 2020.  doi: 10.1038/s41550-019-0958-x.
  * Hurford et al. (2007) T. A. Hurford, P. Helfenstein, G. V. Hoppa, R. Greenberg, and B. G. Bills.  Eruptions arising from tidally controlled periodic openings of rifts on Enceladus.  _Nature_ , 447:292–294, 2007.  doi: 10.1038/nature05821.
  * Kalousová et al. (2016) K. Kalousová, O. Souček, G. Tobie, G. Choblet, and O. Čadek.  Water generation and transport below Europa's strike-slip faults.  _J. Geophys. Res._ , 121:2444–2462, 2016.  doi: 10.1002/2016JE005188.
  * Katz and Worster (2008) R. F. Katz and M. G. Worster.  Simulation of directional solidification, thermochemical convection, and chimney formation in a Hele-Shaw cell.  _J. Comput. Phys._ , 227(23):9823–9840, 2008\.  doi: 10.1016/j.jcp.2008.06.039.
  * Kite and Rubin (2016) E. S. Kite and A. M. Rubin.  Sustained eruptions on Enceladus explained by turbulent dissipation in tiger stripes.  _Proc. Natl. Acad. Sci._ , 113(15):3972–3975, 2016.  doi: 10.1073/pnas.1520507113.
  * Lister (1990) J. R. Lister.  Buoyancy-driven fluid fracture: similarity solutions for the horizontal and vertical propagation of fluid-filled cracks.  _J. Fluid Mech._ , 217:213–239, 1990.  doi: 10.1017/S0022112090000696.
  * Nimmo (2004) F. Nimmo.  Stresses generated in cooling viscoelastic ice shells: Application to Europa.  _J. Geophys. Res._ , 109(E12), 2004.  doi: 10.1029/2004JE002347.
  * Nimmo and Gaidos (2002) F. Nimmo and E. Gaidos.  Strike-slip motion and double ridge formation on Europa.  _J. Geophys. Res._ , 107(E4):5–1–5–8, 2002\.  doi: 10.1029/2000JE001476.
  * Nimmo et al. (2007) F. Nimmo, J. R. Spencer, R. T. Pappalardo, and M. E. Mullen.  Shear heating as the origin of the plumes and heat flux on Enceladus.  _Nature_ , 447(7142):289, 2007.  doi: 10.1038/nature05783.
  * Nimmo et al. (2011) F. Nimmo, B. G. Bills, and P. C. Thomas.  Geophysical implications of the long-wavelength topography of the Saturnian satellites.  _J. Geophys. Res._ , 116(E11), 2011.  doi: 10.1029/2011JE003835.
  * Nimmo et al. (2018) F. Nimmo, A. C. Barr, M. Běhounková, and W. B. McKinnon.  The thermal and orbital evolution of Enceladus: Observational constraints and models.  In P. M. Schenk, R. N. Clark, C. J. A. Howett, A. J. Verbiscer, and J. H. Waite, editors, _Enceladus and the Icy Moons of Saturn_ , pages 79–94. University of Arizona Press, 2018.  ISBN 9780816537075.  doi: 10.2307/j.ctv65sw2b.13.
  * Nye (1955) J. F. Nye.  Comments on Dr. Loewe's letter and notes on crevasses.  _J. Glaciol._ , 2(17):512–514, 1955.  doi: 10.3189/S0022143000032652.
  * Olgin et al. (2011) J. G. Olgin, B. R. Smith-Konter, and R. T. Pappalardo.  Limits of Enceladus's ice shell thickness from tidally driven tiger stripe shear failure.  _Geophys. Res. Lett._ , 38:L02201, 2011.  doi: 10.1029/2010GL044950.
  * Parkinson et al. (2020) J. R. G. Parkinson, D. F. Martin, A. J. Wells, and R. F. Katz.  Modelling binary alloy solidification with adaptive mesh refinement.  _J. Comput. Phys._ , 5:100043, 2020.  doi: 10.1016/j.jcpx.2019.100043.
  * Polashenski et al. (2017) C. Polashenski, K. M. Golden, D. K. Perovich, E. Skyllingstad, Alexandra Arnsten, C. Stwertka, and N. Wright.  Percolation blockage: A process that enables melt pond formation on first year arctic sea ice.  _J. Geophys. Res._ , 122(1):413–440, 2017.  doi: 10.1002/2016JC011994.
  * Rudolph et al. (2022) M. L. Rudolph, M. Manga, M. Walker, and A. R. Rhoden.  Cooling crusts create concomitant cryovolcanic cracks.  _Geophys. Res. Lett._ , 49(5):e2021GL094421, 2022\.  doi: 10.1029/2021GL094421.
  * Schulson and Duval (2009) E. M. Schulson and P. Duval.  _Creep and fracture of ice_.  Cambridge University Press, 2009.
  * Shibley and Laughlin (2021) N. C. Shibley and G. Laughlin.  Do oceanic convection and clathrate dissociation drive Europa’s geysers?  _Planet. Sci. J._ , 2(6):221, 2021.  doi: 10.1029/2022GL098621.
  * Smith (1976) R. A. Smith.  The application of fracture mechanics to the problem of crevasse penetration.  _J. Glaciol._ , 17(76):223–228, 1976.  doi: 10.3198/1976JoG17-76-223-228.
  * Smith-Konter and Pappalardo (2008) B. Smith-Konter and R. T. Pappalardo.  Tidally driven stress accumulation and shear failure of Enceladus's tiger stripes.  _Icarus_ , 198(2):435–451, 2008.  doi: 10.1016/j.icarus.2008.07.005.
  * Spencer and Nimmo (2013) J. R. Spencer and F. Nimmo.  Enceladus: An active ice world in the Saturn system.  _Annu. Rev. Earth Planet. Sci._ , 41(1):693–717, 2013.  doi: 10.1146/annurev-earth-050212-124025.
  * Stevenson (1996) D. J. Stevenson.  Heterogeneous tidal deformation and geysers on Europa.  In _Europa Ocean Conference_ , number 5 in Capistrano conference, pages 69–70, San Juan Capistrano, Calif., 1996. San Juan Capistrano Res. Inst.
  * Tait and Jaupart (1992) S. Tait and C. Jaupart.  Compositional convection in a reactive crystalline mush and melt differentiation.  _J. Geophys. Res._ , 97(B5):6735–6756, 1992.  doi: 10.1029/92JB00016.
  * Timco and Frederking (1996) G. W. Timco and R. M. W. Frederking.  A review of sea ice density.  _Cold Reg. Sci. Technol._ , 24(1):1–6, 1996.  ISSN 0165-232X.  doi: 10.1016/0165-232X(95)00007-X.
  * Walker et al. (2021) C. C. Walker, J. N. Bassis, and B. E. Schmidt.  Propagation of vertical fractures through planetary ice shells: The role of basal fractures at the ice–ocean interface and proximal cracks.  _Planet. Sci. J._ , 2(4):135, 2021.  doi: 10.3847/PSJ/ac01ee.
  * Wells et al. (2019) A. J. Wells, J. R. Hitchen, and J. R. G. Parkinson.  Mushy-layer growth and convection, with application to sea ice.  _Philos. Trans. R. Soc. London, Ser. A_ , 377(2146):20180165, 2019.  doi: 10.1098/rsta.2018.0165.