of 16
PHYSICAL REVIEW FLUIDS
9
, 123802 (2024)
Pattern formation of freezing infiltration in porous media
Nathan D. Jones
,
*
Adrian Moure
,
and Xiaojing Fu
Department of Mechanical and Civil Engineering,
California Institute of Technology
,
1200 E California Blvd, Pasadena, California 91125, USA
(Received 16 May 2024; accepted 13 November 2024; published 2 December 2024)
Gravity-driven infiltration of liquid water into unsaturated porous media can be a spa-
tially 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 be-
lieved 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 nonequilib-
rium 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 an 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. Finally, we identify a regime under high Damköhler number in
which freezing-induced permeability reduction completely clogs flow.
DOI:
10.1103/PhysRevFluids.9.123802
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
5
] and studied in laboratory experiments [
6
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.
*
Contact author: ndjones@caltech.edu
Contact author: amoure@caltech.edu
Contact author: rubyfu@caltech.edu
2469-990X/2024/9(12)/123802(16)
123802-1
©2024 American Physical Society
JONES, MOURE, AND FU
FIG. 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. [
22
] of excavated refrozen structures created after preferential infiltration in
snow. (c) Problem setup for simulations presented in this study. Liquid water at
T
melt
=
0
C infiltrates through
the top of the domain and interacts with the subfreezing porous ice.
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 perme-
ability structure of the porous ice, and could be a precursor to hydrofracturing [
15
]. 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 [
16
18
]. 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 [
19
23
]. 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 [
24
,
25
], both from excavated snow pits [Fig.
1(b)
]aswellas
radar detection [
26
,
27
] 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 [
28
,
29
], and the lasting structures affect the bulk hydraulic, thermodynamic, and structural
properties of the snowpack.
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 [
5
,
18
,
30
] 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 [
30
]or,
more recently, the employment of nonintrusive techniques such as neutron radiography [
31
]or
magnetic resonance imaging [
32
,
33
]. Many numerical studies have been performed that investigate
the formation of refrozen structures (in particular, ice lenses) using one-dimensional models [
34
38
]
without considering the flow instability. Preferential flow in porous ice has also been investigated
through numerical simulations, but these studies mainly consider isothermal conditions [
39
42
]
or are limited by the assumption of local thermodynamic equilibrium [
43
,
44
]. 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.
123802-2
PATTERN FORMATION OF FREEZING INFILTRATION IN ...
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 [
45
] 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 nonequilibrium (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 two-dimensional (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 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. Finally, we identify a regime of very high Damköhler number in which
freezing dominates over thermal advection, and a low-permeability layer develops that completely
clogs flow.
II. GOVERNING EQUATIONS
We consider a porous matrix of ice, water, and air, in which the volume-averaged quantities water
saturation
S
(
x
,
t
) and porosity
φ
(
x
,
t
) determine the phase composition throughout the domain. The
volumetric liquid water content (LWC) is defined as LWC
=
φ
S
and the solid fraction of ice is
1
φ
. Flow of water is governed by a Darcy velocity
u
(
x
,
t
). The temperature of the water and ice
phases are
T
w
(
x
,
t
) and
T
i
(
x
,
t
), respectively. The full details of this model can be found in Ref. [
45
].
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:
(1
φ
)
t
=−
R
m
W
SSA
(
T
int
T
melt
)
,
(1)
(
φ
S
)
t
+
·
u
=
ρ
i
ρ
w
R
m
W
SSA
(
T
int
T
melt
)
,
(2)
where
ρ
i
and
ρ
w
are the densities of ice and water, respectively,
T
int
is the temperature of the ice-
water interface, and
T
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
sol
β
sol
), where
c
p
,
w
is the specific heat capacity of liquid water,
L
sol
is the
latent heat of solidification of ice, and
β
sol
is the kinetic attachment coefficient for ice growth from
liquid water [
46
]. The wet specific surface area
W
SSA
represents the ice-water interfacial area per
unit volume:
W
SSA
=
S
SSA
0
φ
0
ln(
φ
0
)
φ
ln(
φ
)
,
(3)
where SSA
0
is the initial specific surface area. The interfacial temperature
T
int
is an averaged pore-
scale quantity that is defined from the volume integrated Stefan condition:
T
int
T
melt
L
sol
/
c
p
,
w
=−
β
sol
v
n
,
(4)
123802-3
JONES, MOURE, AND FU
where
v
n
is the volume-averaged velocity of the fluid-solid interface and is taken to be positive for
freezing. Following Moure
et al.
[
45
], the interface temperature takes the form:
T
int
(
T
i
,
T
w
)
=
k
i
r
i
T
i
+
k
w
r
w
T
w
c
p
,
w
ρ
w
β
sol
+
k
i
r
i
+
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.
[
45
]). Note that the right-hand side expressions in
Eqs. (
1
) and (
2
) differ by a factor of
ρ
i
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:
[(1
φ
)
T
i
]
t
=
·
(
D
i
(1
φ
)
T
i
)
+
W
SSA
D
i
T
int
T
i
r
i
,
(6)
(
φ
ST
w
)
t
+
·
(
u
T
w
)
=
·
(
D
w
φ
S
T
w
)
+
W
SSA
D
w
T
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 Beljadid, Cueto-Felgueroso, and Juanes [
47
49
]:
u
=−
K
s
(
φ
)
k
r
(
S
)

(
S
)
,
(8)
where
K
s
(
φ
) is the saturated hydraulic conductivity,
k
r
(
S
) is the relative permeability, and

(
S
)is
the total flow potential. We use an empirical expression for the snow hydraulic conductivity [
50
]:
K
s
(
φ
)
=
3
(
d
i
2
)
2
ρ
w
g
μ
w
exp
[
0
.
013
ρ
i
(1
φ
)
]
,
(9)
where
d
i
is the average ice grain diameter,
g
is the gravitational acceleration, and
μ
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 [
51
]:
k
r
(
S
)
=
(
S
S
r
1
S
r
)
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
×
10
3
and
a
=
5. The expression for the flow potential is [
47
]:

(
S
)
=
y
ψ
(
S
)
κ
·
(
κ
S
)
,
(11)
where
y
is the vertical coordinate (increasing with height),
ψ
(
S
) is capillary suction derived from
Leverett scaling, and
κ
(
S
) is the capillary pinning function:
ψ
(
S
)
=
h
cap
S
1
α
{
1
exp
[
β
(
S
ν
e
)
]
(
1
+
β
α
α
1
S
)}
,
(12)
κ
(
S
)
=
h
2
cap
S
0
ψ
(
S
)d
S
=
h
3
cap
α
α
1
S
α
1
α
{
1
exp
[
β
(
S
ν
e
)
]
}
,
(13)
where
h
cap
,
α
,
β
, and
ν
e
are constants that are calibrated from water retention curves for a given
porous material [
18
,
52
,
53
]. These functions along with the values for these parameters are shown in
Fig.
2
. Note that in the case of
κ
=
0, the classical Richards equation for unsaturated flow in porous
media is recovered.
123802-4
PATTERN FORMATION OF FREEZING INFILTRATION IN ...
FIG. 2. Capillary pressure
ψ
and the capillary pinning function
κ
used in this work. Here, we plot
Eqs. (
12
)and(
13
) with
h
cap
=
2
.
5cm,
α
=
4,
β
=
22, and
ν
e
=
1. Parameters are calibrated from experimental
measurements from Ref. [
18
].
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
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:
|
T
0
|=|
T
i
,
0
T
melt
|
is the characteristic
temperature,
U
is the characteristic velocity,
L
=
SSA
1
0
is the length scale, and hence the product
t
c
=
(SSA
0
U
)
1
is the characteristic time scale. We use SSA
1
0
as the length scale as it provides a
more appropriate volume-averaged measurement of the average grain size. Note that SSA
1
0
d
i
FIG. 3. Snapshots of problem unknowns at
t
=
200 min after continuous infiltration for: (a) saturation
S
, (b) porosity
φ
, (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
=
1m, 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.
123802-5
JONES, MOURE, AND FU
for spherical ice grains. The nondimensional equations for mass and thermal conservation of each
phase then take the following form:
∂φ
t
=
Da
W
SSA
(
T
int
T
melt
)
,
(14)
(
φ
S
)
t
+
·
u
=
ρ
i
ρ
w
Da
W
SSA
(
T
int
T
melt
)
,
(15)
[(1
φ
)
T
i
]
t
=
Pe
1
i
[
·
(
(1
φ
)
T
i
)
+
W
SSA
(
T
int
T
i
)
]
,
(16)
(
φ
ST
w
)
t
+
·
(
u
T
w
)
=
Pe
1
w
[
·
(
φ
S
T
w
)
+
W
SSA
(
T
int
T
w
)
]
.
(17)
The Damköhler number Da sets the ratio between the rate of phase change and the rate of
advection:
Da
=
c
p
,
w
|
T
0
|
L
sol
U
β
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,
v
n
, and the rate of advection:
Da
=
Ste
U
β
sol
=
v
n
U
.
(19)
Here, Ste is the Stefan number, which sets the ratio between the amount of sensible heat to the latent
heat of solidification:
Ste
=
c
p
,
w
|
T
0
|
L
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:
Pe
i
=
UL
/
D
i
,
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 Pe
i
and Pe
w
differ only by the thermal diffusivity of each phase, we refer to Pe
=
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 (SSA
1
0
), which correlates with the grain size. However, the original isothermal
problem of gravity-driven infiltration bears its own scaling analysis that yields additional nondi-
mensional 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:
u
=−
K

s
(
φ
)
R
s
k
r


,
(22)


=
y

N
1
gr
ψ

(
S
)
N
nl
(
κ

·
(
κ

S
))
,
(23)
where
ψ

(
S
)
=
ψ
(
S
)
/
h
cap
,
κ

=
κ/
h
3
cap
, and
N
gr
and
N
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
,
48
,
54
]. The influence of these two additional nondimensional groups on the
stability of the flow pattern is investigated in Ref. [
55
], and they are defined:
N
gr
=
L
h
cap
,
N
nl
=
ρ
gL
3
,
(24)
123802-6
PATTERN FORMATION OF FREEZING INFILTRATION IN ...
where
is a coefficient in the nonlocal free energy potential that is associated with an effective
surface tension of the wetting front interface [
49
,
56
]. Following Ref. [
47
], we take
=
ρ
g
κ
, and
recognizing that
κ
h
3
cap
from Eq. (
13
), we arrive at the relation:
N
nl
=
N
3
gr
.
(25)
Increasing
N
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 [
48
,
55
,
57
].
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
gr
or
N
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 low flux ratios between 7
.
9
×
10
6
to 5
.
3
×
10
4
(
U
between 0.4 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 [
58
].
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
×
2m (
L
x
×
L
y
) as shown in Figs.
1(c)
and
3
. Initially, the domain is composed of dry porous ice with a uniform initial ice temperature
T
i
,
0

0
C, uniform porosity
φ
(
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
(
x
,
0)
=
S
r
=
1
×
10
3
. The initial specific surface area is estimated from
Ref. [
59
]tobeSSA
0
=
3514 m
1
. We approximate the kinetic attachment coefficient for ice growth
from Ref. [
46
]as
β
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 water flux
U
and the initial ice phase temperature
T
i
,
0
.
We va r y
U
between 0.4–27 mm/h to explore flow rates near values reported in Ref. [
58
] for natural
systems and consider initial ice temperatures between 0 to
20
C. The domain has 400
×
800
(
N
x
×
N
y
) mesh elements, with the exception of simulations presented in Sec.
VB
, in which we
use a finer mesh resolution of 0.4 mm/element. Our mesh convergence tests show that the solution
converges for mesh resolutions finer than 2.5 mm/element, which corresponds to the coarsest mesh
used in this study. Equations (
1
), (
2
), (
6
), (
7
), and (
8
) are solved using isogeometric analysis [
60
]—a
finite element method that implements B splines as basis functions. The time advancement scheme
is a generalized–
α
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 [
61
] and PetIGA [
62
].
V. RESULTS
A. 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
φ
atafixedtime(
t
f
=
212 min)
123802-7
JONES, MOURE, AND FU
FIG. 4. Snapshots of liquid water content (top) and porosity (bottom) at
t
f
=
212 min for varying
T
i
,
0
.All
simulations have the same imposed flux
U
=
27 mm/h.
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.
To quantify the decrease of the infiltration depth due to freezing, we define the depth of
infiltration
L
infil
over a fixed time
t
f
by measuring the vertical distance from the top of the domain to
the front of infiltration using the horizontally averaged liquid water content
LWC
x
=
(
φ
S
)
x
as shown
in Fig.
5(a)
.InFig.
5(b)
,weplot
L
infil
at
t
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
infil
monotonically increases with increasing
U
. Applying the appropriate scaling to our
data and defining the effective infiltration velocity as
u
=
L
infil
/
t
f
, we find that the curves in Fig.
5(b)
can be collapsed into a master curve prescribing the normalized effective infiltration velocity
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 Da
=
40. When 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
u
/
U
. In this regime, increasing Da reduces
u
/
U
drastically, and the exact
relationship between
u
/
U
and Da is highly nonlinear and weakly dependent on
U
. When Da
>
40,
FIG. 5. Infiltration depths measured at
t
f
=
212 min over different flow and thermal conditions: (a) At a
given time, we calculate the liquid water content LWC
=
φ
S
(left) and compute its horizontal average (right).
The effective infiltration depth
L
infil
is then computed as the depth where the integrated volume of
LWC
x
reaches below 2.5
m
3
m
2
. (b) The effective infiltration depth
L
infil
at
t
f
=
212 min is controlled by both
T
i
,
0
and
the imposed flux
U
.(c)Wedefine
u
=
L
infil
/
t
f
and compute the normalized effective infiltration rate
u
/
U
as a
function of Da.
123802-8
PATTERN FORMATION OF FREEZING INFILTRATION IN ...
FIG. 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.
the rate of freezing dominates the rate of advection, and
u
/
U
becomes controlled solely by the rate
of freezing and can be prescribed by a linear function of Da.
B. 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 Figs.
3(b)
and
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 Secs.
VA
and
VC
, using a domain size of 0
.
1m
×
0
.
5m (
L
x
×
L
y
) with 256
×
1280 elements (
N
x
×
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. [
52
] for coarse snow. We
perform a suite of these simulations by systematically varying
U
and Ste using the same values in
Sec.
VA
. To analyze the results, we consider a cross section of
φ
at
y
=
0
.
375 m [Fig.
6(a)
]ata
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 noninfiltrated region [Fig.
6(b)
]. We note, however, that the ice-water interface eventually
approaches thermal equilibrium as well, and the wall thickness reaches a plateau as the freezing
process halts.
123802-9
JONES, MOURE, AND FU
FIG. 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.
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 inflow velocity
U
. The porosity inside the
finger decreases only slightly with Ste, as shown in Fig.
6(c)
.
C. Formation of secondary fingers
The reduction in the effective infiltration velocity
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
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
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.
The reduction in effective infiltration rate due to the formation of secondary fingers is also
complemented 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
S
y
(
x
)asthe
y
-averaged saturation along the
x
axis. We then define the degree of flow
channelization as the horizontal variance of
S
y
:
σ
2
=
S
2
y
−
S
y

2
.
(26)
123802-10
PATTERN FORMATION OF FREEZING INFILTRATION IN ...
FIG. 8. Evolution of the variance
σ
2
of the vertically averaged saturation
S
y
over time for Ste
{
0
.
031
,
0
.
063
,
0
.
094
,
0
.
126
}
(increasing in the direction of arrows) and for three different Pe. Snapshots in
(a) show the saturation field at various times and demonstrate the decrease in
σ
2
due to the formation of
secondary fingers, and increases in
σ
2
past
t
=
150 min due to secondary finger coalescence.
For a uniform flow, regardless of saturation overshoot in the
y
direction,
S
y
will be constant along
x
and thus
σ
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
σ
2
decreases. We also note that
before the onset of secondary fingers,
σ
2
increases with increasing Ste and can even surpass
σ
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.
VB
), and thus the effective finger width decreases,
leading to sharper gradients in
S
y
(
x
) and therefore a higher variance. Additionally, for 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
σ
2
, as shown in Fig.
8(a)
insets.
D. The complete phase diagram of freezing infiltration
We now present a summary of the freezing infiltration patterns for varying flow rates and thermal
conditions in a Ste–Pe phase diagram in Fig.
9
. Results presented in Secs.
VA
VC
have only
explored the top half of the phase diagram (Pe

5
.
8
×
10
4
or
U

9 mm/h). In particular, the
majority of the phase diagram (red triangular markers) corresponds to simulations of
fingered flow
,
where only the primary fingering instability occurs. The northeast region of the phase diagram
(blue round markers) corresponds to conditions in which
secondary fingers
are observed. Lastly,
the southeast corner of the phase diagram (gray square markers) corresponds to conditions for
clogged flow
, where local permeability undergoes a fast and sharp decrease due to rapid freezing
that completely arrests the flow within the first few centimeters of the domain.
To explain the transition from
fingered flow
to the
secondary fingers
regime, we find that the
compounding effects of phase change, as controlled by Ste, and thermal advection, as controlled by
Pe, collectively predict the formation of secondary fingers. In order for secondary fingers to develop,
123802-11
JONES, MOURE, AND FU
FIG. 9. Ste–Pe phase diagram of flow behavior in log-log space. Alternate axes (
|
T
0
|
and
U
)show
the physical parameters that correspond to the nondimensional axes Ste and Pe, respectively, using SSA
0
=
3514 m
1
as the length scale. Snapshots (right) show saturation and porosity color maps for a typical flow
pattern in each regime: fingered flow (red triangles), secondary fingers (blue dots), and clogged flow (grey
squares). The dotted line corresponds to StePe
=
10
4
and the dashed line corresponds to Da
=
700.
high rates of freezing (high Ste) are necessary to sufficiently lower the local permeability within the
primary fingers as they propagate downwards. Additionally, there must be enough flow available
for transport (high Pe) to route flow towards the regions of thermal equilibrium near the top of the
domain between the primary fingers instead of supplying flow to the primary fingers themselves. If
Ste is not large enough, the primary fingers will continue to propagate downwards as the change in
permeability is not drastic enough to drive the formation of new, easier pathways of flow. Thus, the
product of these two nondimensional numbers, StePe, predicts the emergence of secondary fingers.
Specifically, the transition from the red triangular to blue dotted region is marked by StePe
10
4
(dotted line in Fig.
9
) for the parameter space explored in this study.
Conversely, to explain the transition from
fingered flow
to
clogged flow
in the region of high
freezing rates and low flow rates, we find that the nondimensional ratio Ste/Pe, which is proportional
to the Damköhler number, predicts the clogging behavior. In this context, we can interpret the
Damköhler number as the ratio between the macroscopic freezing front velocity and the imposed
flow velocity at the top boundary [Eq. (
19
)]. A sufficiently high Damköhler number suggests that
the freezing front velocity surpasses the supplied inflow rate, resulting in flow arrest. We find that
this threshold of complete arrest is near Da
700 (dashed line in Fig.
9
).
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 [
45
] 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
123802-12
PATTERN FORMATION OF FREEZING INFILTRATION IN ...
characterized by the freezing Damköhler number (Sec.
VA
). In particular, we identify a regime for
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 Da. We also identify a linear
regime (Da
>
40) where the effects of freezing dominate the infiltration process and the effective
infiltration rate decreases linearly with increasing Da. In Sec.
VB
, 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 [
22
,
23
]. 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.
VC
). By measuring the variance of
S
y
(
x
), we show that the creation
and propagation of secondary fingers reduces the degree of channelization of the flow field. We find
that the patterns of freezing infiltration can be well categorized in a Ste–Pe phase diagram marked
by three distinct regimes:
fingered flow
,
secondary fingers
, and
clogged flow
. To produce secondary
fingers, both high rates of freezing and adequate thermal transport, as described by a sufficiently
large nondimensional product StePe, is necessary. To enter the
clogged flow
regime, where flow
is completely arrested by the formation of a low-permeability region at the top of the domain,
as similarly described in Ref. [
63
], a sufficiently high Damköhler number is required. However,
to rigorously derive the transition conditions between the different regimes, theoretical techniques
such as linear stability analysis may be required. The stability of gravity-driven preferential flow
with freezing has not yet been investigated, and is the topic of future work.
An important assumption we make in this study is that local freezing does not alter the capillary
pressure function
ψ
(
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 [
18
,
52
,
53
]. In
freezing flows through porous media, both the volume-averaged ice grain diameter
d
i
and the ice
fraction (1
φ
) 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 numer-
ical simulations at the pore-scale and/or laboratory experiments would be valuable in determining
how the freezing process affects
ψ
(
S
).
We note that all results presented in this work are based on 2D simulations. The isothermal
gravity fingering model has been investigated in three dimensions in Ref. [
64
], which shows similar
fingering dynamics between 2D and 3D simulations. However, we do expect some differences
between 2D and 3D simulations as a result of a more realistic description of the macroscopic
air-water interface curvature in 3D. As perturbations in this macroscopic interface lead to the initial
fingering instability, differences may be seen in the time scales of the onset of the initial stability
and the finger spacing. These differences in the finger spacing and 3D morphology may also affect
the magnitude of the downward flux, and the threshold values for nondimensional groups that
predict flow behavior, such as describing normalized infiltration rate (Fig.
5
) and the emergence
of secondary fingers (Fig.
9
) may differ from those predicted by 2D simulations. Further, as natural
snowpack is a stratified, heterogeneous porous medium, the interaction between gravity-driven flow
and discretely layered porous media with varying hydraulic and capillary properties, will have a
profound effect on the effective infiltration rate compared to an initially homogeneous system. In
this case, both laboratory experiments and 3D simulations may also prove useful in describing the
ponding effect observed at the interface between these layers and the formation of macroscopic ice
lenses and slabs as found in nature, and this will be a topic of future studies.
123802-13
JONES, MOURE, AND FU
ACKNOWLEDGMENT
The authors acknowledge the partial support from the Resnick Sustainability Institute at Cal-
ifornia Institute of Technology and from the U.S. National Science Foundation under Grant No.
EAR-2243631.
[1] J. Bouma and L. W. Dekker, A case study on infiltration into dry clay soil I. Morphological observations,
Geoderma
20
, 27 (1978)
.
[2] K. Beven and P. Germann, Macropores and water flow in soils,
Water Resour. Res.
18
, 1311 (1982)
.
[3] 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)
.
[4] 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)
.
[5] 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)
.
[6] 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)
.
[7] J. Selker, T. Steenhuis, and J.-Y. Parlange, Wetting front instability in homogeneous sandy soils under
continuous infiltration,
Soil Sci. Soc. Am. J.
56
, 1346 (1992)
.
[8] T. Yao and J. Hendrickx, Stability analysis of the unsaturated water flow equation: 2. Experimental
verification,
Water Resour. Res.
37
, 1875 (2001)
.
[9] Y. Wei, C. M. Cejas, R. Barrois, R. Dreyfus, and D. J. Durian, Morphology of rain water channeling in
systematically varied model sandy soils,
Phys. Rev. Appl.
2
, 044004 (2014)
.
[10] P. Marsh and M.-K. Woo, Wetting front advance and freezing of meltwater within a snow cover: 1.
Observations in the Canadian Arctic,
Water Resour. Res.
20
, 1853 (1984)
.
[11] R. J. Glass, G. H. Oosting, and T. S. Steenhuis, Preferential solute transport in layered homogeneous sands
as a consequence of wetting front instability,
J. Hydrology
110
, 87 (1989)
.
[12] 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 Sci. Soc. Am. J.
64
, 1290 (2000)
.
[13] 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,
Hydrol. Earth Syst. Sci.
21
, 1741 (2017)
.
[14] 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,
Nat. Hazards
97
, 517 (2019)
.
[15] Y. Meng, R. Culberg, and C.-Y. Lai, Vulnerability of firn to hydrofracture: poromechanics modeling,
J. Glaciol. (2024), doi:
10.1017/jog.2024.47
.
[16] N. J. Kinar and J. W. Pomeroy, Measurement of the physical properties of the snowpack,
Rev. Geophys.
53
, 481 (2015)
.
[17] 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)
.
[18] T. Katsushima, S. Yamaguchi, T. Kumakura, and A. Sato, Experimental analysis of preferential flow in
dry snowpack,
Cold Regions Sci. Technol.
85
, 206 (2013)
.
[19] M.-K. Woo, R. Heron, and P. Marsh, Basal ice in high arctic snowpacks,
Arctic Alpine Res.
14
, 251
(1982)
.
[20] M. W. Williams, M. Rikkers, and W. T. Pfeffer, Ice columns and frozen rills in a warm snowpack, Green
Lakes Valley, Colorado, USA,
Hydrol. Research
31
, 169 (2000)
.
[21] N. F. Humphrey, J. T. Harper, and W. T. Pfeffer, Thermal tracking of meltwater retention in Greenland’s
accumulation area,
J. Geophys. Res. Earth Surface
117
, F01010 (2012)
.
123802-14
PATTERN FORMATION OF FREEZING INFILTRATION IN ...
[22] 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 Resour. Res.
56
, e2019WR026348
(2020)
.
[23] B. Bouchard, D. F. Nadeau, F. Domine, N. Wever, A. Michel, M. Lehning, and P.-E. Isabelle, Impact of
intercepted and sub-canopy snow microstructure on snowpack response to rain-on-snow events under a
boreal canopy,
The Cryosphere
18
, 2783 (2024)
.
[24] C. S. Benson, Stratigraphic studies in the snow and firn of the Greenland ice sheet, Ph.D. thesis, California
Institute of Technology, 1960.
[25] K. Echelmeyer, W. Harrison, T. Clarke, and C. Benson, Surficial glaciology of Jakobshavns Isbræ, West
Greenland: Part II. Ablation, accumulation and temperature,
J. Glaciol.
38
, 169 (1992)
.
[26] 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)
.
[27] R. Culberg, D. M. Schroeder, and W. Chu, Extreme melt season ice layers reduce firn permeability across
Greenland,
Nat. Commun.
12
, 2336 (2021)
.
[28] W. J. J. van Pelt, J. Oerlemans, C. H. Reijmer, V. A. Pohjola, R. Pettersson, and J. H. van Angelen,
Simulating melt, runoff and refreezing on Nordenskiöldbreen, Svalbard, using a coupled snow and energy
balance model,
The Cryosphere
6
, 641 (2012)
.
[29] 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,
Front. Earth Sci.
4
, 102 (2016)
.
[30] 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,
The Cryosphere
10
,
2013 (2016)
.
[31] 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,
Ann. Glaciol., 1 (2023)
.
[32] 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,
Front. Earth Sci.
8
, 179 (2020)
.
[33] S. Yamaguchi, S. Adachi, and S. Sunako, A novel method to visualize liquid distribution in snow:
Superimposition of MRI and X-ray CT images,
Ann. Glaciol., 1 (2023)
.
[34] 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)
.
[35] 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)
.
[36] C. R. Meyer and I. J. Hewitt, A continuum model for meltwater flow through compacting snow,
The
Cryosphere
11
, 2799 (2017)
.
[37] 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)
.
[38] M. A. Shadab, S. Adhikari, A. Rutishauser, C. Grima, and M. A. Hesse, A mechanism for ice layer
formation in glacial firn,
Geophys. Res. Lett.
51
, e2024GL109893 (2024)
.
[39] H. Hirashima, S. Yamaguchi, and T. Katsushima, A multi-dimensional water transport model to reproduce
preferential flow in the snowpack,
Cold Regions Sci. Technol.
108
, 80 (2014)
.
[40] H. Hirashima, F. Avanzi, and N. Wever, Wet-snow metamorphism drives the transition from preferential
to matrix flow in snow,
Geophys. Res. Lett.
46
, 14548 (2019)
.
[41] 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 Resour. Res.
55
, 236
(2019)
.
[42] N. R. Leroux, C. B. Marsh, and J. W. Pomeroy, Simulation of preferential flow in snow with a
2D non-equilibrium Richards model and evaluation against laboratory data,
Water Resour. Res.
56
,
e2020WR027466 (2020)
.
123802-15