# Pattern Formation of Freezing Infiltration in Porous Media

Nathan D. Jones<EMAIL_ADDRESS>Adrian Moure<EMAIL_ADDRESS>Xiaojing
Fu<EMAIL_ADDRESS>Department of Mechanical and Civil Engineering,
California Institute of Technology, 1200 E California Blvd, Pasadena CA, 91125

###### Abstract

Gravity-driven infiltration of liquid water into unsaturated porous media can
be a spatially heterogeneous process due to the gravity fingering instability.
When such infiltration occurs in a subfreezing porous medium, liquid water can
readily freeze, leading to both the removal of liquid water available for
transport and a reduction in local permeability. As a result of the coupling
between gravity fingering and freezing, macroscopic frozen structures can form
that record the shape and history of the wetting front. These structures have
been observed in the field in terrestrial snowpack and glacial firn layers and
are believed to have profound impacts on how liquid water and its accompanying
thermal content distribute during infiltration. However, a more detailed
physics-based understanding of freezing infiltration has been missing. In this
work, we use a thermodynamic nonequilibrium infiltration model to investigate
the emergence of refrozen structures during water infiltration into an
initially homogeneous and subfreezing porous medium. From scaling analysis, we
recover the relevant nondimensional groups that govern the physics of the
freezing infiltration process. We identify two key mechanisms caused by
freezing that reduce the effective infiltration rate, calculated as the
maximum depth of infiltration per elapsed time. In the first mechanism, the
effective infiltrate rate decreases because a portion of the liquid water is
consumed due to freezing, and such effect can be well quantified by the
freezing Damköhler number. For the second mechanism, we report on a new
phenomenon termed secondary fingering, where new flow paths are established in
between the primary infiltration channels. We find that secondary fingering
reduces the degree of flow channelization and thus weakens the effective rate
of infiltration via flow field homogenization.

##  I Introduction

When a viscous fluid infiltrates into an initially dry porous medium under
gravity, an instability can readily occur that results in the formation of
preferential flow pathways. In the absence of macro pores or other types of
heterogeneities in the porous medium [1, 2], pore-scale competition amongst
gravitational, viscous, and capillary forces at the wetting front gives rise
to this instability. This phenomenon, often referred to as the gravity
fingering instability, is widely observed in the field [3, 4, 5] and studied
in laboratory experiments [6, 7, 8, 9] for various porous media such as sand,
soil, or snow. Preferential flow in unsaturated media will enhance the
effective infiltration rate of the wetting-phase fluid compared to stable flow
[10], and has shown to accelerate the arrival of solutes or contaminants [11,
12] into deeper regions of the domain.

Gravity fingering has recently gathered a renewed interest in the field of
snow and firn hydrology. In this context, water and its associated thermal
content enter the snowpack from radiation-induced surface melting or rain-on-
snow events [13, 14] and actively modify the thermal budget and permeability
structure of the porous ice. Experiments of gravity-driven water percolation
both in the field, as shown in Fig. 1(a), and in the laboratory in 3D snow
columns have demonstrated that the process is intrinsically unstable, giving
rise to preferential pathways [15, 16, 17]. Meanwhile, as snow exists below
the freezing point of water, heat transfer and phase change processes readily
occur during the infiltration process. In contrast to isothermal gravity
fingering, here the effect of phase change will hinder the vertical
infiltration rate of meltwater by converting the amount of liquid water
available for transport into solid ice and consequently lowering the local
permeability in the snow. This refreezing of melt within the snowpack forms
cylindrical frozen structures known as ice pipes on the order of centimeters
in width at the site of preferential flow pathways [18, 19, 20, 21].
Additionally, during contact with discrete snow layers below the surface,
vertical flow can be stopped and redirected laterally, forming horizontal
frozen structures commonly referred to as ice lenses [5]. Evidence of these
refrozen structures has been widely observed in the field [22, 23], both from
excavated snow pits [Fig. 1(b)] as well as radar detection [24, 25] for
larger-scale frozen structures such as ice lenses or slabs. The formation of
refrozen structures is known to affect the thermal budget of the snowpack due
to latent heat release [26, 27], and the lasting structures affect the bulk
hydraulic, thermodynamic, and structural properties of the snowpack.

Figure 1:  Preferential flow during infiltration of water into subfreezing
porous ice: (a) Photograph of a red dye tracer experiment from Ref. [10] a
snow pit exhibiting preferential flow. Photograph courtesy of Philip Marsh.
(b) Field image from Ref. [20] of excavated refrozen structures created after
preferential infiltration in snow. (c) Problem setup for simulations presented
in this study. Liquid water at $T_{\mathrm{melt}}=0~{}^{\circ}$C infiltrates
through the top of the domain and interacts with the subfreezing porous ice.

Motivated by these observations, there have been a series of experimental and
modeling efforts to better understand the infiltration process of water
through permeable ice such as snow or firn. Experimental studies that
investigate gravity-driven preferential flow in porous ice typically operate
at isothermal conditions [17, 28, 5] in which no phase change is permitted to
occur. Additionally, due to the opacity of porous ice, directly imaging the
macroscopic infiltration process is difficult, and requires the dissection of
experiments a posteriori to investigate the flow patterns [28] or, more
recently, the employment of non-intrusive techniques such as neutron
radiography [29] or magnetic resonance imaging [30, 31]. Many numerical
studies have been performed that investigate the formation of refrozen
structures (in particular, ice lenses) using one-dimensional models [32, 33,
34, 35] without considering the flow instability. Preferential flow in porous
ice has also been investigated through numerical simulations, but these
studies mainly consider isothermal conditions [36, 37, 38, 39] or are limited
by the assumption of local thermodynamic equilibrium [40, 41]. Despite the
efforts mentioned above, a systematic investigation of the mechanical details
of infiltration into subfreezing porous media and the resulting mass transfer
dynamics is still lacking.

In this work, we investigate gravity-driven unstable infiltration of liquid
water into an initially dry porous medium composed of ice and air. We use a
recently developed model [42] that couples a two-phase flow model of unstable
infiltration with nonequilibrium thermodynamics of ice-water heat transfer and
phase change. In this model, water flow is prescribed with a Darcy velocity
that takes contributions from gravity, capillary pressure, and macroscopic
surface tension effects. The model captures phase change between liquid water
and ice by assigning independent temperature fields to the water and ice
phases and thus allowing local thermodynamic non-equilibrium (LTNE). The LTNE
assumption allows us to impose phase change rates based on thermodynamic
properties and local temperatures of the porous media instead of fixed rates
of phase change. We perform high-resolution numerical simulations of the model
on a 2D rectilinear domain and focus our analysis on the nonlinear regime past
the onset of the fingering instability. Our simulations capture the formation
of vertical refrozen structures during infiltration over a wide range of flow
rates and thermal conditions. Quantitative analysis of the simulations shows
that the effective infiltration rate coupled with phase change can be well
predicted by the freezing Damköhler number. We also show that increased rates
of both freezing and water flux can decrease the heterogeneity of the flow
field by forming additional preferential pathways. We report on a new
phenomenon we term as “secondary fingering”, where new infiltration channels
form between pre-established ones. This phenomenon occurs only under the
regime of high rates of both freezing and flow, and arises primarily due to
the permeability reduction from freezing. The secondary fingering process is
in contrast to classical gravity-driven infiltration, in which the
infiltrating fluid tends to travel through pre-existing pathways.

##  II Governing Equations

We consider a porous matrix of ice, water, and air, in which the volume-
averaged quantities water saturation $S(\mathbf{x},t)$ and porosity
$\phi(\mathbf{x},t)$ determine the phase composition throughout the domain.
The volumetric liquid water content (LWC) is defined as $\mathrm{LWC}=\phi S$
and the solid fraction of ice is $1-\phi$. Flow of water is governed by a
Darcy velocity $\bm{u}(\mathbf{x},t)$. The temperature of the water and ice
phases are $T_{w}(\mathbf{x},t)$ and $T_{i}(\mathbf{x},t)$, respectively. The
full details of this model can be found in Ref. [42]. Here, we provide a brief
summary of the governing equations and their underlying assumptions.

The conservation equations for ice and water mass take the form:

$\displaystyle\frac{\partial(1-\phi)}{\partial
t}=-R_{m}W_{\mathrm{SSA}}\left(T_{\mathrm{int}}-T_{\mathrm{melt}}\right),$ (1)
$\displaystyle\frac{\partial(\phi S)}{\partial
t}+\nabla\cdot\bm{u}=\frac{\rho_{i}}{\rho_{w}}R_{m}W_{\mathrm{SSA}}(T_{\mathrm{int}}-T_{\mathrm{melt}}),$
(2)

where $\rho_{i}$ and $\rho_{w}$ are the densities of ice and water,
respectively, $T_{\mathrm{int}}$ is the temperature of the ice-water
interface, and $T_{\mathrm{melt}}=0$ ∘C is the melting point of ice. The phase
change rate coefficient $R_{m}$ is defined as
$R_{m}=c_{p,w}/(L_{\mathrm{sol}}\beta_{\mathrm{sol}})$, where $c_{p,w}$ is the
specific heat capacity of liquid water, $L_{\mathrm{sol}}$ is the latent heat
of solidification of ice, and $\beta_{\mathrm{sol}}$ is the kinetic attachment
coefficient for ice growth from liquid water [43]. The wet specific surface
area $W_{\mathrm{SSA}}$ represents the ice-water interfacial area per unit
volume:

$W_{\mathrm{SSA}}=S\frac{\mathrm{SSA}_{0}}{\phi_{0}\mathrm{ln}(\phi_{0})}\phi\mathrm{ln}(\phi),$
(3)

where $\mathrm{SSA}_{0}$ is the initial specific surface area. The interfacial
temperature $T_{\mathrm{int}}$ is an averaged pore-scale quantity that is
defined from the volume integrated Stefan condition:

$\frac{T_{\mathrm{int}}-T_{\mathrm{melt}}}{L_{\mathrm{sol}}/c_{p,w}}=-\beta_{\mathrm{sol}}\overline{v_{n}},$
(4)

where $\overline{v_{n}}$ is the volume averaged velocity of the fluid-solid
interface and is taken to be positive for freezing. Following Moure  _et al._
[42], the interface temperature takes the form:

$T_{\mathrm{int}}(T_{i},T_{w})=\frac{\frac{k_{i}}{r_{i}}T_{i}+\frac{k_{w}}{r_{w}}T_{w}}{\frac{c_{p,w}\rho_{w}}{\beta_{\mathrm{sol}}}+\frac{k_{i}}{r_{i}}+\frac{k_{w}}{r_{w}}},$
(5)

where $k_{i}$ and $k_{w}$ are the thermal conductivity of the ice and water
phases, and $r_{i}$ and $r_{w}$ are thermal length scales of the ice and water
phases, respectively. The parameters $r_{i}$ and $r_{w}$ represent thermal
diffusion lengths that are determined from pore-scale simulations of the
generalized Stefan problem for solidification (further details in Moure  _et
al._ [42]). Note that the right-hand side expressions in Eqs. (1) and (2)
differ by a factor of $\rho_{i}/\rho_{w}$ and have opposite signs, which
ensures mass conservation during phase change.

The model assumes local thermodynamic nonequilibrium by tracking the
temperatures of the ice and water phases independently:

$\displaystyle\frac{\partial[(1-\phi)T_{i}]}{\partial
t}=\nabla\cdot\left(D_{i}(1-\phi)\nabla
T_{i}\right)+W_{\mathrm{SSA}}D_{i}\frac{T_{\mathrm{int}}-T_{i}}{r_{i}},$ (6)
$\displaystyle\frac{\partial(\phi ST_{w})}{\partial
t}+\nabla\cdot\left(\bm{u}T_{w}\right)=\nabla\cdot\left(D_{w}\phi S\nabla
T_{w}\right)+W_{\mathrm{SSA}}D_{w}\frac{T_{\mathrm{int}}-T_{w}}{r_{w}},$ (7)

where and $D_{i}$ and $D_{w}$ are the thermal diffusivities of ice and water,
respectively.

To capture unstable gravity-driven infiltration in porous media, here we adopt
the unsaturated flow model proposed by Cueto-Felgueroso and Juanes [44, 45,
46]:

$\bm{u}=-K_{s}(\phi)k_{r}(S)\nabla\Pi(S),$ (8)

where $K_{s}(\phi)$ is the saturated hydraulic conductivity, $k_{r}(S)$ is the
relative permeability, and $\Pi(S)$ is the total flow potential. We use an
empirical expression for the snow hydraulic conductivity [47]:

$K_{s}(\phi)=3\left(\frac{d_{i}}{2}\right)^{2}\frac{\rho_{w}g}{\mu_{w}}\mathrm{exp}\left[-0.013\rho_{i}(1-\phi)\right],$
(9)

where $d_{i}$ is the average ice grain diameter, $g$ is the gravitational
acceleration, and $\mu_{w}$ is the dynamic viscosity of water, all of which we
assume to be constant. The relative permeability takes the form of a convex
function of saturation [48]:

$k_{r}(S)=\left(\frac{S-S_{r}}{1-S_{r}}\right)^{a},$ (10)

where $S_{r}$ is the irreducible saturation and the exponent $a>1$ is
dependent on the porous media. Here, we take $S_{r}=1\times 10^{-3}$ and
$a=5$. The expression for the flow potential is [44, 45]:

$\Pi(S)=y-\psi(S)-\sqrt{\kappa}\nabla\cdot(\sqrt{\kappa}\nabla S),$ (11)

where $y$ is the vertical coordinate (increasing with height), $\psi(S)$ is
capillary suction derived from Leverett scaling, and $\kappa(S)$ is the
capillary pinning function:

$\displaystyle\psi(S)=h_{\mathrm{cap}}S^{-\frac{1}{\alpha}}\left\\{1-\mathrm{exp}\left[\beta(S-\nu_{e})\right]\left(1+\beta\frac{\alpha}{\alpha-1}S\right)\right\\},$
(12)
$\displaystyle\kappa(S)=h_{\mathrm{cap}}^{2}\int_{0}^{S}\psi(S)\,\mathrm{d}S=h_{\mathrm{cap}}^{3}\frac{\alpha}{\alpha-1}S^{\frac{\alpha-1}{\alpha}}\left\\{1-\mathrm{exp}\left[\beta(S-\nu_{e})\right]\right\\},$
(13)

where $h_{\mathrm{cap}}$, $\alpha$, $\beta$, and $\nu_{e}$ are constants that
are calibrated from water retention curves for a given porous material [49,
50, 17]. These functions along with the values for these parameters are shown
in Fig. 2. Note that in the case of $\kappa=0$, the classical Richards
equation for unsaturated flow in porous media is recovered.

Figure 2:  Capillary pressure $\psi$ and the capillary pinning function
$\kappa$ used in this work. Here, we plot Eq. (12) and (13) with
$h_{\mathrm{cap}}=2.5$ cm, $\alpha=4$, $\beta=22$, and $\nu_{e}=1$. Parameters
are calibrated from experimental measurements from Ref. [17].

##  III Problem Setup and Scaling analysis

We apply the above model to the problem of uniform infiltration in 2D, in
which we introduce liquid water into a subfreezing 2D porous medium through
the top boundary at a constant flow rate $U$, as illustrated in Fig. 3. The
incoming water temperature is prescribed as $T_{\mathrm{melt}}$, and the
initial temperature of the porous ice is $T_{i,0}$. A no-flux boundary
condition is imposed on the lateral boundaries and bottom boundary of the
domain for both flow and thermal energy. We identify the following
characteristic scales in our system: $|\Delta
T_{0}|=|T_{i,0}-T_{\mathrm{melt}}|$ is the characteristic temperature, $U$ is
the characteristic velocity, $L=\mathrm{SSA}_{0}^{-1}$ is the length scale,
and hence the product $t_{c}=(\mathrm{SSA}_{0}U)^{-1}$ is the characteristic
time scale. We use $\mathrm{SSA}_{0}^{-1}$ as the length scale as it provides
a more appropriate volume-averaged measurement of the average grain size. Note
that $\mathrm{SSA}_{0}^{-1}\sim d_{i}$ for spherical ice grains. The
nondimensional equations for mass and thermal conservation of each phase then
take the following form:

$\displaystyle\frac{\partial\phi}{\partial t}$
$\displaystyle=\mathrm{Da}W_{\mathrm{SSA}}(T_{\mathrm{int}}-T_{\mathrm{melt}}),$
(14) $\displaystyle\frac{\partial(\phi S)}{\partial t}+\nabla\cdot\bm{u}$
$\displaystyle=\frac{\rho_{i}}{\rho_{w}}\mathrm{Da}W_{\mathrm{SSA}}(T_{\mathrm{int}}-T_{\mathrm{melt}}),$
(15) $\displaystyle\frac{\partial[(1-\phi)T_{i}]}{\partial t}$
$\displaystyle=\mathrm{Pe}_{i}^{-1}\left[\nabla\cdot\left((1-\phi)\nabla
T_{i}\right)+W_{\mathrm{SSA}}(T_{\mathrm{int}}-T_{i})\right],$ (16)
$\displaystyle\frac{\partial(\phi ST_{w})}{\partial
t}+\nabla\cdot\left(\bm{u}T_{w}\right)$
$\displaystyle=\mathrm{Pe}_{w}^{-1}\left[\nabla\cdot\left(\phi S\nabla
T_{w}\right)+W_{\mathrm{SSA}}(T_{\mathrm{int}}-T_{w})\right],$ (17)

Figure 3:  Snapshots of problem unknowns at $t=200$ min after continuous
infiltration for: (a) saturation $S$, (b) porosity $\phi$, (c) ice temperature
$T_{i}$, and (d) water temperature $T_{w}$. For this simulation, $U=27$ mm/h,
$T_{i,0}=-10$ ºC, $L_{x}=1$ m, and $L_{y}=2$ m. Porosity is decreased locally
during freezing, leaving behind vertical frozen structures that show the
history of the gravity fingers. Thermal equilibrium in the ice phase
temperature $T_{i}$ is established behind the fingering front, and phase
change does not continue to occur in this region. The water temperature
$T_{w}$ also reaches thermal equilibrium behind the wetting front and exhibits
local undercooling at the fingering front where phase change is occurring.

The Damköhler number Da sets the ratio between the rate of phase change and
the rate of advection:

$\mathrm{Da}=\frac{c_{p,w}|\Delta
T_{0}|}{L_{\mathrm{sol}}U\beta_{\mathrm{sol}}}.$ (18)

Applying Eq. (4) to the above definition allows us to rewrite Da as the ratio
between the freezing velocity at the ice-water interface, $\overline{v_{n}}$,
and the rate of advection:

$\mathrm{Da}=\frac{\mathrm{Ste}}{U\beta_{\mathrm{sol}}}=\frac{\overline{v_{n}}}{U}.$
(19)

Here, $\mathrm{Ste}$ is the Stefan number, which sets the ratio between the
amount of sensible heat to the latent heat of solidification:

$\mathrm{Ste}=\frac{c_{p,w}|\Delta T_{0}|}{L_{\mathrm{sol}}}.$ (20)

As Da increases, the rate of freezing becomes comparable to the rate of
infiltration. The thermal Péclet numbers associated with ice and water set the
ratio between advection and thermal diffusion:

$\displaystyle\mathrm{Pe}_{i}=UL/D_{i},$
$\displaystyle\mathrm{Pe}_{w}=UL/D_{w}.$ (21)

The thermal Péclet numbers control how quickly heat is transferred by
diffusion both ahead of the infiltration front as well as in the spanwise
direction. When Pe is sufficiently small, thermal diffusion dominates and any
thermal nonequilibrium is due to the kinetic limitations of phase change. As
Pei and Pew differ only by the thermal diffusivity of each phase, we refer to
$\mathrm{Pe}=\mathrm{Pe}_{i}$ in this work.

The above scaling analysis introduces the nondimensional groups that control
the heat transfer and phase change dynamics of the problem. The analysis
assumes the natural length scale of the porous media
($\mathrm{SSA}_{0}^{-1}$), which correlates with the grain size. However, the
original isothermal problem of gravity-driven infiltration bears its own
scaling analysis that yields additional nondimensional groups that control the
length scale and onset of the fingering pattern. Applying the same scaling
analysis to Eqs. (8) and (11), we arrive at the nondimensional equations for
velocity and flow potential:

$\displaystyle\bm{u}=-\frac{K^{\prime}_{s}(\phi)}{R_{s}}k_{r}\nabla\Pi^{\prime},$
(22)
$\displaystyle\Pi^{\prime}=y^{\prime}-N_{\mathrm{gr}}^{-1}\psi^{\prime}(S)-N_{\mathrm{nl}}\left(\sqrt{\kappa^{\prime}}\nabla\cdot(\sqrt{\kappa^{\prime}}\nabla
S)\right),$ (23)

where $\psi^{\prime}(S)=\psi(S)/h_{\mathrm{cap}}$,
$\kappa^{\prime}=\kappa/h_{\mathrm{cap}}^{3}$, and $N_{\mathrm{gr}}$ and
$N_{\mathrm{nl}}$ are the gravity number and nonlocal number, respectively.
Note that in the absence of phase change, $K_{s}$ is constant, and $U/K_{s}$
is defined as the flux ratio $R_{s}$ [6, 51, 44]. The influence of these two
additional nondimensional groups on the stability of the flow pattern is
investigated in Ref. [52], and they are defined:

$\displaystyle N_{\mathrm{gr}}=\frac{L}{h_{\mathrm{cap}}},$ $\displaystyle
N_{\mathrm{nl}}=\frac{\Gamma}{\rho gL^{3}},$ (24)

where $\Gamma$ is a coefficient in the nonlocal free energy potential that is
associated with an effective surface tension of the wetting front interface
[53, 45]. Following Ref. [46], we take $\Gamma=\rho g\kappa$, and recognizing
that $\kappa\sim h^{3}_{\mathrm{cap}}$ from Eq. (13), we arrive at the
relation:

$N_{\mathrm{nl}}=N_{\mathrm{gr}}^{-3}.$ (25)

Increasing $N_{\mathrm{gr}}$ has been shown to decrease the intrinsic length
scale of the gravitational fingering instability, leading to thinner fingers.
In addition to the hydraulic properties of the porous media such as the
average grain size diameter $d_{i}$ or the capillary height $h_{cap}$, both
the flux ratio $R_{s}$ and the initial saturation $S_{0}$ of the porous media
have been shown to influence the onset and nonlinear regimes of the pattern
formation process [54, 44, 52].

As we focus on the effects of the supplied flux $U$ and the rate of freezing
on the pattern formation process, in this work we do not vary
$N_{\mathrm{gr}}$ nor $N_{\mathrm{nl}}$. Although the entry flux $U$ has been
shown to influence the size and spacing of fingers in an isothermal scenario,
in this work, we explore a more narrow range of flux ratios between $1.2\times
10^{-4}$ to $5.3\times 10^{-4}$ ($U$ between 6.3 to 27 mm/h), and thus the
variability of the finger width and spacing is negligible. Larger input fluxes
may result in significantly wider fingers in these simulations, however here
we chose the range of $U$ to reflect appropriate values of surface melt rates
for natural snowpack reported in the literature [55].

##  IV Simulation Setup and Numerical Methods

We conduct high-resolution numerical simulations of the model to investigate
the dynamics and patterns for flow and porosity structures during water
infiltration into subfreezing porous ice. The simulations are performed on a
2D rectilinear domain of 1 m $\times$ 2 m ($L_{x}\times L_{y}$) as shown in
Fig. 1(c) and Fig. 3. Initially, the domain is composed of dry porous ice with
a uniform initial ice temperature $T_{i,0}\leq 0$ ∘C, uniform porosity
$\phi(\mathbf{x},0)=0.456$, and average ice grain diameter $d_{i}=1.5$ mm
which mimics the hydraulic properties of coarse snow. We assume the porous
media is initially at the irreducible saturation
$S(\mathbf{x},0)=S_{r}=1\times 10^{-3}$. The initial specific surface area is
estimated from Ref. [56] to be $\mathrm{SSA}_{0}=3514$ mm-1. We approximate
the kinetic attachment coefficient for ice growth from Ref. [43] as
$\beta_{\mathrm{sol}}=800$ s/m.

To explore the effect of the relative strength of water freezing rate on the
overall infiltration pattern, we systematically vary the surface meltwater
flux $U$ and the initial ice phase temperature $T_{i,0}$. We vary $U$ between
6.3–27 mm/h to explore flow rates near values reported in Ref. [55] for
natural systems and consider initial ice temperatures between 0 to -20 ∘C. The
domain has $400\times 800$ ($N_{x}\times N_{y}$) mesh elements, with the
exception of simulations presented in Sec. V.2, in which we use a finer mesh
resolution. Equations (1), (2), (6), (7), and (8) are solved using
isogeometric analysis [57]—a finite element method that implements B-splines
as basis functions. The time advancement scheme is a generalized–$\alpha$
method with an adaptive time step based on the number of Newton-Raphson
iterations. To implement the above methods, we develop a code using the open
source libraries PETSc [58] and PetIGA [59].

##  V Results

###  V.1 Impact of freezing on infiltration depth

Freezing can directly hinder gravity-driven infiltration through two
mechanisms. In the first mechanism, mass in the liquid phase is lost as it is
converted into ice at the site of freezing and thus the amount of total water
available for transport is reduced. In the second mechanism, an increase in
the ice volume fraction leads to a reduction in local permeability which
lowers the Darcy velocity of infiltration. Both of these effects will slow the
effective rate of infiltration of the wetting phase into the porous medium
compared to the isothermal problem without freezing. In Fig. 4, we demonstrate
this effect by plotting the liquid water content LWC and porosity $\phi$ at a
fixed time ($t_{\mathrm{f}}=212$ min) after constant infiltration of $U=27$
mm/h for varying $T_{i,0}$. With decreasing $T_{i,0}$, the wetting front
extends shallower into the domain; at the same time, the local porosity
reduction due to freezing becomes more pronounced.

Figure 4:  Snapshots of liquid water content (top) and porosity (bottom) at
$t_{\mathrm{f}}=212$ min for varying $T_{i,0}$. All simulations have the same
imposed flux $U=27$ mm/h.

To quantify the decrease of the infiltration depth due to freezing, we define
the depth of infiltration $L_{\mathrm{infil}}$ over a fixed time
$t_{\mathrm{f}}$ by measuring the vertical distance from the top of the domain
to the front of infiltration using the horizontally-averaged liquid water
content $\overline{\mathrm{LWC}}_{x}=(\overline{\phi S})_{x}$ as shown in Fig.
5(a). In Fig. 5(b), we plot $L_{\mathrm{infil}}$ at $t_{\mathrm{f}}=212$ min
as a function of both $T_{i,0}$ and $U$. We find that, for a fixed $U$, the
depth of infiltration decreases monotonically with decreasing $T_{i,0}$;
meanwhile, at a fixed $T_{i,0}$, $L_{\mathrm{infil}}$ monotonically increases
with increasing $U$. Applying the appropriate scaling to our data and defining
the effective infiltration velocity as
$\overline{u}=L_{\mathrm{infil}}/t_{\mathrm{f}}$, we find that the curves in
Fig. 5(b) can be collapsed into a master curve prescribing the normalized
effective infiltration velocity $\overline{u}/U$ as a function of Da, as shown
in Fig. 5(c). In particular, the master curve reveals two distinct regimes of
infiltration behavior identified by the threshold $\mathrm{Da}=40$. When
$\mathrm{Da}<40$, the rate of freezing is comparable to the flow rate, and the
nonlinear interplay between infiltration and freezing leads to complex
behavior in $\overline{u}/U$. In this regime, increasing $\mathrm{Da}$ reduces
$\overline{u}/U$ drastically, and the exact relationship between
$\overline{u}/U$ and $\mathrm{Da}$ is highly nonlinear and weakly dependent on
$U$. When $\mathrm{Da}>40$, the rate of freezing dominates the rate of
advection, and $\overline{u}/U$ becomes controlled solely by the rate of
freezing and can be prescribed by a linear function of $\mathrm{Da}$.

Figure 5:  Infiltration depths measured at $t_{\mathrm{f}}=212$ min over
different flow and thermal conditions: (a) At a given time, we calculate the
liquid water content $\mathrm{LWC}=\phi S$ (left) and compute its horizontal
average (right). The effective infiltration depth $L_{\mathrm{infil}}$ is then
computed as the depth where the integrated volume of
$\overline{\mathrm{LWC}}_{x}$ reaches below 2.5
$\frac{\mathrm{m}^{3}}{\mathrm{m}^{2}}$. (b) The effective infiltration depth
$L_{\mathrm{infil}}$ at $t_{\mathrm{f}}=212$ min is controlled by both
$T_{i,0}$ and the imposed flux $U$. (c) We define
$\overline{u}=L_{\mathrm{infil}}/t_{\mathrm{f}}$ and compute the normalized
effective infiltration rate $\overline{u}/U$ as a function of $\mathrm{Da}$.

###  V.2 Porosity structure of a freezing finger

Our simulation results thus far demonstrate that the initial fingering pattern
creates heterogeneous regions of low porosity in an initially homogeneous
porous medium, as demonstrated in Fig. 3(b) and Fig. 4. The reduction in
porosity is directly proportional to the rate of freezing and is focused
predominantly at the macroscopic ice-water interface where thermal gradients
that drive freezing are the strongest (Eq. (1)). In this section, we study the
detailed porosity structure of a freezing finger by considering the formation
and infiltration of a single finger into unsaturated permeable ice with
subfreezing temperature $T_{i,0}$. We conduct finer resolution simulations
than presented in Sec. V.1 and V.3, using a domain size of 0.1 m $\times$ 0.5
m ($L_{x}\times L_{y}$) with 256 $\times$ 1280 elements ($N_{x}\times N_{y}$).
For computational efficiency, we impose a thin “pinning” layer at the top of
the domain with a lower permeability to promote the fingering instability more
quickly. For the hydraulic properties of the porous medium, we use
approximated parameter values from Ref. [49] for coarse snow.

We perform a suite of these simulations by systematically varying $U$ and
$\mathrm{Ste}$ using the same values in Sec. V.1. To analyze the results, we
consider a cross-section of $\phi$ at $y=0.375$ m [Fig. 6(a)] at a
sufficiently late time in which local thermal equilibrium is achieved. The
porosity change inside the finger is small due to the fast downward advection
of thermal content which quickly brings the inner region of the finger into
thermal equilibrium. Meanwhile, the outer edges of the finger see sustained
and significant freezing, as its nonequilibrium state is maintained by the
diffusive heat exchange with the non-infiltrated region [Fig. 6(b)]. We note
however, that the ice-water interface eventually approaches thermal
equilibrium [Fig. 3(c)] as well, allowing the wall thickness to approach a
plateau value as the freezing process halts.

Figure 6:  Freezing infiltration of a single finger: (a) Porosity snapshot of
a single finger after infiltrating into porous ice initially at $T_{i,0}=-10$
∘C. (b) Porosity along a cutline at 0.125 m below the top of the domain at
multiple times. Inset shows a zoomed-in region of the frozen wall formed over
time until it reaches an equilibrium value. Freezing begins at the macroscopic
ice-water interface and propagates inwards where there is available water to
freeze. The Stefan number Ste alone controls the (c) minimum porosity for both
the frozen wall region (bottom curve) and the interior of the finger (top
curve), as well as (d) the thickness of the resulting frozen wall.

Because we only explore a range of relatively low $U$ that does not influence
the finger width, we observe that $U$ only affects the downward finger
velocity and does not modify the finger geometry or spanwise heat transfer. As
a result, the wall thickness and changes in porosity are only affected by Ste.
The equilibrium configuration of the porosity field exhibits an annulus of a
low-porosity region surrounding a still-porous region in which water is still
able to percolate through [Fig. 6(a)]. We find that both the wall thickness
[Fig. 6(d)] and the porosity of the ice walls [Fig. 6(c)] are only functions
of Ste and are independent of the imposed meltwater velocity $U$. The porosity
inside the finger decreases only slightly with Ste, as shown in Fig. 6(c).

###  V.3 Formation of secondary fingers

The reduction in the effective infiltration velocity $\overline{u}$ caused by
freezing can be explained, to the first order, by the consumption of melt and
the reduction in local permeability. Here, we highlight an additional
mechanism, which we term secondary fingering, that weakens the localized
nature of the infiltration flow paths and decreases $\overline{u}$.

During unstable infiltration, an initial set of fingers with roughly uniform
spacing forms, which we denote as the primary fingers as shown in Fig. 7(a) at
$t=80$ min. The majority of freezing occurs at the ice-water interface of the
primary fingers due to high thermal gradients across the interface (Fig. 6).
The latent heat release and accompanying diffusive heat transfer from this
freezing warms the unsaturated region between each finger, leading to a region
of uniform thermal equilibrium behind the fingering front [Fig. 3(c)]. As the
regions between each finger have remained unsaturated, local permeability is
higher than the wetted regions that have undergone freezing. With flow still
supplied at the top of the domain, new preferential pathways, which we term
the secondary fingers, emerge in these regions of higher permeability and
local thermodynamic equilibrium, as shown in Fig. 7(a) for $t>80$ min. Since
the conditions for flow are more favorable in these regions, the secondary
fingers quickly propagate towards the front near the primary fingers until
they reach the front and begin to undergo freezing, as shown in Fig. 7(b). The
formation and propagation of secondary fingers will slow down the overall
propagation of the wetting front, as water available for transport is now
diverted into the unsaturated regions in between the primary fingers.

Figure 7:  Dynamics of secondary finger formation: (a) Secondary fingers form
between primary fingers in regions of local thermal equilibrium. (b) Secondary
fingers propagate more quickly than primary fingers, and reduce the
infiltration velocity of primary fingers as they redirect water content
available for transport. In this simulation, secondary fingers are formed at
$t=93$ min (dashed black line). A moving average filter is applied to smooth
the average finger velocity curves.

The reduction in effective infiltration rate due to the formation of secondary
fingers is also corroborated by a reduction in the degree of channelization of
the flow structure. To demonstrate this, we quantify the lateral flow path
distribution by measuring the variance in saturation. We introduce
$\overline{S}_{y}(x)$ as the $y$-averaged saturation along the $x$-axis. We
then define the degree of flow channelization as the horizontal variance of
$\overline{S}_{y}$:

$\sigma^{2}=\langle{\overline{S}}^{2}_{y}\rangle-\langle\overline{S}_{y}\rangle^{2}.$
(26)

For a uniform flow, regardless of saturation overshoot in the $y$-direction,
$\overline{S}_{y}$ will be constant along $x$ and thus $\sigma^{2}$ will be
zero. In Fig. 8, we demonstrate that, as the formation of secondary fingers
redistributes flow into regions where the primary fingers do not visit, the
area of the unsaturated region decreases as more channels occupy the domain,
and thus $\sigma^{2}$ decreases. We also note that before the onset of
secondary fingers, $\sigma^{2}$ increases with increasing $\mathrm{Ste}$ and
can even surpass $\sigma^{2}$ of flow configurations without secondary
fingers. In these scenarios, more freezing occurs at the ice-water interface
of each finger for higher Ste (Sec. V.2), and thus the effective finger width
decreases, leading to sharper gradients in $\overline{S}_{y}(x)$ and therefore
a higher variance. Additionally, for $\mathrm{Ste}>0.1$ ($T_{i,0}<-7.5$ ∘C),
some fingers have the potential to meander and begin to propagate diagonally.
Only the secondary fingers exhibit this meandering behavior, as the porosity
structure at the site of their formation is heterogeneous due to freezing from
the primary set of fingers. This can cause slightly more freezing at one side
of the meandering finger due to differing local thermodynamic conditions on
either side. For an initially homogeneous medium, we postulate that this small
amount of phase change causes a “runaway” effect, causing the meandering
finger to continue propagating diagonally until it combines with another
finger. The coalescence of two fingers here will also raise $\sigma^{2}$, as
shown in Fig. 8(a) insets.

Figure 8:  Evolution of the variance $\sigma^{2}$ of the vertically-averaged
saturation $\overline{S}_{y}$ over time for Ste
$\in\\{0.031,0.063,0.094,0.126\\}$ (increasing in direction of arrows) and for
three different Pe. Snapshots in (a) show the saturation field at various
times and demonstrate the decrease in $\sigma^{2}$ due to the formation of
secondary fingers, and increases in $\sigma^{2}$ past $t=150$ min due to
secondary finger coalescence.

In Sec. V.1, we demonstrate that the competition between the rate of phase
change and the rate of advection, or the Damköhler number, controls the
normalized effective infiltration rate $\overline{u}/U$. Here, we find that
the compounding effects of phase change, as controlled by Ste, and thermal
advection, as controlled by Pe, predict the formation of secondary fingers. In
Fig. 9, we present a Ste–Pe phase diagram that predicts the emergence of
secondary fingers. Both sufficiently high rates of freezing around the primary
fingers (high Ste) as well as sufficient flow available for transport (high
Pe) are required to form secondary fingers. Thus, the product of these
nondimensional numbers, $\mathrm{Ste}\mathrm{Pe}$ controls the emergence of
secondary fingers. As shown in Fig. 9, we find that secondary fingers only
appear for $\mathrm{Ste}\mathrm{Pe}>10^{-4}$ (dashed curve).

Figure 9:  Phase diagram predicting the emergence of secondary fingers. The
compounding effects of sufficiently high thermal advection and rates of
freezing, as quantified by the nondimensional product StePe, controls the
emergence of secondary fingers.

##  VI Conclusions

In this paper, we study the effects of freezing on the unstable infiltration
of water into porous media. We use a continuum model [42] to describe gravity-
driven infiltration coupled with thermal transport and ice-water phase change
in an initially unsaturated and subfreezing porous media. We present high-
resolution numerical simulations of the model applied to various thermal
conditions and flow rates and identify key nondimensional parameters that
govern the patterns of flow and the resulting frozen structures. Our results
demonstrate that the effective infiltration rate can be well characterized by
the freezing Damköhler number (Sec. V.1). In particular, we identify a regime
for $\mathrm{Da}<40$ in which thermal advection and the freezing kinetics are
in strong competition, leading to a nonlinear decrease in the effective
infiltration rate with increasing $\mathrm{Da}$. We also identify a linear
regime ($\mathrm{Da}>40$) where the effects of freezing dominate the
infiltration process and the effective infiltration rate decreases linearly
with increasing Da. In Sec. V.2, we focus on the characteristics of vertical
structures, or ice pipes, and their formation. We find that the ice pipes are
initially created by freezing at the macroscopic ice-water interface while the
interior region of the finger remains porous and continues to facilitate
downward transport of water. Upon future cycles of environmental forcing, this
saturated region has the possibility to fully freeze which could result in a
solid column of ice, as observed in the field [20, 21]. In addition to the
consumption of liquid water and the decrease in local permeability, we
identify a third mechanism—the emergence of secondary fingers—that reduces the
effective infiltration rate by diverting water into new flow channels behind
the initial fingering front (Sec V.3). We find that both sufficient rates of
freezing and adequate thermal transport, as described by the nondimensional
product $\mathrm{Ste}\mathrm{Pe}$, is required to produce secondary fingers.
By measuring the variance of $\overline{S}_{y}(x)$, we show that the creation
and propagation of secondary fingers reduces the degree of channelization of
the flow field.

An important assumption we make in this study is that local freezing does not
alter the capillary pressure function $\psi(S)$. This function, also called
the water retention curve, has been documented only for a handful of various
snow types, but only under isothermal flow conditions [49, 50, 17]. In
freezing flows through porous media, both the volume-averaged ice grain
diameter $d_{i}$ and the ice fraction $(1-\phi)$ would increase at the site of
phase change, and the capillary suction force may strengthen or weaken locally
and modify the infiltration dynamics. Future studies involving numerical
simulations at the pore-scale and/or laboratory experiments would be valuable
in determining how the freezing process affects $\psi(S)$. Further, as natural
snowpack is often a stratified, heterogeneous porous medium, the interaction
between gravity-driven flow and discretely layered porous media with varying
hydraulic and capillary properties, and in particular the formation of ice
lenses, will be a topic of future studies.

## References

  * Bouma and Dekker [1978] J. Bouma and L. W. Dekker, A case study on infiltration into dry clay soil I. Morphological observations, Geoderma 20, 27 (1978).
  * Beven and Germann [1982] K. Beven and P. Germann, Macropores and water flow in soils, Water Resources Research 18, 1311 (1982).
  * Ritsema  _et al._ [1996] C. J. Ritsema, T. S. Steenhuis, J.-Y. Parlange, and L. W. Dekker, Predicted and observed finger diameters in field soils, Geoderma 70, 185 (1996).
  * Campbell  _et al._ [2006] F. M. A. Campbell, P. W. Nienow, and R. S. Purves, Role of the supraglacial snowpack in mediating meltwater delivery to the glacier system as inferred from dye tracer investigations, Hydrological Processes 20, 969 (2006).
  * Clerx  _et al._ [2022] N. Clerx, H. Machguth, A. Tedstone, N. Jullien, N. Wever, R. Weingartner, and O. Roessler, In situ measurements of meltwater flow through snow and firn in the accumulation zone of the SW Greenland Ice Sheet, The Cryosphere 16, 4379 (2022).
  * Glass and Nicholl [1996] R. J. Glass and M. J. Nicholl, Physics of gravity fingering of immiscible fluids within porous media: An overview of current understanding and selected complicating factors, Geoderma 70, 133 (1996).
  * Selker  _et al._ [1992] J. Selker, T. Steenhuis, and J.-Y. Parlange, Wetting front instability in homogeneous sandy soils under continuous infiltration, Soil Science Society of America Journal 56, 1346 (1992).
  * Yao and Hendrickx [2001] T. Yao and J. Hendrickx, Stability analysis of the unsaturated water flow equation: 2. Experimental verification, Water Resources Research 37, 1875 (2001).
  * Wei  _et al._ [2014] Y. Wei, C. M. Cejas, R. Barrois, R. Dreyfus, and D. J. Durian, Morphology of rain water channeling in systematically varied model sandy soils, Physical Review Applied 2, 1 (2014).
  * Marsh and Woo [1984] P. Marsh and M.-K. Woo, Wetting front advance and freezing of meltwater within a snow cover: 1. Observations in the Canadian Arctic, Water Resources Research 20, 1853 (1984).
  * Glass  _et al._ [1989a] R. J. Glass, G. H. Oosting, and T. S. Steenhuis, Preferential solute transport in layered homogeneous sands as a consequence of wetting front instability, Journal of Hydrology 110, 87 (1989a).
  * Kung  _et al._ [2000] K.-J. Kung, T. Steenhuis, E. Kladivko, T. Gish, G. Bubenzer, and C. Helling, Impact of preferential flow on the transport of adsorbing and non-adsorbing tracers, Soil Science Society of America Journal 64, 1290 (2000).
  * Würzer  _et al._ [2017] S. Würzer, N. Wever, R. Juras, M. Lehning, and T. Jonas, Modelling liquid water transport in snow under rain-on-snow conditions – considering preferential flow, Hydrology and Earth System Sciences 21, 1741 (2017).
  * Abermann  _et al._ [2019] J. Abermann, M. Eckerstorfer, E. Malnes, and B. U. Hansen, A large wet snow avalanche cycle in West Greenland quantified using remote sensing and in situ observations, Natural Hazards 97, 517 (2019).
  * Kinar and Pomeroy [2015] N. J. Kinar and J. W. Pomeroy, Measurement of the physical properties of the snowpack, Reviews of Geophysics 53, 481 (2015).
  * Waldner  _et al._ [2004] P. A. Waldner, M. Schneebeli, U. Schultze-Zimmermann, and H. Flühler, Effect of snow structure on water flow and solute transport, Hydrological Processes 18, 1271 (2004).
  * Katsushima  _et al._ [2013] T. Katsushima, S. Yamaguchi, T. Kumakura, and A. Sato, Experimental analysis of preferential flow in dry snowpack, Cold Regions Science and Technology 85, 206 (2013).
  * Woo  _et al._ [1982] M.-K. Woo, R. Heron, and P. Marsh, Basal ice in high arctic snowpacks, Arctic and Alpine Research 14, 251 (1982).
  * Williams  _et al._ [2000] M. W. Williams, M. Rikkers, and W. T. Pfeffer, Ice columns and frozen rills in a warm snowpack, Green Lakes Valley, Colorado, USA, Hydrology Research 31, 169 (2000).
  * Miller  _et al._ [2020] O. Miller, D. K. Solomon, C. Miège, L. Koenig, R. Forster, N. Schmerr, S. R. M. Ligtenberg, A. Legchenko, C. I. Voss, L. Montgomery, and J. R. McConnell, Hydrology of a perennial firn aquifer in Southeast Greenland: An overview driven by field data, Water Resources Research 56, 10.1029/2019WR026348 (2020).
  * Bouchard  _et al._ [2024] B. Bouchard, D. F. Nadeau, F. Domine, N. Wever, A. Michel, M. Lehning, and P.-E. Isabelle, Impact of rain-on-snow events on snowpack structure and runoff under a boreal canopy, EGUsphere 2024, 1 (2024).
  * Benson [1960] C. S. Benson,  _Stratigraphic studies in the snow and firn of the Greenland ice sheet_ , Ph.D. thesis, California Institute of Technology (1960).
  * Echelmeyer  _et al._ [1992] K. Echelmeyer, W. Harrison, T. Clarke, and C. Benson, Surficial glaciology of Jakobshavns Isbræ, west Greenland: Part II. Ablation, accumulation and temperature, Journal of Glaciology 38, 169 (1992).
  * Heilig  _et al._ [2018] A. Heilig, O. Eisen, M. MacFerrin, M. Tedesco, and X. Fettweis, Seasonal monitoring of melt and accumulation within the deep percolation zone of the greenland ice sheet and comparison with simulations of regional climate modeling, The Cryosphere 12, 1851 (2018).
  * Culberg  _et al._ [2021] R. Culberg, D. M. Schroeder, and W. Chu, Extreme melt season ice layers reduce firn permeability across Greenland, Nature Communications 12, 2336 (2021).
  * Van Pelt  _et al._ [2012] W. Van Pelt, J. Oerlemans, C. Reijmer, V. Pohjola, R. Pettersson, and J. Van Angelen, Simulating melt, runoff and refreezing on Nordenskiöldbreen, Svalbard, using a coupled snow and energy balance model, The Cryosphere 6, 641 (2012).
  * Van Pelt  _et al._ [2016] W. J. Van Pelt, V. A. Pohjola, and C. H. Reijmer, The changing impact of snow conditions and refreezing on the mass balance of an idealized Svalbard glacier, Frontiers in Earth Science 4, 102 (2016).
  * Avanzi  _et al._ [2016] F. Avanzi, H. Hirashima, S. Yamaguchi, T. Katsushima, and C. De Michele, Observations of capillary barriers and preferential flow in layered snow during cold laboratory experiments, Cryosphere 10, 2013 (2016).
  * Lombardo  _et al._ [2023] M. Lombardo, P. Lehmann, A. Kaestner, A. Fees, A. Van Herwijnen, and J. Schweizer, A method for imaging water transport in soil–snow systems with neutron radiography, Annals of Glaciology , 1 (2023).
  * Adachi  _et al._ [2020] S. Adachi, S. Yamaguchi, T. Ozeki, and K. Kose, Application of a magnetic resonance imaging method for nondestructive, three-dimensional, high-resolution measurement of the water content of wet snow samples, Frontiers in Earth Science 8, 179 (2020).
  * Yamaguchi  _et al._ [2023] S. Yamaguchi, S. Adachi, and S. Sunako, A novel method to visualize liquid distribution in snow: superimposition of MRI and X-ray CT images, Annals of Glaciology , 1 (2023).
  * De Michele  _et al._ [2013] C. De Michele, F. Avanzi, A. Ghezzi, and C. Jommi, Investigating the dynamics of bulk snow density in dry and wet conditions using a one-dimensional model, The Cryosphere 7, 433 (2013).
  * Wever  _et al._ [2016] N. Wever, S. Würzer, C. Fierz, and M. Lehning, Simulating ice layer formation under the presence of preferential flow in layered snowpacks, The Cryosphere 10, 2731 (2016).
  * Meyer and Hewitt [2017] C. R. Meyer and I. J. Hewitt, A continuum model for meltwater flow through compacting snow, The Cryosphere 11, 2799 (2017).
  * Webb  _et al._ [2021] R. W. Webb, K. Jennings, S. Finsterle, and S. R. Fassnacht, Two-dimensional liquid water flow through snow at the plot scale in continental snowpacks: simulations and field data comparisons, The Cryosphere 15, 1423 (2021).
  * Hirashima  _et al._ [2014] H. Hirashima, S. Yamaguchi, and T. Katsushima, A multi-dimensional water transport model to reproduce preferential flow in the snowpack, Cold Regions Science and Technology 108, 80 (2014).
  * Hirashima  _et al._ [2019] H. Hirashima, F. Avanzi, and N. Wever, Wet-snow metamorphism drives the transition from preferential to matrix flow in snow, Geophysical Research Letters 46, 14548 (2019).
  * Leroux and Pomeroy [2019] N. R. Leroux and J. W. Pomeroy, Simulation of capillary pressure overshoot in snow combining trapping of the wetting phase with a nonequilibrium Richards equation model, Water Resources Research 55, 236 (2019).
  * Leroux  _et al._ [2020] N. R. Leroux, C. B. Marsh, and J. W. Pomeroy, Simulation of preferential flow in snow with a 2-D non-equilibrium Richards model and evaluation against laboratory data, Water Resources Research 56, e2020WR027466 (2020).
  * Illangasekare  _et al._ [1990] T. H. Illangasekare, R. J. Walter Jr, M. F. Meier, and W. T. Pfeffer, Modeling of meltwater infiltration in subfreezing snow, Water Resources Research 26, 1001 (1990).
  * Leroux and Pomeroy [2017] N. R. Leroux and J. W. Pomeroy, Modelling capillary hysteresis effects on preferential flow through melting and cold layered snowpacks, Advances in Water Resources 107, 250 (2017).
  * Moure  _et al._ [2023] A. Moure, N. Jones, J. Pawlak, C. Meyer, and X. Fu, A thermodynamic nonequilibrium model for preferential infiltration and refreezing of melt in snow, Water Resources Research , e2022WR034035 (2023).
  * Libbrecht [2017] K. G. Libbrecht, Physical dynamics of ice crystal growth, Annual Review of Materials Research 47, 271 (2017).
  * Cueto-Felgueroso and Juanes [2008] L. Cueto-Felgueroso and R. Juanes, Nonlocal interface dynamics and pattern formation in gravity-driven unsaturated flow through porous media, Physical Review Letters 101, 244504 (2008).
  * Cueto-Felgueroso and Juanes [2009a] L. Cueto-Felgueroso and R. Juanes, A phase field model of unsaturated flow, Water resources research 45 (2009a).
  * Beljadid  _et al._ [2020] A. Beljadid, L. Cueto-Felgueroso, and R. Juanes, A continuum model of unstable infiltration in porous media endowed with an entropy function, Advances in Water Resources 144, 10.1016/j.advwatres.2020.103684 (2020).
  * Calonne  _et al._ [2012] N. Calonne, C. Geindreau, F. Flin, S. Morin, B. Lesaffre, S. Rolland du Roscoat, and P. Charrier, 3-d image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy, The Cryosphere 6, 939 (2012).
  * Brooks and Corey [1966] R. H. Brooks and A. T. Corey, Properties of porous media affecting fluid flow, Journal of the irrigation and drainage division 92, 61 (1966).
  * Yamaguchi  _et al._ [2010] S. Yamaguchi, T. Katsushima, A. Sato, and T. Kumakura, Water retention cuasdfasdfrve of snow with different grain sizes, Cold Regions Science and Technology 64, 87 (2010).
  * Yamaguchi  _et al._ [2012] S. Yamaguchi, K. Watanabe, T. Katsushima, A. Sato, and T. Kumakura, Dependence of the water retention curve of snow on snow characteristics, Annals of Glaciology 53, 6 (2012).
  * Wang  _et al._ [1998] Z. Wang, J. Feyen, and D. E. Elrick, Prediction of fingering in porous media, Water Resources Research 34, 2183 (1998).
  * Cueto-Felgueroso and Juanes [2009b] L. Cueto-Felgueroso and R. Juanes, Stability analysis of a phase-field model of gravity-driven unsaturated flow through porous media, Physical Review E 79, 036301 (2009b).
  * Weitz  _et al._ [1987] D. Weitz, J. Stokes, R. Ball, and A. Kushnick, Dynamic capillary pressure in porous media: Origin of the viscous-fingering length scale, Physical review letters 59, 2967 (1987).
  * Glass  _et al._ [1989b] R. Glass, T. Steenhuis, and J.-Y. Parlange, Wetting front instability: 2. Experimental determination of relationships between system parameters and two-dimensional unstable flow field behavior in initially dry porous media, Water Resources Research 25, 1195 (1989b).
  * Hürkamp  _et al._ [2017] K. Hürkamp, S. Tafelmeier, and J. Tschiersch, Influence of melt-freeze-cycles on the radionuclide transport in homogeneous laboratory snowpack, Hydrological Processes 31, 1360 (2017).
  * Domine  _et al._ [2007] F. Domine, A. Taillandier, and W. R. Simpson, A parameterization of the specific surface area of seasonal snow for field use and for models of snowpack evolution, Journal of Geophysical Research: Earth Surface 112, 2006JF000512 (2007).
  * Hughes  _et al._ [2005] T. J. Hughes, J. A. Cottrell, and Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer methods in applied mechanics and engineering 194, 4135 (2005).
  * Balay  _et al._ [2023] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang, PETSc Web page, https://petsc.org/ (2023).
  * Dalcin  _et al._ [2016] L. Dalcin, N. Collier, P. Vignal, A. Côrtes, and V. Calo, PetIGA: A framework for high-performance isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 308, 151 (2016).