# A unified kinematic wave theory for melt infiltration into firn

Mohammad Afzal Shadab1,2 Anja Rutishauser3 Cyril Grima2 Marc Andre Hesse4 1
Oden Institute for Computational Engineering and Sciences, University of Texas
at Austin, USA
2 University of Texas Institute for Geophysics, The University of Texas at
Austin, USA
3 Department of Glaciology and Climate, Geological Survey of Denmark and
Greenland, Denmark
4 Department of Earth and Planetary Sciences, Jackson School of Geosciences,
University of Texas at Austin, USA  [

###### Abstract

Motivated by the refreezing of melt water in firn we revisit the one-
dimensional percolation of liquid water and non-reactive gas in porous ice. We
analyze the dynamics of infiltration in the absence of capillary forces and
heat conduction to understand the coupling between advective heat and mass
transport in firn. In this limit, we formulate a kinematic wave theory that
results in a 2$\times$2-system of hyperbolic partial differential equations
(PDEs) corresponding to the conservation of composition and enthalpy. For
simple initial conditions (Riemann problems) this system admits self-similar
solutions that illuminate the structure of melting/refreezing fronts and
analytical solutions are provided for 12 basic cases of physical relevance
encountered in the literature. Further we develop an extended kinematic theory
that encompasses the cases when the firn saturates completely to form a
perched water table governed by elliptic PDE so that the model is no longer
fully hyperbolic (local). These solutions provide benchmarks for numerical
models of melt infiltration into firn. They also provide insight into
important physical processes such as the formation of frozen fringes, the
perching of melt water on pre-existing low porosity layers and the conditions
required for impermeable ice lens formation. Lastly, these analytic solutions
can be utilized to improve and compare the performance of the firn hydrology,
ice-sheet and Earth system models. Our analysis provides a theoretical
framework to understand these important processes in firn which affect the
partitioning between meltwater infiltration and surface runoff and therefore
determine the surface mass loss from ice sheets and its contribution to sea
level rise.

###### keywords:

Kinematic wave theory, infiltration, impermeable ice lens, frozen fringe,
perched water table

Mohammad Afzal<EMAIL_ADDRESS>
##  1 Introduction

Modeling the infiltration of water from both rainfall and surface melting into
firn is still a challenging problem firn2024firn. A better understanding of
this process is needed to improve its representation in the firn hydrology,
ice sheet system and Earth system models clark2017analytical. An understanding
of firn processes can be developed through the analysis of simplified but
physics-based models. Such models are fast and cost-efficient but they only
provide approximate results. Furthermore, the analytical results for synthetic
test problems are useful for the verification and validation of more complex
firn hydrological models that currently show significant deviations
clark2017analytical; vandecrux2020firn. These analytic solutions can also help
to reduce the uncertainties in the existing models and numerical algorithms
such as bucket-schemes to match the field observations by introducing more
physics-based constraints vandecrux2020firn; stevens2020community;
ashmore2020meltwater.

colbeck1971one; colbeck1972theory and colbeck1973water spearheaded the
mathematical modeling of melt infiltration into temperate snow by adapting
models for two-phase flow in porous media to snow and firn Scheidegger1960;
Morel-Seytoux1969. colbeck1974capillary demonstrates that flow due to gravity
is the dominant process while capillary suction/diffusion plays a minor role.
Neglecting capillary diffusion results in a non-linear kinematic wave model
that describes the evolution of the water/melt saturation in the firn
colbeck1972theory. This kinematic wave model results in a hyperbolic partial
differential equation that allows analytic solutions using the method of
characteristics (MOC) lighthill1955kinematic. These analytic solutions capture
the main features of field data from infiltration in the Seward glacier firn
on the St. Elias mountains in Canada sharp1951meltwater.

Later Colbeck extended this kinematic theory to analyze the effects of
impermeable basal boundary colbeck1974water, the retention of water in snow
colbeck1976analysis, studied the effects of layering and heterogeneity on melt
percolation colbeck1979water; Colbeck1991 and determined the permeability of
snow using lysimeter data colbeck1982permeability. Subsequently,
singh1997kinematic use the kinematic theory to explore the interaction of
drying and wetting fronts and then propose the cases of constant, increasing,
and decreasing (drainage) rainfall rates with time (for review see also
singh1997kinematicbook).

colbeck1976analysis also added thermodynamics to the kinematic theory to study
infiltration into cold firn that requires refreezing of the melt at the
wetting front. He shows that the retardation of the wetting front is
relatively minor in snow because the latent heat of fusion is large. Further,
clark2017analytical develop analytic solutions for the kinematic model that
describe the interaction of the wetting and drying fronts in the cold and
temperate snow, following similar solutions for soils Charbeneau1984. They
provide both vertical profiles of water saturation in the snow as well as the
outflow hydrograph and present detailed comparisons with numerical solutions.

Variably saturated flows due to infiltration in firn have not been considered
in the kinematic wave models. The conventional kinematic wave theory for
infiltration fails in a fully-saturated region as the model equation is not
valid anymore. Therefore, the kinematic wave theory needs an extension to
capture fully-saturated regions. Perched water tables may likely form when the
low porosity, initial or refrozen (also known as frozen fringe
miller1972freezing) medium is unable to accommodate the volumetric flux of
water. shadab2022analysis proposes a kinematic wave theory for transitional
infiltration in heterogeneous temperate soil (or firn) which leads to a rising
perched water table. It is a crucial pre-requisite to understand the formation
of frozen fringe or ice lens since the porosity may drop enough to cause
complete saturation.

Recently, there has been a consistent effort to find ice lenses (¡10 cm
thick), layers (10 cm-1 m thick) and slabs (¿1 m thick) in firn from radar
based measurements culberg2021extreme, firn core studies samimi2020meltwater
and satellite data miller2022empirical. These impermeable barriers greatly
affect the surface mass balance by making the lower firn inaccessible to the
meltwater and therefore, increase the runoff especially in the percolation
zone tc-9-1203-2015; Machguth2016; culberg2021extreme. However, there is still
a dearth of literature on understanding the mechanism of their formation and
evolution. A better understanding of the meltwater infiltration process will
help constrain its ultimate partitioning into refreeze, liquid storage, and
runoff meyer2017continuum. Additionally, this improvement will facilitate the
development of more accurate physics-based models and enable comparison with
existing firn hydrology models that currently exhibit significant deviations
vandecrux2020firn. This comparison will be achieved through verification and
validation against analytic solutions for simple benchmark problems.

Therefore in this paper we propose a unified kinematic wave theory for the
infiltration of liquid water and gas in porous ice in the limit of negligible
capillary forces by coupling the mass and enthalpy balances. We extend the
theory to account for the formation and evolution of perched water table where
hyperbolic nature of PDEs fails. It can help construct simple analytical
solutions which can prove to be useful for infiltration in temperate as well
as cold firn where the medium can also form fully saturated regions possibly
due to the formation/existence of a frozen fringe. Finally, we derive
conditions for an impermeable ice lens formation due to meltwater advection
and refreezing. The remainder of this paper is divided into four sections.
Section 2 presents the model formulation. Section 3 considers the problem of
melt transport across a discontinuity and documents twelve nature-inspired
cases with their analytical solutions. Section 4 applies this theory to study
a multilayered firn leading to formation of a perched firn aquifer and
validates it with the numerical solution. Finally, Section 5 concludes the
paper.

##  2 Continuum Model Formulation

We consider a one-dimensional transport problem of three phase, two component
system comprising of water (H2O) in liquid ($w$) and solid ice ($i$) phases
along with non-reactive gas phase ($g$). Ice is considered to be a porous
material when its volume fraction $\phi_{i}\in(0,1)$, where volume fraction
$\phi_{\alpha}$ is the ratio of volume (m3) of the phase
$\alpha\in\\{w,i,g\\}$ to the representative elemental volume (m3). In this
section, we first define the conserved quantities, then introduce the
governing equations and constitutive models, and finally provide the resulting
dimensionless continuum model. The related assumptions will be introduced in
this work when required.

###  2.1 Conserved Quantities

The first conserved variable is the water ($\mathrm{H}_{2}\mathrm{O}$)
composition (kg/m3), $C$, defined as the total mass of water components per
unit representative elemental volume (REV) as

$\displaystyle
C:=\rho_{i}\phi_{i}X_{\mathrm{H}_{2}\mathrm{O},i}+\rho_{w}\phi_{w}X_{\mathrm{H}_{2}\mathrm{O},w}+\rho_{g}\phi_{g}X_{\mathrm{H}_{2}\mathrm{O},g}$
(1)

where $\rho_{\alpha}$ is the density (kg/m3), $\phi_{\alpha}$ refers to volume
fraction of the phase $\alpha\in\\{w,i,g\\}$. The symbol
$X_{\mathrm{H}_{2}\mathrm{O},\alpha}$ represents the mole fraction of
$\mathrm{H}_{2}\mathrm{O}$ component present in phase $\alpha$. Assuming that
water vapor component is negligible in the gas phase
($X_{\mathrm{H}_{2}\mathrm{O},g}=0$) and moreover ice and water phases are
considered pure, i.e.,
$X_{\mathrm{H}_{2}\mathrm{O},i}=X_{\mathrm{H}_{2}\mathrm{O},w}=1$, we can
simplify Equation (1) as

$\displaystyle C=\rho_{i}\phi_{i}+\rho_{w}\phi_{w}.$ (2)

Assuming same density for ice and water, $\rho_{i}\approx\rho_{w}=\rho$,
Equation (2) further simplifies to

$\displaystyle C=\rho(\phi_{i}+\phi_{w}).$ (3)

Next, the volume constraint gives

$\displaystyle\phi_{i}+\phi_{w}+\phi_{g}$ $\displaystyle=1\quad\textrm{and}$
(4) $\displaystyle\Rightarrow C$ $\displaystyle=\rho(1-\phi_{g}).$ (5)

The second conserved variable is related to the thermodynamics of the system.
Since phase change is involved when ice melts, temperature does not accurately
represent all states of the system considered aschwanden2012enthalpy;
alexiades2018mathematical. Therefore the second conserved variable is chosen
to be the enthalpy of the system (J/m3), $H$, which is defined as

$\displaystyle H$
$\displaystyle:=\rho_{i}\phi_{i}h_{i}(T)+\rho_{w}\phi_{w}h_{w}(T)+\rho_{g}\phi_{g}h_{g}(T),$
(6)

where $h_{\alpha}$ is the specific enthalpy (J/kg) of the phase $\alpha$,
which is a piecewise linear function of temperature (K), $T$. For simplicity
we fix the reference enthalpy at the solidus to be $H=0$ where the system is
at the melting temperature, $T=T_{m}$. The specific enthalpy $h_{\alpha}$ of
each phase $\alpha\in\\{w,i,g\\}$ can then be defined as

$\displaystyle h_{i}(T)$
$\displaystyle=\begin{cases}c_{p,i}(T-T_{m}),&T<T_{m}\quad(\textrm{or }H<0)\\\
0,&T\geq T_{m}\quad(\textrm{or }H\geq 0)\end{cases},$ (7) $\displaystyle
h_{w}(T)$ $\displaystyle=\begin{cases}0,&T<T_{m}\quad(\textrm{or }H<0)\\\
c_{p,w}(T-T_{m})+L,&T\geq T_{m}\quad(\textrm{or }H\geq
0)\end{cases}\quad\text{and}$ (8) $\displaystyle h_{g}(T)$
$\displaystyle=c_{p,g}(T-T_{m}).$ (9)

Here $c_{p,\alpha}$ is the specific heat capacity at constant pressure
(J/kg$\cdot$K) for phase $\alpha$, $T_{m}$ is the melting temperature (K) and
$L$ is the latent heat of fusion of water (J/kg). The density and specific
heat capacity of gas are much lower than those of liquid water or ice, i.e.,
$\rho_{g}\ll\rho$ and $c_{p,g}<c_{p,i}$ or $c_{p,w}$ (see Table 1). It helps
make a simplification that the gas phase contribution to the total enthalpy of
the system is negligible. After substituting Equations (7) and (8) in Equation
(6), $H$ can be ultimately formulated as

$\displaystyle H$ $\displaystyle=\begin{cases}\rho
c_{p,i}\phi_{i}(T-T_{m}),&{T<T_{m}}\quad(\textrm{or }H\leq 0)\\\
\rho\phi_{w}L,&{T=T_{m}}\quad(\textrm{or }0<H<CL)\\\
\rho\phi_{w}\left(c_{p,w}(T-T_{m})+L\right),&{T>T_{m}}\quad(\textrm{or }H\geq
CL)\end{cases}.$ (10)

Here the maximum enthalpy limit for the three-phase region is $CL$ because it
is the enthalpy of the system at liquidus (see Figure 1 _a_). The boundaries
of the three-phase region are not included in the region $0<H<CL$ because it
strictly refers to the three phase region. However, no ice (no liquid water)
state can exist at liquidus (solidus), which occurs at $T=T_{m}$. From the
above formulation we classify three regions, where region 1 ($H\leq 0$)
comprises of ice and gas, region 2 ($0<H<CL$) contains all three phases and
region 3 ($H\geq CL$) corresponds to no-matrix state consisting of only water
and gas phases. This classification of regions will be referred to in the
subsequent sections of this paper. In reality, the solid skeleton breaks down
when porosity is less than a critical threshold and the two-phase medium
becomes a slurry of non-contiguous solid grains suspended in liquid
katz2022dynamics.

The temperature and volume fractions of water, ice and gas phases can be
evaluated from composition, $C$, and enthalpy, $H$, as shown in Figures 1 _a_
-1 _d_ respectively. As shown in Equation (5), the volume fraction of gas,
$\phi_{g}$, only depends on the composition, $C$, as shown in Figure 1 _d_.
Next we formulate the governing equations for this model corresponding to the
two conserved variables.

Figure 1: The dependence of temperature and volume fractions on dimensional
and dimensionless enthalpy and composition, ($C,H$) and
($\mathcal{C},\mathcal{H}$) respectively. Dimensional $C,H$: (a) temperature
and volume fractions of (b) water, (c) ice and (d) gas phases. Dimensionless
$\mathcal{C},\mathcal{H}$: (e) scaled temperature and volume fractions of (f)
water, (g) ice and (h) gas phases. The contours are restricted to
$T\in[-100^{\circ}\textrm{C},100^{\circ}\textrm{C}]$ to avoid phase change at
boiling as well as keep the contour levels consistent. Solid black lines are
the level-sets whereas the dashed lines show the boundaries of the regions,
i.e., solidus ($H=0$ or $\mathcal{H}=0$) and liquidus ($H=CL$ or
$\mathcal{H}=\mathcal{C}$).

###  2.2 Transport Model

In a reference frame moving with ice, the conservation equations for water
composition and system enthalpy are respectively given as

$\displaystyle\frac{\partial C}{\partial t}+\nabla\cdot(\textbf{q}\rho
X_{\mathrm{H}_{2}\mathrm{O},w}+\textbf{q}_{g}\rho_{g}X_{\mathrm{H}_{2}\mathrm{O},g}-\phi\rho
D\nabla X_{\mathrm{H}_{2}\mathrm{O},w})$ $\displaystyle=0\quad\textrm{and}$
(11) $\displaystyle\frac{\partial H}{\partial t}+\nabla\cdot(\textbf{q}\rho
h_{w}-\overline{\kappa}\nabla T)$ $\displaystyle=0,$ (12)

where q ($\textbf{q}_{g}$) is the volumetric flux of water (gas) phase
(m3/m${}^{2}\cdot$s) relative to ice phase. Diffusion in the solid has been
neglected in Equation (11). Additionally, $D$ is the effective hydrodynamic
dispersion of the porous medium (m2/s). The effective thermal conductivity of
the mixture (W/m$\cdot$K), $\overline{\kappa}$, is defined as a volume
fraction weighted sum of thermal conductivities (W/m$\cdot$K),
$\kappa_{\alpha}$, of the individual phases $\alpha$ as

$\overline{\kappa}=\kappa_{w}\phi_{w}+\kappa_{i}\phi_{i}+\kappa_{g}\phi_{g},$
(13)

where gas component can again be neglected since gas is relatively a thermal
insulator as $\kappa_{g}\ll\kappa_{w}$ or $\kappa_{i}$ (see Table 1).
Therefore, the effective thermal conductivity can be simplified as

$\overline{\kappa}=\kappa_{w}\phi_{w}+\kappa_{i}\phi_{i}.$ (14)

Since $X_{\mathrm{H}_{2}\mathrm{O},g}=0$ and
$X_{\mathrm{H}_{2}\mathrm{O},w}=1$, Equation (11) can be simplified as

$\displaystyle\frac{\partial C}{\partial t}+\nabla\cdot(\textbf{q}\rho)$
$\displaystyle=0$ (15)

###  2.3 Constitutive Relations

The volumetric flux of water relative to ice, q, can be written using extended
Darcy’s law,

$\displaystyle\textbf{q}=-\frac{\mathrm{k}(\varphi)k_{rw}(s)}{\mu}(\nabla
p-\rho\textbf{g})$ (16)

where $\mathrm{k}$ is the absolute permeability (m2) which is a function of
porosity $\varphi$ $(\varphi=\phi_{w}+\phi_{g}=1-\phi_{i})$, the ratio of void
volume to the bulk volume, $p$ is water pressure (Pa), $\mu$ is the viscosity
of water (Pa$\cdot$s) and g is the acceleration due to gravity vector (m2/s).
The relative permeability for multi-phase flow, $k_{rw}$, display complex
hysteresis blunt2017multiphase, but here we only consider the simplest case
with power law dependence. The relative permeability of water phase, $k_{rw}$,
is a function of the water saturation, $s$, which is the ratio of water phase
volume to void volume $\left(s=\frac{\phi_{w}}{1-\phi_{i}}\right)$. We assume
that the water phase becomes immobile below a certain residual water
saturation, $s_{wr}$. Similarly the gas phase becomes immobile below the
residual gas saturation, $s_{gr}$. As a result, the two-phase fluid flow of
both gas and water phases is restricted to regions where $s_{wr}<s<1-s_{gr}$.
We will refer to regions with $s=1-s_{gr}$ as saturated in the remainder of
this paper. The residual water saturation during drainage comes out to be
$s_{wr}\sim 0.07$ from lysimeter colbeck1976analysis and calorimeter
coleou1998irreducible techniques. However, as the ice is water-wet with a
near-zero contact angle at the ice-water-air interface knight1971experiments,
the residual water saturation during saturation rise (imbibition) is zero due
to hysteresis in the relative permeability-capillary pressure curve
carlson1981simulation; blunt2017multiphase. The phenomenon of hysteresis has
been neglected in the firn hydrology literature but it will affect the speeds
of the meltwater fronts.

Next we assume the problem is gravity dominated in the unsaturated mediums
such that the spatial variations in the difference of fluid phases’ pressures
(capillary pressure) is negligible at the problem length scales
colbeck1972theory. See Smith1983; shadab2022analysis; shadab2022hyperbolic for
a more detailed discussion on neglecting the pressure term in context of soils
using scaling analysis. As a result, the pressure of the water phase in the
unsaturated regions becomes a constant, equal to the reference gas pressure,
i.e., $p=0$ colbeck1972theory; shadab2022analysis. Plugging it in Equation
(16) eliminates the diffusive, pressure term. The volumetric flux of water, q,
finally takes the gravity-driven form

$\displaystyle\textbf{q}=\frac{\mathrm{k}(\varphi)k_{rw}(s)}{\mu}\rho\textbf{g}.$
(17)

The absolute permeability of ice (m2), $\mathrm{k}$, and the relative
permeability of water, $k_{rw}$, are assumed to be power laws kozeny1927uber;
carman1937fluid; brooks1964hydraulic; bear2013dynamics; meyer2017continuum
defined as

$\displaystyle\mathrm{k}(\varphi)$
$\displaystyle=\mathrm{k}_{0}\varphi^{m}=\mathrm{k}_{0}(1-\phi_{i})^{m},$ (18)
$\displaystyle k_{rw}(s)$
$\displaystyle=k^{0}_{rw}s^{n}=k^{0}_{rw}\left(\frac{\phi_{w}}{1-\phi_{i}}\right)^{n},$
(19)

where $\mathrm{k}_{0}$ is a model constant, which can be considered as an
absolute permeability (m2) when there no ice matrix, and $k^{0}_{rw}$ is end
point relative permeability of water phase. Here we have also assumed that the
residual saturations of both water and gas phases are zero, i.e.,
$s_{wr}=s_{gr}=0$. It will provide accurate speeds for the wetting fronts
moving into dry firn, due to hysteresis in the relative permeability. Plugging
Equations (18) and (19) in Equation (17) finally gives

$\displaystyle\textbf{q}(\phi_{i},\phi_{w})=\begin{cases}\textbf{0},&H\leq
0\\\ \frac{\mathrm{k}_{0}k^{0}_{rw}}{\mu}\rho
g(1-\phi_{i})^{m}\left(\frac{\phi_{w}}{1-\phi_{i}}\right)^{n}\hat{\textbf{g}}=K_{h}(1-\phi_{i})^{m}\left(\frac{\phi_{w}}{1-\phi_{i}}\right)^{n}\hat{\textbf{g}},&0<H<CL\textrm{
\& }{C<\rho}\end{cases}$ (20)

where the acceleration due to gravity vector is $\textbf{g}=g\hat{\textbf{g}}$
with $\hat{\textbf{g}}$ being the unit vector in the direction of gravity.
Moreover, 0 is a zero vector in three-dimensions, i.e.,
$\textbf{0}=\left[0,0,0\right]^{T}$. The symbol
$K_{h}=\frac{\mathrm{k}_{0}k^{0}_{rw}}{\mu}\rho g$ is a model constant
equivalent to the saturated hydraulic conductivity of water at unity porosity
(m/s). It is a known constant which can be considered as the maximum gravity-
dependent volumetric flux of water, $|\textbf{q}|$, at unity porosity. Note
that the dynamics at unity porosity is not Darcy-type as the constitutive
relationships (16-18) are only valid for a porous medium. Therefore we will
restrict our analysis to the porous media where $\varphi<1$.

Next we will non-dimensionalize the model in order to make it scale
independent and also find dominant terms which govern the crucial physics of
the problem.

###  2.4 Scaling

The model is scaled using dimensionless variables for composition,
$\mathcal{C}$, enthalpy, $\mathcal{H}$, temperature, $\mathcal{T}$, depth,
$\zeta$, and time, $\tau$, which are defined as

$\mathcal{C}=\frac{C}{\rho},\mathcal{H}=\frac{H}{\rho L},\hskip
2.84526pt\mathcal{T}=\frac{T-T_{m}}{T_{m}},\hskip
2.84526pt\zeta=\frac{z}{\delta},\hskip 2.84526pt\text{ and }\hskip
2.84526pt\tau=\frac{tK_{h}}{\delta}.$ (21)

Here the spatial coordinates (for example, the depth coordinate $z$) are non-
dimensionalized by length scale of heterogeneity or the REV scale of the
problem, $\delta$. Time variable is scaled by the shortest time of water
seepage across the characteristic length through a medium with unity porosity,
i.e., $\delta/K_{h}$. The definitions of conserved quantities $C$ and $H$,
given in Equations (3) and (10) respectively, thus transform into the
dimensionless forms

$\displaystyle\mathcal{C}$ $\displaystyle=\phi_{i}+\phi_{w}=1-\phi_{g},$ (22)
$\displaystyle\mathcal{H}$ $\displaystyle=\begin{cases}\mathcal{C}\hskip
1.42262pt\mathrm{Ste}\hskip
1.42262pt\mathfrak{c}_{p,r}\mathcal{T},&\mathcal{H}\leq 0\\\
\phi_{w},&0<\mathcal{H}<\mathcal{C}\\\ \mathcal{C}(\mathrm{Ste}\hskip
0.85358pt\mathcal{T}+1),&\mathcal{H}\geq\mathcal{C}\end{cases},$ (23)

where Ste is the Stefan number defined as ratio of sensible heat of water at
melting temperature to the latent heat of fusion of
$\mathrm{H}_{2}\mathrm{O}$, i.e., Ste$=c_{p,w}T_{m}/L$ and
$\mathfrak{c}_{p,r}=c_{p,i}/c_{p,w}$ is the ratio of specific heat of ice to
that of water. From the formulations of dimensionless enthalpy (23) and
dimensional specific enthalpy of water phase (8), the dimensionless
temperature, $\mathcal{T}$, and dimensionless specific enthalpy of water
phase, $\mathfrak{h}_{w}=h_{w}/L$, can be derived as

$\displaystyle\mathcal{T}(\mathcal{C},\mathcal{H})$
$\displaystyle=\begin{cases}\frac{\mathcal{H}}{\mathcal{C}\hskip
0.85358pt\textrm{Ste}\hskip 1.42262pt},&{\mathcal{H}\leq 0}\\\
0,&0<\mathcal{H}<\mathcal{C}\\\
\frac{1}{\textrm{Ste}}\left(\frac{\mathcal{H}}{\mathcal{C}}-1\right),&{\mathcal{H}\geq\mathcal{C}}\end{cases},\qquad\textrm{
and
}\quad\mathfrak{h}_{w}(\mathcal{C},\mathcal{H})=\begin{cases}0,&\mathcal{H}\leq
0\\\ 1,&0<\mathcal{H}<\mathcal{C}\\\
\frac{\mathcal{H}}{\mathcal{C}},&\mathcal{H}\geq\mathcal{C}\\\ \end{cases}.$
(24)

Subsequently the volume fractions of the phases, $\phi_{\alpha}$, and the
porosity of the medium, $\varphi$, can be rewritten as functions of
$\mathcal{C}$ and $\mathcal{H}$ as

$\displaystyle\phi_{w}(\mathcal{C},\mathcal{H})$
$\displaystyle=\begin{cases}0,&\mathcal{H}\leq 0\\\
\mathcal{H},&0<\mathcal{H}<\mathcal{C}\\\
\mathcal{C},&\mathcal{H}\geq\mathcal{C}\\\
\end{cases},\quad\phi_{i}(\mathcal{C},\mathcal{H})=\begin{cases}\mathcal{C},&\mathcal{H}\leq
0\\\ \mathcal{C}-\mathcal{H},&0<\mathcal{H}<\mathcal{C}\\\
0,&\mathcal{H}\geq\mathcal{C}\\\
\end{cases},\quad\phi_{g}(\mathcal{C})=1-\mathcal{C}$ (25)
$\displaystyle\textrm{ and }\quad\varphi(\mathcal{C},\mathcal{H})$
$\displaystyle=\begin{cases}1-\mathcal{C},&\mathcal{H}\leq 0\\\
1-\mathcal{C}+\mathcal{H},&0<\mathcal{H}<\mathcal{C}\\\
1,&\mathcal{H}\geq\mathcal{C}\\\ \end{cases}.$ (26)

The scaled temperature and volume fractions of water, ice and gas phases can
be evaluated from dimensionless composition, $\mathcal{C}$, and dimensionless
enthalpy, $\mathcal{H}$, as shown in Figures 1 _e_ -1 _h_ respectively. As
shown in Equation (5), the volume fraction of gas, $\phi_{g}$, only depends on
the dimensionless composition as illustrated in Figure 1 _h_.

The composition and enthalpy transport equations (15 and 12) thus take the
dimensionless form

$\displaystyle\frac{\partial\mathcal{C}}{\partial\tau}+\nabla\cdot\left((1-\phi_{i})^{m}\left(\frac{\phi_{w}}{1-\phi_{i}}\right)^{n}\hat{\textbf{g}}\right)=0,$
(27)
$\displaystyle\frac{\partial\mathcal{H}}{\partial\tau}+\nabla\cdot\left(\mathfrak{h}_{w}(1-\phi_{i})^{m}\left(\frac{\phi_{w}}{1-\phi_{i}}\right)^{n}\hat{\textbf{g}}-(\phi_{w}+\kappa\phi_{i})\frac{\textrm{Ste}}{Pe_{\mathcal{H}}}\nabla\mathcal{T}\right)=0.$
(28)

Here $\kappa=\kappa_{i}/\kappa_{w}$ is the ratio of the heat conductivities of
ice to water and $\alpha_{T}={\kappa_{w}}/{\rho c_{p,w}}$ (m2/s) is the
thermal diffusivity of water. The ratio of heat convected to heat diffused is
defined as the Peclet number for enthalpy equation,
$Pe_{\mathcal{H}}=K_{h}\delta/\alpha_{T}$. Moreover, the divergence and
gradient operators are now scaled with inverse of characteristic depth,
$1/\delta$.

Using thermodynamic and fluid flow parameters from Table 1, the value of
Peclet number, $Pe_{\mathcal{H}}$, comes out to be 3758.6 for $\delta=1$ m.
This value indicates a three orders of magnitude higher heat convection
compared to heat diffusion for gravity-driven infiltration in firn. Therefore,
we can neglect the second order heat conduction term in Equation (28), which
is also the necessary condition in the three phase region
($0<\mathcal{H}<\mathcal{C}$ or $\mathcal{T}=0$) as $\nabla\mathcal{T}=0$.
Assuming local thermodynamic equilibrium (LTE), in the limit
$Pe_{\mathcal{H}}\to\infty$, the system of dimensionless governing equations
(27 and 28) then reduces to quasi-linear system of coupled hyperbolic
equations,

$\displaystyle\frac{\partial\textbf{u}}{\partial\tau}+\nabla\cdot\textit{{f}}(\textbf{u})=\textbf{0}$
(29)

where $\textbf{u}=[\mathcal{C},\mathcal{H}]^{T}$ is the vector of
dimensionless conserved variables and
$\textit{{f}}(\textbf{u})=[\textit{{f}}_{\mathcal{C}},\textit{{f}}_{\mathcal{H}}]^{T}$
is the vector of their corresponding nonlinear flux vectors. Here the flux
vector functions for the dimensionless composition and enthalpy are given as

$\displaystyle\textit{{f}}:=\textit{{f}}_{\mathcal{C}}=\textit{{f}}_{\mathcal{H}}=\begin{cases}\textbf{0},&\mathcal{H}\leq
0,\\\
(1-\mathcal{C}+\mathcal{H})^{m}\left(\frac{\mathcal{H}}{1-\mathcal{C}+\mathcal{H}}\right)^{n}\hat{\textbf{g}},&0<\mathcal{H}<\mathcal{C}\end{cases}.$
(30)

Figure 2: The dimensionless flux of composition or enthalpy in
$\mathcal{C}\mathcal{H}$ phase space for $m=3$ and $n=2$. In region 1
consisting of water and gas region ($\mathcal{H}\leq 0$) as well as region 3
comprising of three phase region ($0<\mathcal{H}<\mathcal{C}$), the flux of
dimensionless enthalpy and composition are identical, i.e.,
$f_{\mathcal{C}}=f_{\mathcal{H}}$. Region 3 with water and gas
($\mathcal{H}\geq\mathcal{C}$) is not considered in the present work.

The above analysis shows the distinct system behaviors in the different
regions based on fluxes (see Figure 2). In this work, region 3 with no solid
matrix is not considered since the constitutive relation for volumetric flux
(Darcy’s law) is not valid anymore. From Equation (30), it can be observed
that the fluxes of enthalpy and composition in the system of governing
equations (29) are identical in regions 1 and 2 defined by the symbol f for
brevity. This is the result of scaling owing to the fact that the composition
changes only when water infiltrates or convects while carrying the enthalpy in
form of latent heat (and specific heat) along with it.

In the next section, we will consider a simple problem of melt transport
across a discontinuity to utilize the method of characteristics (MOC) for
solving the system of hyperbolic partial differential equations
lighthill1955kinematic; leveque1992numerical given in Equation (29).

Table 1: A summary of simplified thermodynamic properties as well as flow properties of water in porous ice used in present work Parameter | Value | Units
---|---|---
$\rho$ | 1000 | kg/m3
$c_{p,w}$ | 4186 | J/(kg K)
$c_{p,i}$ | 2106.1 | J/(kg K)
$\kappa_{w}$ | 0.606 | W/(m K)
$\kappa_{i}$ | 2.25 | W/(m K)
$L$ | 333.55 | kJ/kg
$T_{m}$ | 273.16 | K
$\alpha_{T}$ | 1.45$\times 10^{-7}$ | m2/s
$\mathfrak{c}_{p,r}$ | 0.503 | -
Ste | 3.428 | -
$\mathrm{k}_{0}$ |  $5.56\times 10^{-11}$ meyer2017continuum | m2
$k_{rw}^{0}$ | 1.0 | -
$m$ | 3.0 | -
$n$ | 2.0 | -
$g$ | 9.81 | m/s2
$\mu$ | $10^{-3}$ | Pa s
$K_{h}$ | $5.45\times 10^{-4}$ | m/s

##  3 Melt transport across a discontinuity

Figure 3: One-dimensional Riemann problem: (a) Schematic representation of u
across a discontinuity within or at the boundary of a porous firn. Initial
conditions for the Riemann problem for conserved variables (b) $\mathcal{C}$
and (c) $\mathcal{H}$ plotted against dimensionless depth coordinate, $\zeta$.

This section considers the reaction front arising from melt flow across a
discontinuity in dimensionless composition and enthalpy at a depth, say
$\zeta=0$, as shown in Figure 3. The crucial dynamics of such problems can be
understood using hyperbolic analysis of the coupled system of partial
differential equations in one dimension lighthill1955kinematic;
leveque1992numerical; venkatraman2014analytical; jordan2015reactive;
ghaderi2019theory.

###  3.1 General structure of reaction fronts

Consider the following one-dimensional initial value problem with two constant
states, known as Riemann problem. See leveque1992numerical for a pedagogical
introduction to Riemann problems and their analysis. Let the spatial dimension
be the direction of gravity, $\hat{\textbf{g}}$, which aligns with the
[dimensionless] depth coordinate, $z$ [$\zeta$]. In that case, the
dimensionless flux vector (30) reduces to

$\displaystyle f:=\begin{cases}0,&\mathcal{H}\leq 0\\\
(1-\mathcal{C}+\mathcal{H})^{m}\left(\frac{\mathcal{H}}{1-\mathcal{C}+\mathcal{H}}\right)^{n},&0<\mathcal{H}<\mathcal{C}\end{cases}$
(31)

The system of dimensionless composition and enthalpy conservation equations
(29) can be written in one-dimensional depth coordinates, $\zeta$, as

$\displaystyle\textbf{u}_{\tau}+\textbf{f}(\textbf{u})_{\zeta}=\textbf{0},\quad\tau\in\mathbb{R}^{+},\quad\zeta\in\mathbb{R},$
(32)

with initial conditions

$\displaystyle\textbf{u}=\begin{cases}\textbf{u}_{l},&\zeta<0\\\
\textbf{u}_{r},&\zeta>0\end{cases}$ (33)

where the flux vector in $\zeta$ direction is
$\textbf{f}(\textbf{u})=[f,f]^{T}$ and the subscripts $\tau$ and $\zeta$ refer
to the partial derivatives with respect to the dimensionless time and depth
respectively. The symbol $\mathbb{R}$ refers to the field of real numbers
$(-\infty,\infty)$ whereas $\mathbb{R}^{+}$ refers to the set of positive real
numbers $(0,\infty)$. The subscripts $l$ and $r$ refer to the state of the
system on the left and right sides of a discontinuity. An example of an
initial discontinuity with left state, $\textbf{u}_{l}$, and right state,
$\textbf{u}_{r}$, is shown in Figures 3 _b_ and 3 _c_. The flow is towards the
direction of gravity, assumed to be in $+\zeta$-direction. The solution to the
Riemann problem for well-behaved systems of two coupled nonlinear partial
differential equations is characterized by the formation of an intermediate
state, $\textbf{u}_{i}$, bounded by two waves $\mathscr{W}_{1}$ and
$\mathscr{W}_{2}$ leveque1992numerical. This solution structure, observed in
Figure 4, can be represented as

$\displaystyle\textbf{u}_{l}\xrightarrow{\mathscr{W}_{1}}\textbf{u}_{i}\xrightarrow{\mathscr{W}_{2}}\textbf{u}_{r}$
(34)

Figure 4: Solution of the Riemann problem introduced in Figure 3. (a)
Evolution of dimensionless composition, $\mathcal{C}$, in space, $\zeta$, and
time, $\tau$. (b) The same self-similar solution plotted as a function of
similarity variable $\eta=\zeta/\tau$.

In the context of reactive meltwater transport, the waves $\mathscr{W}_{1}$
and $\mathscr{W}_{2}$ are the reaction fronts and the intermediate state,
$\textbf{u}_{i}$, corresponds to a state between the fronts. The system (32)
can be recasted into a quasilinear form using chain rule as

$\displaystyle{\textbf{u}}_{\tau}+\nabla_{\textbf{u}}\textbf{f}(\textbf{u})\hskip
1.42262pt\textbf{u}_{\zeta}=\textbf{0}$ (35)

where $\nabla_{\textbf{u}}\textbf{f}(\textbf{u})$ is the gradient of flux,
$\textbf{f}(\textbf{u})$, with respect to the conserved variables, u, which
takes the matrix form

$\nabla_{\textbf{u}}\textbf{f}(\textbf{u})=\begin{bmatrix}f_{,\mathcal{C}}&f_{,\mathcal{H}}\\\
f_{,\mathcal{C}}&f_{,\mathcal{H}}\end{bmatrix}$ (36)

where the subscripts $,\mathcal{C}$ and $,\mathcal{H}$ refer to the partial
derivatives with respect to dimensionless composition and enthalpy
respectively. The derivatives of the flux gradient above can be evaluated
explicitly, which are given in Appendix A. The system of advection equations
(35) results in waves (fronts) propagating with their characteristic
velocities. These fronts have self-similar stretching patterns, because of
their own characteristic velocities given by the flux gradient.

###  3.2 Self-similarity of reaction fronts

The recognition of the constant stretching morphology of the reaction fronts
from an initial jump condition allows the introduction of the similarity
variable

$\eta=\frac{\zeta}{\tau}$ (37)

Physically, $\eta$ describes the dimensionless propagation velocity of the
reaction front. The solution generally collapses into a single profile when
plotted as a function of $\eta$ (see Figure 4b). Therefore, the system of
partial differential equations (35) can be transformed into a system of
ordinary differential equations by considering the nonlinear eigenvalue
problem

$\displaystyle(\textbf{A}-\lambda_{p}\textbf{I})~{}\textbf{r}_{p}=\textbf{0},\quad
p\in\\{1,2\\}$ (38)

where the flux gradient is
$\textbf{A}=\nabla_{\textbf{u}}\textbf{f}(\textbf{u})$ and the eigenvector is
$\textbf{r}_{p}=\textrm{d}\textbf{u}/\textrm{d}\eta$ corresponding to the
eigenvalue $\lambda_{p}$. Here the eigenvalues $\lambda_{1}$ and $\lambda_{2}$
are the characteristic propagation speeds of the waves $\mathscr{W}_{1}$ and
$\mathscr{W}_{2}$ respectively. The associated eigenvectors,
$\textbf{r}_{p}=\textrm{d}\textbf{u}/\textrm{d}\eta$, give the pathways
through $\mathcal{C}\mathcal{H}$ plane, also referred to as hodograph plane,
that satisfy the conservation equations (see Figure 5b for example). Any
constant state u that satisfies $\textrm{d}\textbf{u}/\textrm{d}\eta=0$
satisfies the conservation law (35) trivially. Continuous solutions to
Equation (38) satisfy the conservation equations (35).

The opposite case is of the self-sharpening wave front instead of self-
stretching as the information that propagates through the characteristics
crosses each other leveque1992numerical. Here the continuous solution fails
and a moving or a stationary discontinuity appears where the local state
changes abruptly. The propagation velocity of this discontinuity in the
solution satisfies the Rankine-Hugoniot (R-H) jump condition
leveque1992numerical, which is derived from the discrete conservation, in this
case of both mass and enthalpy, around the discontinuity as

$\displaystyle\Lambda_{\mathscr{S}}(\textbf{u}_{+},\textbf{u}_{-})=\frac{\textbf{f}(\textbf{u}_{+})-\textbf{f}(\textbf{u}_{-})}{\textbf{u}_{+}-\textbf{u}_{-}}=\frac{[\textbf{f}(\textbf{u})]}{[\textbf{u}]}$
(39)

where $\Lambda_{\mathscr{S}}$ is the shock speed and $[\,.\,]$ refers to the
jump condition across the shock and subscripts $+$ and $-$ refer to the state
on the left and right sides of a shock wave. Note that the left ($l$) and
right ($r$) states might not necessarily be the left ($-$) and right ($+$)
sides of a shock front, due to the presence of intermediate state(s).

###  3.3 Construction of solution in $\mathcal{C}\mathcal{H}$ hodograph plane

The self-similar solutions are constructed by identifying directions in
$\mathcal{C}\mathcal{H}$ hodograph plane that satisfy conservation laws and
the equation of state. One such direction allows a continuous variation in u,
which can be found by integrating the eigenvectors of the flux gradient.
Another set of direction is determined by the nonlinear algebraic system of
equations arising from R-H jump condition (39), described by shock fronts.
Firstly we consider a system in the disjointed fashion where both left state,
$\textbf{u}_{l}=[\mathcal{C}_{l},\mathcal{H}_{l}]^{T}$, and right state,
$\textbf{u}_{r}=[\mathcal{C}_{r},\mathcal{H}_{r}]^{T}$, reside in the same
region. Then we investigate more complicated cases where left and right states
can lie in different regions. Lastly we discuss the cases where a fully-
saturated region forms, leading to the failure of the current hyperbolic PDE
analysis. This theory is then further extended to analyze the formation and
evolution of fully-saturated regions which are governed by a different,
elliptic PDE shadab2022analysis. Although there can be at most one moving wave
for simple cases where the medium remains unsaturated ($C<\rho$,
$\mathcal{C}<1$ or $\phi_{g}>0$), there can be two moving waves when a fully-
saturated region appears, i.e., $\mathcal{C}=1$.

####  3.3.1 Region 2 only (Three-phase region)

Region 2 ($0<\mathcal{H}<\mathcal{C}$) consists of all three-phases, which is
relevant for temperate glaciers where $\mathcal{T}=0$ or $T=T_{m}$. In three-
phase region, the eigenvalues of the flux gradient (36) are

$\displaystyle\lambda_{1}=0\quad\textrm{and}\quad\lambda_{2}=f_{,\mathcal{C}}+f_{,\mathcal{H}}=n\mathcal{H}^{n-1}(1-\mathcal{C}+\mathcal{H})^{m-n}.$
(40)

where the subscripts ${\mathcal{C}}$ and ${\mathcal{H}}$ refer to the partial
derivatives with respect to the corresponding conserved variable.

Figure 5: Contour plots of (a) second eigenvalue, $\lambda_{2}$, and (b)
propagation velocity of second front with respect to the melt given by
$\phi_{w}\lambda_{2}$ in hodograph plane ($\mathcal{C}\mathcal{H}$ phase
space) for $m=3$ and $n=2$. The slow path $\textbf{r}_{1}$ and fast path
$\textbf{r}_{2}$ are shown with dashed and solid lines respectively in panel
b.

The first reaction front is a stationary contact discontinuity as
$\lambda_{1}=0$. The eigenvalue $\lambda_{2}$ gives a dimensionless
propagation speed of the second reaction front $\mathscr{W}_{2}$ as a function
of $\mathcal{C}$ and $\mathcal{H}$, as plotted in Figure 5 _a_ for $m=3$ and
$n=2$. Therefore, all solutions governed by hyperbolic PDEs (32) will have at
most a single moving reaction front. Due to variable porosity, the
dimensionless system is scaled with respect to the largest saturated hydraulic
conductivity, $K_{h}$, which corresponds to the volumetric flux of water
(Darcy’s flux), rather than the melt velocity. Hence, the propagation speed of
the second reaction front relative to the melt is $\phi_{w}\lambda_{2}$, as
shown in Figure 5 _b_. As expected, the propagation velocity
$\phi_{w}\lambda_{2}$ has significantly retarded near the solidus, as shown by
more separated contours. These eigenvalues yield two corresponding, linearly
independent eigenvectors in $\mathcal{C}\mathcal{H}$ hodograph plane given by

$\displaystyle\textbf{r}_{1}=\begin{bmatrix}\frac{-f_{,\mathcal{H}}}{f_{,\mathcal{C}}}\\\
1\end{bmatrix}=\begin{bmatrix}\frac{(n(1-\mathcal{C})+m\mathcal{H})}{(m-n)\mathcal{H}}\\\
1\end{bmatrix}\quad\textrm{and}\quad\textbf{r}_{2}=\begin{bmatrix}\frac{f_{,\mathcal{H}}}{f_{,\mathcal{H}}}\\\
1\end{bmatrix}=\begin{bmatrix}1\\\ 1\end{bmatrix}.$ (41)

These eigenvectors can be used to find the integral curves using the system of
ODEs

$\displaystyle\frac{\textrm{d}\textbf{u}}{\textrm{d}\eta}$
$\displaystyle=\frac{1}{\nabla_{\textbf{u}}\lambda_{p}\cdot\textbf{r}_{p}}\textbf{r}_{p}$
(42)

which can be further integrated to obtain the solution pathways
$\textbf{u}(\textbf{u}_{0},\eta)$ as

$\displaystyle\textbf{u}(\textbf{u}_{0},\eta)$
$\displaystyle=\textbf{u}_{0}+\int_{\lambda_{p}(\textbf{u}_{0})}^{\eta}\frac{1}{\nabla_{\textbf{u}}\lambda_{p}\cdot\textbf{r}_{p}}\textbf{r}_{p}~{}\textrm{d}\eta^{\prime}.$
(43)

These paths in $\mathcal{C}\mathcal{H}$ hodograph plane comprise the set of
states that can be connected to the initial state $\textbf{u}_{0}$ by a
reaction front with a continuous variation in u. In three phase region, the
family of integral curves corresponding to first eigenvector $\textbf{r}_{1}$
is referred to as slow path as $\lambda_{1}<\lambda_{2}$ and is given by

$\displaystyle\mathcal{C}$
$\displaystyle=1+\mathcal{H}+\mathfrak{C}\mathcal{H}^{-\frac{n}{m-n}}.$ (44)

where $\mathfrak{C}$ is the constant of integration, which can be found for
the initial point $\textbf{u}_{0}$. The slow path lines are the same as
constant flux lines, as shown in Figure 5b. Next, the family of integral
curves corresponding to the second eigenvector is known as fast path and is
given by

$\displaystyle\mathcal{C}$ $\displaystyle=\mathcal{H}+\mathfrak{C}.$ (45)

The speed of second characteristic $\lambda_{2}$ is non-negative and increases
monotonically in the direction of integral curves corresponding to the second
eigenvector $\textbf{r}_{2}$ (fast paths) in three-phase region for $n>1$ (see
Lemma Appendix 1.1 and Figure 5 _b_). Additionally, the slow path corresponds
to constant flux contours as $\lambda_{1}=0$ and the fast path corresponds to
constant porosity, $\varphi$, contours (Lemma Appendix 1.2).

#### Solutions in three-phase region

We will now discuss the different scenarios involving three-phase region only.
The discussion below assumes that $\textbf{u}_{0}$ is the left state,
$\textbf{u}_{l}$, and describes the set of permissible right states
$\textbf{u}_{r}$.

##### (a.) Stationary linear reaction front (Case I):

The integral curves associated with $\lambda_{1}$ and $\textbf{r}_{1}$, known
as the first characteristic field ($\lambda_{1},\textbf{r}_{1}$), are constant
flux lines in the $\mathcal{C}\mathcal{H}$ phase space (or hodograph plane).
Any right state $\textbf{u}_{r}$, along the integral curve, connected to
$\textbf{u}_{0}$ by a stationary discontinuity is a weak solution of Equation
(35). Because $\lambda_{1}=0$, the first wave is a stationary contact
discontinuity $\mathscr{C}_{1}$. The fluxes of $\mathcal{C}$ and $\mathcal{H}$
on both sides are the same so that the melt transport does not change
$\mathcal{C}$ and $\mathcal{H}$ and the front does not evolve. So, a contact
discontinuity is the solution for the left and right states lying on the slow
path (44) (constant flux lines) satisfying

$\displaystyle\frac{\varphi_{l}}{\varphi_{r}}$
$\displaystyle=\frac{1-\mathcal{C}_{l}+\mathcal{H}_{l}}{1-\mathcal{C}_{r}+\mathcal{H}_{r}}=\left(\frac{\mathcal{H}_{l}}{\mathcal{H}_{r}}\right)^{-\frac{n}{m-n}}\quad\textrm{or}\quad
f(\textbf{u}_{l})=f(\textbf{u}_{r}).$ (46)

The final solution in this case takes the form

u $\displaystyle=\begin{cases}\textbf{u}_{l},\quad&\zeta<0\\\
\textbf{u}_{r},\quad&\zeta>0\end{cases}.$ (47)

Figures 6a-c (green color) illustrate an example of a system that results in a
stationary contact discontinuity. This system corresponds to a steady
meltwater flux of $f=0.112$ inside temperate porous firn with a jump in
porosity from 70% at $\zeta<0$ to 55.3% in $\zeta>0$ (coarse-to-fine
transition in firn) leading to liquid water contents of $0.4$ and $0.45$
respectively in these regions. These values correspond to
$\textbf{u}_{l}=[\mathcal{C}_{l},\mathcal{H}_{l}]^{T}=[0.7,0.4]^{T}$ and right
state $\textbf{u}_{r}=[0.897,0.45]^{T}$. The resulting system illustrated by
the volume fractions of the three phases highlights a porosity jump in the
temperate firn that has a constant steady meltwater flux on both sides, as
shown in Figures 6d-f at different times. This case is summarized in Table 2.

Figure 6: The simple solutions of Riemann problem leading to a contact
discontinuity $\mathscr{C}$ (green), rarefaction $\mathscr{R}$ (blue), and
shock waves $\mathscr{S}$ (red). (a) Construction of solution in hodograph
plane and their corresponding self-similar analytical solutions for (b)
dimensionless composition and (c) dimensionless enthalpy with dimensionless
velocity $\eta$. The evolution of the volume fractions of the three phases in
the system for different configurations at dimensionless times
$\tau=0,0.5,1.0$: (d-f) Case I - contact discontinuity, (g-i) Case II - drying
front/rarefaction wave, and (j-l) Case III - wetting front/shock wave.

##### (b.) Moving nonlinear reaction front:

The integral curves associated with the second characteristic field
($\lambda_{2},\textbf{r}_{2}$) are the constant porosity $\varphi$ contours.
Thus, a nonlinear characteristic wave is the solution for the left and right
states lying on the fast path (45) satisfying

$\displaystyle\mathcal{C}_{r}$
$\displaystyle=\mathcal{C}_{l}+(\mathcal{H}_{r}-\mathcal{H}_{l}).$ (48)

###### (i.) Rarefaction wave (Case II):

If the characteristic speed $\lambda_{2}$ varies smoothly from left to right
state, any right state u along the fast path is connected to left state
$\textbf{u}_{0}$ by a continuously varying reaction front (Figure 6a, blue
lines). The propagation velocity, $\lambda_{2}$, along these continuous
reaction fronts increases monotonically such that the reaction front spreads
with time. These self-smoothening drying fronts are referred to as rarefaction
waves, denoted by the symbol $\mathscr{R}$. The term drying front refers to
drying due to meltwater drainage instead of refreezing clark2017analytical.
Rarefaction waves are a weak solution of Equation (38) if the resultant
profile of u is single-valued. This condition is satisfied if u lies on the
branch of the integral curve $\textbf{r}_{2}$ emanating from $\textbf{u}_{0}$
in the direction of increasing $\lambda_{2}$ (see Figures 5 and 6a, blue
lines). The analytical solution concerning the self-similar variable $\eta$
for a rarefaction wave on an integral curve can be evaluated from Equation
(43), which comes out to be

$\displaystyle\mathcal{H}$
$\displaystyle=\sqrt[n-1]{\frac{\eta}{n(1-\mathcal{C}_{l}+\mathcal{H}_{l})^{m-n}}}\quad\text{and}$
(49) $\displaystyle\mathcal{C}$
$\displaystyle=\mathcal{C}_{l}-\mathcal{H}_{l}+\sqrt[n-1]{\frac{\eta}{n(1-\mathcal{C}_{l}+\mathcal{H}_{l})^{m-n}}}.$
(50)

The final solution in this case takes the form

u
$\displaystyle=\begin{cases}\textbf{u}_{l},\quad&\eta<\lambda_{2}(\textbf{u}_{l})\\\
\begin{bmatrix}\mathcal{C}_{l}-\mathcal{H}_{l}+\sqrt[n-1]{\frac{\eta}{n(1-\mathcal{C}_{l}+\mathcal{H}_{l})^{m-n}}}\\\
\sqrt[n-1]{\frac{\eta}{n(1-\mathcal{C}_{l}+\mathcal{H}_{l})^{m-n}}}\end{bmatrix},\quad&\lambda_{2}(\textbf{u}_{l})<\eta<\lambda_{2}(\textbf{u}_{r})\\\
\textbf{u}_{r},\quad&\eta>\lambda_{2}(\textbf{u}_{r})\end{cases}.$ (51)

where the speed of the second characteristic $\lambda_{2}(\cdot)$ is evaluated
from Equation (40). An example of this case is when meltwater flux instantly
drops on a glacier, leading to a smoothening drainage front. Figures 6a-c
(blue line) show a moving rarefaction developed inside a 70% porous firn due
to an instantaneous reduction in meltwater flux, captured by lower (40%) and
higher (55%) liquid water content layers on top and bottom respectively. These
numbers translate to left (top) $\textbf{u}_{l}=[0.7,0.4]^{T}$ and right
(bottom) states being $\textbf{u}_{r}=[0.85,0.55]^{T}$. The resulting
evolution of volume fractions is shown in Figures 6g-i. The leading edge of
the rarefaction front moves faster than the trailing edge connected by a
gradual, linear smoothening. In this case, it is a linear profile (straight
line) as the power law exponents in Equation (51) are $m=3$ and $n=2$.
Drying/rarefaction front has been studied by clark2017analytical and was
observed in models studying the Dye-2 site in Greenland on 12 August 2016
after the meltwater flux ceased samimi2020meltwater; vandecrux2020firn;
colliander2022ice.

###### (ii.) Shock wave (Case III):

If the right state u lies on the opposite branch of the integral curve (shown
by red line in Figure 6a for example), a continuous reaction front would
result in unphysical solutions as the characteristics will cross each other.
Therefore, in this case u is connected to left state $\textbf{u}_{0}$ by a
discontinuous reaction front that propagates with a velocity,
$\Lambda_{\mathscr{S}}(\textbf{u}_{0},\textbf{u})$, which can be calculated
from R-H jump condition (39) using initial conditions
($\textbf{u}_{l},\textbf{u}_{r}$). Such fronts are referred to as
wetting/shock fronts, denoted by the symbol $\mathscr{S}$. The set of
permissible right states u that can be connected to the left state
$\textbf{u}_{0}$ by shocks lie on the segment of the Hugoniot-locus that
satisfies the entropy condition. In the three-phase region for the system of
equations considered, the Hugoniot-locus is the same as the integral curve,
which is found to be the fast path $\textbf{r}_{2}$ (45) from the Hugoniot
jump condition (39), since the flux of enthalpy and composition are the same.
The dimensionless velocity of the shock from Equation (39) is then

$\displaystyle\Lambda_{\mathscr{S}}(\textbf{u}_{+},\textbf{u}_{-})=\frac{\textrm{d}\zeta}{\textrm{d}\tau}=\frac{f(\textbf{u}_{+})-f(\textbf{u}_{-})}{\mathcal{H}_{+}-\mathcal{H}_{-}}=\frac{f(\textbf{u}_{+})-f(\textbf{u}_{-})}{\mathcal{C}_{+}-\mathcal{C}_{-}},$
(52)

where the symbol $-$ is the left (top) state and $+$ is the right (bottom)
state for this particular configuration where
$f(\textbf{u}_{l})>f(\textbf{u}_{r})$. The final solution thus takes the form

u
$\displaystyle=\begin{cases}\textbf{u}_{l},\quad&\zeta/\tau<\Lambda_{\mathscr{S}}(\textbf{u}_{r},\textbf{u}_{l})\\\
\textbf{u}_{r},\quad&\zeta/\tau>\Lambda_{\mathscr{S}}(\textbf{u}_{r},\textbf{u}_{l})\end{cases}.$
(53)

An example of this case is when meltwater flux instantly increases on a
temperate glacier, it leads to a sharp wetting front propagating in the
direction of gravity. Figures 6j-k show a moving shock front developed inside
a 70% porous firn due to an instantaneous increase in meltwater flux, captured
by higher (40%) and lower (20%) liquid water content $\phi_{w}$ in top and
bottom layers respectively. These numbers translate to left (top) state being
$\textbf{u}_{l}=[0.7,0.4]^{T}$ and right (bottom) state being
$\textbf{u}_{r}=[0.5,0.2]^{T}$ sketched analytically in Figures 6a-c. The
wetting front has been discussed (colbeck1972theory; humphrey2012thermal;
meyer2017continuum; clark2017analytical, e.g.) and observed in models
(vandecrux2020firn; samimi2021time; colliander2022ice, e.g.) and field
observations including at the Dye-2 site in Greenland on 9 August 2016 when
meltwater percolates in temperate firn heilig2018seasonal;
samimi2020meltwater.

##### (c.) Two fronts with an intermediate state:

The solution profile contains a single reaction front if $\textbf{u}_{l}$ and
$\textbf{u}_{r}$ share the same integral curve or Hugoniot-locus. In all other
cases in the three-phase region without complete saturation, a different state
than left or right state forms which is referred to as an intermediate state,
$\textbf{u}_{i}$, in the hodograph plane. At this intermediate state, the
solution switches from the first characteristic field
$(\lambda_{1},\textbf{r}_{1})$ to the second $(\lambda_{2},\textbf{r}_{2})$.
In other words, at the intermediate state $\textbf{u}_{i}$, the solution
changes from the stationary front along constant flux lines (slow path) to the
advancing reaction front along the path of constant porosity contours (fast
path) that also mimic the liquidus ($\mathcal{H}=\mathcal{C}$). The two
possible intermediate states are given by the intersections of the integral
curves emanating from $\textbf{u}_{l}$ and $\textbf{u}_{r}$ (see Figure 7a for
example). Only one intersection yields a physically realistic single-value
solution. The correct intersection is selected by requiring that the
propagation speed increases monotonically from $\textbf{u}_{l}$ and
$\textbf{u}_{r}$. A single-valued solution is ensured if and only if
$\textbf{u}_{l}$ and $\textbf{u}_{i}$ lie on the slow path first and then
$\textbf{u}_{i}$ is connected to $\textbf{u}_{r}$ along the fast path. The
reactive melt transport system considered here only allows two solutions for
this case:

$\displaystyle\textbf{u}_{l}\xrightarrow{\mathscr{C}_{1}}\textbf{u}_{i}\xrightarrow{\mathscr{R}_{2}}\textbf{u}_{r}\quad\textrm{and}\quad\textbf{u}_{l}\xrightarrow{\mathscr{C}_{1}}\textbf{u}_{i}\xrightarrow{\mathscr{S}_{2}}\textbf{u}_{r},$
(54)

because the first characteristic is linearly degenerate and the reaction front
along the slow path is always a contact discontinuity $\mathscr{C}_{1}$. The
reactive melt transport across an initial discontinuity is characterized by
the formation of a reacted zone corresponding to $\textbf{u}_{i}$ that is
bounded between a stationary front $\mathscr{C}_{1}$, and an advancing front
that is either a rarefaction wave $\mathscr{R}_{2}$ or a shock wave
$\mathscr{S}_{2}$. Below we will discuss these two cases.

Figure 7: Formation of an intermediate state $\textbf{u}_{i}$ or
$\textbf{u}^{*}_{i}$ for Case IV - $\mathscr{C}_{1}\mathscr{R}_{2}$ or Case V
- $\mathscr{C}_{1}\mathscr{S}_{2}$ respectively. An asterisk is used to
differentiate the two intermediate states corresponding to the two cases. (a)
Construction of solution in hodograph plane and their corresponding self-
similar analytical solutions for (b) dimensionless composition and (c)
dimensionless enthalpy with dimensionless velocity $\eta$. The evolution of
the volume fractions of the three phases in the system at dimensionless times
$\tau=0,0.5,1.0$ for the two configurations: (d-f) Case IV -
$\mathscr{C}_{1}\mathscr{R}_{2}$ and (g-i) Case V -
$\mathscr{C}_{1}\mathscr{S}_{2}$.

###### (i.) 1-Contact discontinuity and 2-Rarefaction (Case IV):

In this case $f(\textbf{u}_{l})<f(\textbf{u}_{r})$ and the resulting first
characteristic wave is a contact discontinuity, $\mathscr{C}_{1}$, which
satisfies Equation (44) for left state $\textbf{u}_{l}$ and the intermediate
state $\textbf{u}_{i}$. The second characteristic wave is a rarefaction which
is governed by Equations (45), (49) and (50) for intermediate state
$\textbf{u}_{i}$ and right state $\textbf{u}_{r}$. Combining all these
equations results in a nonlinear algebraic equation to evaluate
$\mathcal{H}_{i}$ that is given by

$\displaystyle{1+\mathcal{H}_{r}-\mathcal{C}_{r}-\mathcal{H}_{i}^{-\frac{1}{n-1}}{\mathcal{H}_{l}}^{\frac{n}{n-1}}{\left(\frac{1-\mathcal{C}_{l}+\mathcal{H}_{l}}{1-\mathcal{C}_{r}+\mathcal{H}_{r}}\right)^{\frac{m-n}{n-1}}}+\mathcal{H}_{i}}$
$\displaystyle=\left({1+\mathcal{H}_{l}-\mathcal{C}_{l}}\right)\left(\frac{\mathcal{H}_{i}}{\mathcal{H}_{l}}\right)^{-\frac{n}{m-n}}$
(55)

which can be re-written in terms of porosities, $\varphi$, as

$\displaystyle{\varphi_{r}-\mathcal{H}_{i}^{-\frac{1}{n-1}}{\mathcal{H}_{l}}^{\frac{n}{n-1}}{\left(\frac{\varphi_{l}}{\varphi_{r}}\right)^{\frac{m-n}{n-1}}}+\mathcal{H}_{i}}$
$\displaystyle={\varphi_{l}}\left(\frac{\mathcal{H}_{i}}{\mathcal{H}_{l}}\right)^{-\frac{n}{m-n}}.$
(56)

Next, the composition at intermediate state, $\mathcal{C}_{i}$, can be
computed from Equations (45) and (40) which corresponds to the fast path and
the speed of the second characteristic, $\lambda_{2}$, respectively. The final
solution in this case takes the form

u $\displaystyle=\begin{cases}\textbf{u}_{l},\quad&\zeta<0\\\
\textbf{u}_{i},\quad&0<\zeta/\tau<\lambda_{2}(\textbf{u}_{i})\\\
\begin{bmatrix}\mathcal{C}_{i}-\mathcal{H}_{i}+\sqrt[n-1]{\frac{\eta}{n(1-\mathcal{C}_{i}+\mathcal{H}_{i})^{m-n}}}\\\
\sqrt[n-1]{\frac{\eta}{n(1-\mathcal{C}_{i}+\mathcal{H}_{i})^{m-n}}}\end{bmatrix},\quad&\lambda_{2}(\textbf{u}_{i})<\eta<\lambda_{2}(\textbf{u}_{r})\\\
\textbf{u}_{r},\quad&\eta>\lambda_{2}(\textbf{u}_{r})\end{cases}.$ (57)

An example of this case is a sudden drop in meltwater flux inside temperate
firn where porosity also reduces with depth (Figure 7d). Physically, the first
contact discontinuity represents the constant flux of water entering the
bottom, low porosity layer (Figures 7e,f). The second rarefaction shows the
drainage of the wetter firn due to gravity. Figure 7 _a_ shows the
construction of the solution $\mathscr{C}_{1}\mathscr{R}_{2}$ in the blue line
for the left (top) state which is more porous $\varphi_{l}=80\%$, and has less
water content (LWC) $\phi_{w,l}=0.1$. The right (bottom) state is less porous
$\varphi_{r}=58.0\%$ but has more liquid water content $\phi_{w,r}=0.528$.
These values correspond to left and right states being
$\textbf{u}_{l}=(0.3,0.1)^{T}$ and $\textbf{u}_{r}=(0.948,0.528)^{T}$
respectively shown in Figures 7a-c. The intermediate state comes out to be
$\textbf{u}_{i}=(0.538,0.117)^{T}$ which corresponds to LWC
$\phi_{w,i}=11.7\%$ and the same porosity as the right state, i.e.,
$\varphi_{i}=58.0\%$. Figures 7 _b_ and 7 _c_ show the corresponding self-
similar analytical solutions for composition and enthalpy respectively for
this case with blue lines which only depend on the dimensionless velocity. The
rarefaction moves down with a characteristic velocity that can be computed
analytically. Figures 7 _d_ -_f_ shows the resulting evolution of volume
fraction of each phase at different times showing self-similar expansion of
the rarefaction wave.

###### (ii.) 1-Contact discontinuity and 2-Shock (Case V):

This case is similar to Case IV but has $f(\textbf{u}_{l})>f(\textbf{u}_{r})$.
As a result, the first characteristic wave is a contact discontinuity,
$\mathscr{C}_{1}$, that satisfies Equation (44) for the left state
$\textbf{u}_{l}$ and intermediate state $\textbf{u}_{i}$. Since the second
characteristic lies on the Hugoniot locus, the result is a shockwave,
$\mathscr{S}_{2}$, which satisfies the Hugoniot-jump condition (39) for the
intermediate state and right state. Combining Equations (44) and (45) gives a
simple relation for dimensionless enthalpy at intermediate state,
$\mathcal{H}_{i}$, to be

$\displaystyle\mathcal{H}_{i}$
$\displaystyle=\mathcal{H}_{l}{\left(\frac{\varphi_{l}}{\varphi_{r}}\right)}^{\frac{m-n}{n}}=\mathcal{H}_{l}{\left(\frac{1+\mathcal{H}_{l}-\mathcal{C}_{l}}{1+\mathcal{H}_{r}-\mathcal{C}_{r}}\right)}^{\frac{m-n}{n}}.$
(58)

Then Equation (45) for the intermediate state, $\textbf{u}_{i}$, and right
state, $\textbf{u}_{r}$, provides the value of dimensionless composition at
intermediate state, $\mathcal{C}_{i}$. The final solution in this case takes
the form

u $\displaystyle=\begin{cases}\textbf{u}_{l},\quad&\zeta<0\\\
\textbf{u}_{i},\quad&0<\zeta/\tau<\Lambda_{\mathscr{S}}(\textbf{u}_{i},\textbf{u}_{r})\\\
\textbf{u}_{r},\quad&\zeta/\tau>\Lambda_{\mathscr{S}}(\textbf{u}_{i},\textbf{u}_{r})\end{cases}.$
(59)

An example of this case is a sudden rise in meltwater flux inside temperate
firn where porosity also reduces with depth (Figure 7g). Physically, the first
contact discontinuity represents the increased, constant flux of water
entering the bottom, low porosity layer. The second shock shows the wetting
front advancing the water content to dryer firn because of gravity (Figures
7h,i). Figure 7 _a_ shows the construction of the solution
$\mathscr{C}_{1}\mathscr{S}_{2}$ with red line for the left state which is
less porous, i.e., $\varphi_{l}=58.0\%$, but has more water content
($\phi_{w,l}=52.8\%$), and the right state is more porous and less wet
corresponding to $\varphi_{r}=80\%$ and LWC $\phi_{w,r}=10\%$. These values
correspond to $\textbf{u}_{l}=(0.948,0.528)^{T}$ and
$\textbf{u}_{r}=(0.3,0.1)^{T}$ in Figures 7a-c. The intermediate state comes
out to be $\textbf{u}^{*}_{i}=(0.648,0.448)^{T}$ which corresponds to the
wetter intermediate region behind the wetting front with LWC
$\phi_{w,i}=44.8\%$ but the same porosity as the right state, i.e.,
$\varphi_{i}=80\%$ (see Figures 7 _g_ -_i_). Figures 7 _b_ and 7 _c_ show the
corresponding self-similar analytical solutions for composition and enthalpy
respectively for this case with red lines which only depend on the
dimensionless velocity.

##### (d.) Two fronts and a jump with two intermediate states (formation of a
saturated region, Case VI):

The initial conditions of this case are similar to Case V
($\mathscr{C}_{1}\mathscr{S}_{2}$). However, in this case, the slow path
emanating from the left state does not intersect with the fast path from the
right state in the three-phase region where $\mathcal{C}<1$. In other words,
the bottom layer is unable to accommodate the flux of meltwater from the top
layer, as discussed for soils in shadab2022analysis. If the intermediate state
leads to complete saturation, i.e., $\mathcal{C}_{i}=1$, then the proposed
hyperbolic PDE solution framework breaks down. This happens because inside the
fully-saturated region, the dynamics of the water phase change from gravity-
driven to pressure-driven as the governing model changes from hyperbolic
(local) to elliptic (global) partial differential equation (see
shadab2022analysis for a detailed analysis). In this case, the dynamics
becomes complicated as the complete solution cannot be directly interpreted
from the hodograph plane because the flux in the saturated region may not
simply be the hydraulic conductivity (30) anymore.

Nevertheless, we can still construct a full analytical solution to this
problem using the extended kinematic wave approximation proposed by
shadab2022analysis. In this case, the solution consists of three waves
including a backfilling shock moving upwards, denoted by symbol
$\mathscr{S}_{1}^{*}$, a stationary jump at the initial location of the jump,
denoted by $\mathscr{J}_{2}$, and a downward moving shock (wetting front),
$\mathscr{S}_{3}$, into the less porous and temperate layer. Note that the
jump $\mathscr{J}_{2}$ lies in the saturated region and therefore does not
represent any hyperbolic wave. Therefore, it is represented by a broken arrow
in the full solution given by

$\displaystyle\textbf{u}_{l}\xrightarrow{\mathscr{S}^{*}_{1}}\textbf{u}_{i_{1}}\mathrel{\mathchoice{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar
0\relax\symAMSa
04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar
0\relax\symAMSa
04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptscriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar
0\relax\symAMSa
04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar
0\relax\symAMSa
04B}{}}}\textbf{u}_{i_{2}}\xrightarrow{\mathscr{S}_{3}}\textbf{u}_{r}.$

The solution is explained in detail as follows. First, the backfilling shock
on the fast path (constant porosity line) connects to the first intermediate
state, $i_{1}$, which lies on the line $\mathcal{C}=1$. Therefore, the first
intermediate state variables are

$\displaystyle\mathcal{C}_{i_{1}}=1\quad\text{and}\quad\mathcal{H}_{i_{1}}=1-\mathcal{C}_{l}+\mathcal{H}_{l}.$
(60)

This state is observed right next to the left state and can be considered as
the rising perched water table in the region $\zeta<0$. The speed of this
backfilling, upper shock in dimensionless form is again given by Rankine-
Hugoniot jump condition leveque1992numerical as

$\displaystyle\Lambda_{\mathscr{S}^{*}_{1}}=\frac{\textrm{d}\zeta_{U}}{\textrm{d}\tau}=\frac{f(\textbf{u}_{l})-q_{s}(\tau)}{\mathcal{H}_{l}-\mathcal{H}_{i_{1}}}=\frac{f(\textbf{u}_{l})-q_{s}(\tau)}{\mathcal{C}_{l}-\mathcal{C}_{i_{1}}}<0,$
(61)

where $\zeta_{U}$ is the location of the upper shock. The shock moves upwards
due to choking as the numerator is positive since flux $f(\textbf{u}_{l})$ is
more than the time-dependent dimensionless flux in the saturated region,
$q_{s}(\tau)$, which is also scaled by $K_{h}$. Note that the flux in the
saturated region $q_{s}(\tau)$ is not the saturated hydraulic conductivity but
instead, it is governed by the dynamics of the saturated region.

Next, the first intermediate state $i_{1}$ is connected to the second
intermediate state $i_{2}$ through a stationary jump $\mathscr{J}_{2}$ at the
location of initial jump at $\zeta=0$. Both intermediate states lie in the
fully-saturated region and therefore the jump $\mathscr{J}_{2}$ does not
represent a hyperbolic wave. The flux inside the isothermal saturated region
$q_{s}(\tau)$ is found to be uniform in this case shadab2022analysis.
Therefore, the flux between the two intermediate states is also same, equal to
$q_{s}(\tau)$. The state variables $\mathcal{C}$ and $\mathcal{H}$ for the
second intermediate state, $i_{2}$, are provided by the right state. The
second intermediate state also lies at $\mathcal{C}=1$ on the fast path (45)
(Hugoniot locus) emanating from right state, i.e.,
$\mathcal{C}_{r}-\mathcal{H}_{r}=\mathcal{C}_{i_{2}}-\mathcal{H}_{i_{2}}$.
Therefore, the second intermediate state variables are simply,

$\displaystyle\mathcal{C}_{i_{2}}=1\quad\text{and}\quad\mathcal{H}_{i_{2}}=1-\mathcal{C}_{r}+\mathcal{H}_{r}.$
(62)

Similarly, the velocity of the downward-moving lower shock (wetting front) is

$\displaystyle\Lambda_{\mathscr{S}_{3}}=\frac{\textrm{d}\zeta_{L}}{\textrm{d}\tau}=\frac{q_{s}(\tau)-{f{(\textbf{u}_{r}})}}{\mathcal{C}_{i_{2}}-\mathcal{C}_{r}}=\frac{q_{s}(\tau)-f{(\textbf{u}_{r}})}{1-\mathcal{C}_{r}}>0$
(63)

where $\zeta_{L}$ is the dimensionless location of the lower shock. Similar to
a two-layer soil shadab2022analysis, the dimensionless flux in the saturated
region is a depth-based harmonic mean of the dimensionless saturated hydraulic
conductivities at the two intermediate states given by

$\displaystyle
q_{s}(\tau)=\frac{\zeta_{U}(\tau)-\zeta_{L}(\tau)}{\frac{\zeta_{U}(\tau)}{K_{i_{1}}}-\frac{\zeta_{L}(\tau)}{K_{i_{2}}}},$
(64)

where $K_{i_{1}}$ and $K_{i_{2}}$ are the dimensionless saturated hydraulic
conductivities at the first and second intermediate states given by

$\displaystyle
K_{i_{1}}=\varphi_{i_{1}}^{m}=(1-\mathcal{C}_{l}+\mathcal{H}_{l})^{m}\quad\text{and}\quad
K_{i_{2}}=\varphi_{i_{2}}^{m}=(1-\mathcal{C}_{r}+\mathcal{H}_{r})^{m}$ (65)

using Equations (60) and (62). Note that the porosities at the first and
second intermediate states are the porosities of the left and right state
respectively. Solving the system of coupled ordinary differential equations
(63) & (61) along with the definition of flux in the saturated region (64)
gives the location of the shocks and the flux in the saturated region. Similar
to the case of two-layered soils in shadab2022analysis, to find the analytic
value of $q_{s}(\tau)$, the ratio of shock speeds can be considered a constant
as an ansatz given by

$\displaystyle\frac{\Lambda_{{\mathscr{S}^{*}_{1}}}}{\Lambda_{{\mathscr{S}_{3}}}}=\frac{\zeta_{U}(\tau)}{\zeta_{L}(\tau)}=\mathfrak{R}\quad\mbox{for}\quad
0<\tau\leq\tau_{p}.$ (66)

where $\mathfrak{R}$ is the constant ratio of shock speeds which is negative
and $\tau_{p}$ is the dimensionless time of ponding when the upward moving
shock reaches to the surface. In the special case when this jump condition
happens to exist at the surface, $\mathfrak{R}=0$, $\tau_{p}=0$ and
$q_{s}=K_{i_{2}}$ which is saturated hydraulic conductivity of the second
intermediate state. Otherwise, by substituting equations from the shock speed
ratio definition (66), shock speeds (61) and (63), and flux in the saturated
region (64), $\mathfrak{R}$ comes out to be the solution of a quadratic
equation

$\displaystyle\mathfrak{R}(\textbf{u}_{l},\textbf{u}_{r})$
$\displaystyle=\frac{-b-\sqrt{b^{2}-4ac}}{2a}\quad\text{where}$ (67)
$\displaystyle a$
$\displaystyle=\left(\frac{1-\mathcal{C}_{l}}{1-\mathcal{C}_{r}}\right)\left[1-\left(\frac{1-\mathcal{C}_{r}+\mathcal{H}_{r}}{1-\mathcal{C}_{l}+\mathcal{H}_{l}}\right)^{m}\left(\frac{\mathcal{H}_{r}}{1-\mathcal{C}_{r}+\mathcal{H}_{r}}\right)^{n}\right],$
$\displaystyle b$
$\displaystyle=-\left(\frac{1-\mathcal{C}_{l}}{1-\mathcal{C}_{r}}\right)\left[1-\left(\frac{\mathcal{H}_{r}}{1-\mathcal{C}_{r}+\mathcal{H}_{r}}\right)^{n}\right]+\left(\frac{\mathcal{H}_{l}}{1-\mathcal{C}_{l}+\mathcal{H}_{l}}\right)^{n}-1\quad\text{
and }$ $\displaystyle c$
$\displaystyle=1-\left(\frac{1-\mathcal{C}_{l}+\mathcal{H}_{l}}{1-\mathcal{C}_{r}+\mathcal{H}_{r}}\right)^{m}\left(\frac{\mathcal{H}_{l}}{1-\mathcal{C}_{l}+\mathcal{H}_{l}}\right)^{n}.$

Subsequently, substituting Equations (65) and (66) in (64) gives the time-
invariant dimensionless flux in the saturated region $q_{s}$ during
$0<\tau\leq\tau_{p}$ as

$\displaystyle
q_{s}=\frac{\mathfrak{R}-1}{\mathfrak{R}/(1-\mathcal{C}_{l}+\mathcal{H}_{l})^{m}-1/(1-\mathcal{C}_{r}+\mathcal{H}_{r})^{m}}.$
(68)

The time and space invariant flux $q_{s}$ leads to the result using Equations
(61) and (63) that the shock speeds are constant for $0<\tau\leq\tau_{p}$.
Therefore, the shock locations before ponding vary linearly with time. Lastly,
the full solution for this case can be summarized as

u
$\displaystyle=\begin{cases}\textbf{u}_{l},\quad&\zeta<\Lambda_{\mathscr{S}^{*}_{1}}\\\
\textbf{u}_{i_{1}}=\begin{bmatrix}1\\\
1-\mathcal{C}_{l}+\mathcal{H}_{l}\end{bmatrix},\quad&\Lambda_{\mathscr{S}^{*}_{1}}<\zeta/\tau<0\\\
\textbf{u}_{i_{2}}=\begin{bmatrix}1\\\
1-\mathcal{C}_{r}+\mathcal{H}_{r}\end{bmatrix},\quad&0<\zeta/\tau<\Lambda_{\mathscr{S}_{3}}\\\
\textbf{u}_{r},\quad&\zeta/\tau>\Lambda_{\mathscr{S}_{3}}\end{cases}.$ (69)

Here the shock speeds $\Lambda_{\mathscr{S}^{*}_{1}}$ and
$\Lambda_{\mathscr{S}_{3}}$ depend on both left and right states due to
formation of the saturated region as the flux is governed by an elliptic PDE.
This is the case of a temperate firn with more porous and more wet layer lies
on top of less porous and less wet layer (Figure 8d). In this case, the lower
layer is unable to accommodate the flux of water from the upper layer and
therefore leads to a rising perched water table as well as a wetting front
(Figures 8e,f). Figure 8 shows an example of the same where the top layer has
50% porosity and 40% liquid water content ($\phi_{w,l}$) and the bottom layer
is 30% porous and has only 10% LWC. These values correspond to left and right
states being $\textbf{u}_{l}=(0.9,0.4)^{T}$ and $\textbf{u}_{r}=(0.8,0.1)^{T}$
respectively (Figures 8a-c). Consequently, the first and second intermediate
states, $\textbf{u}_{i_{1}}=(1,0.5)^{T}$ and $\textbf{u}_{i_{2}}=(1,0.3)^{T}$
respectively, lie in the expanding saturated region (Figures 8 _d_ -_f_). The
speeds of both fronts $\mathscr{S}_{1}^{*}$ and $\mathscr{S}_{3}$ are constant
before ponding occurs.

Figure 8: Formation of a fully-saturated region in temperate firn (Case VI):
(a) Construction of solution in hodograph plane and their corresponding self-
similar analytical solutions for (b) dimensionless composition and (c)
dimensionless enthalpy with dimensionless velocity $\eta$. The result shown
with dark blue line consists of a backfilling shock, $\mathscr{S}^{*}_{1}$, a
jump, $\mathscr{J}_{2}$, and another wetting shock, $\mathscr{S}_{3}$, along
with two intermediate states
$\textbf{u}_{i_{1}}=(\mathcal{C}_{i_{1}},\mathcal{H}_{i_{1}})^{T}$ and
$\textbf{u}_{i_{2}}=(\mathcal{C}_{i_{2}},\mathcal{H}_{i_{2}})^{T}$. The left
and right states are $\textbf{u}_{l}=(0.9,0.4)^{T}$ and
$\textbf{u}_{r}=(0.8,0.1)^{T}$ respectively. The first and second intermediate
states, $\textbf{u}_{i_{1}}=(1,0.5)^{T}$ and $\textbf{u}_{i_{2}}=(1,0.3)^{T}$
respectively. The evolution of the volume fractions of the three phases in the
resulting system at dimensionless times (d) $\tau=0$, (e) $0.5$, (f) $1.0$.

####  3.3.2 Region 1 only (Ice and gas region), Case VII

In region 1 ($\mathcal{H}\leq 0$), since both fluxes are zero the system is
not strictly hyperbolic and leads to a single wave. As the characteristic
speed, $\lambda_{p}$, is constant, this wave is linearly degenerate and since
$\lambda_{1}=\lambda_{2}=0$, the characteristic is stationary. The resulting
wave is a stationary contact discontinuity $\mathscr{C}$. There won’t be any
transport between the two states since the fluxes of both composition and
enthalpy are zero on either side. In other words, $\textbf{u}_{l}$ will be
connected to $\textbf{u}_{r}$ by a stationary contact discontinuity
$\mathscr{C}$ as

$\displaystyle\textbf{u}_{l}\xrightarrow{\mathscr{C}}\textbf{u}_{r}.$

The resulting solution thus takes the form

u $\displaystyle=\begin{cases}\textbf{u}_{l},\quad&\zeta<0\\\
\textbf{u}_{r},\quad&\zeta>0\end{cases}.$ (70)

Figure 9: Solutions when the left state lies in Region 1 (ice and gas): either
only a contact discontinuity appears (Case VII) or an intermediate state,
$\textbf{u}_{i}$, along with a rarefaction wave $\mathscr{R}_{2}$ also forms
in a $\mathscr{C}_{1}\mathscr{R}_{2}$ fashion (Case VIII). (a) Construction of
solution in hodograph plane and their corresponding self-similar analytical
solutions for (b) dimensionless composition and (c) dimensionless enthalpy
with dimensionless velocity $\eta$. Blue and red lines respectively show the
solutions when the right states are in Regions 1 and 2 respectively. The
evolution of the volume fractions of the three phases in the resulting system
at dimensionless times $\tau=0,0.5,1.0$ for the two configurations: (d-f) Case
VII - $\mathscr{C}$ and (g-i) Case VIII - $\mathscr{C}_{1}\mathscr{R}_{2}$.

This corresponds to the case when a more porous, cold firn lies on top of a
less, cold porous firn (or glacial ice) or vice versa which results in a dry,
static system in this model. An example for a temperate, 40% porous firn on
top of a cold ($T=-19.8^{\circ}$C), 20% porous firn corresponding to
$\textbf{u}_{l}=(0.6,0.0)^{T}$ and $\textbf{u}_{r}=(0.8,-0.1)^{T}$, as shown
in Figures 9a-c (blue lines) and 9d-f. However, in reality the firn may
compact due to overburden cuffey2010physics which is not considered in the
present model. Note that the temperatures (not shown) in the two layers are
the same as the initial condition with sharp transition at $\zeta=0$ due to
the absence of heat diffusion.

###  3.4 Inter-regional cases

####  3.4.1 From Region 1 (ice and gas) to Region 2 (three-phase region), Case
VIII

In this case, the left state lies in region 1 corresponding to the cold firn
and the right state resides in the three-phase region (region 2). The result
is a contact discontinuity $\mathscr{C}_{1}$ onto the solidus
($\mathcal{H}=0$) where the intermediate state lies (see Figure 9a for
example), i.e.,

$\displaystyle\mathcal{H}_{i}=0.$ (71)

Moreover, the intermediate state, $\mathcal{C}_{i}$, lies on the fast path
(constant porosity line) in region 2 satisfying Equation (45) which gives

$\displaystyle\mathcal{C}_{i}$
$\displaystyle=\mathcal{C}_{r}-\mathcal{H}_{r}.$ (72)

Simultaneously, the intermediate state is connected to the right state on the
fast path, resulting in a rarefaction wave $\mathscr{R}_{2}$. It is important
to note that the second wave, $\mathscr{W}_{2}$, is supposed to be faster than
the first and that is why the slow path is avoided in the three-phase region
(region 2). The final solution to this case is the same as given in Equation
(57).

As an example, the left state corresponds to a 30% porous cold layer at
T$=-22.63^{\circ}$C lying on top of the layer corresponding to wet, temperate
firn similar to liquid storage with 80% porosity and 60% LWC ($\phi_{w,r}$),
as shown in Figure 9g. This configuration corresponds to
$\textbf{u}_{l}=(0.7,-0.1)^{T}$ and $\textbf{u}_{r}=(0.8,0.6)^{T}$ with
intermediate state $\textbf{u}_{i}=(0.2,0.0)^{T}$ (Figure 9a, red line). As a
result, the porosity jump remains stationary but the liquid storage drains
downwards due to gravity forming a self-similar rarefaction wave as shown in
Figures 9 _b_ -_c_ (red lines) and 9 _g_ -_i_. In a nutshell, this case
describes the evolution of a more saturated firn layer below a previously
formed less permeable, cold frozen fringe.

####  3.4.2 From Region 2 (three-phase region) to Region 1 (ice and gas)

This case corresponds to porous and wetter firn at melting temperature with
liquid meltwater above the cold firn initially (see Figures 10-12 for
example). In this case, the left state corresponding to the top layer lies in
the three-phase region (region 2) and the right state corresponding to the
bottom layer lies in region 1. Note that the temperature remains subzero only
in the cold region (region 1, $\mathcal{H}<0$) which lies only in the right
state with a sharp transition from the left or intermediate state as heat
diffusion is not considered in this model. There can be four results
corresponding to this initial condition:

Figure 10: Solutions when the left state lies in region 2 (three-phase) and
the right state lies in region 1 without formation of saturated regions:
either only a refreezing front $\mathscr{S}$ (Case IX) appears or a contact
discontinuity, $\mathscr{C}_{1}$, an intermediate state, $\textbf{u}_{i}$,
along with a refreezing front $\mathscr{S}_{2}$ forms in a
$\mathscr{C}_{1}\mathscr{S}_{2}$ (Case X) fashion. (a) Construction of
solution in hodograph plane and their corresponding self-similar analytical
solutions for (b) dimensionless composition and (c) dimensionless enthalpy
with dimensionless velocity $\eta$. Red and blue lines respectively show the
solutions when the right states are in region 2 and region 1 respectively. The
evolution of the volume fractions of the three phases in the resulting system
at dimensionless times $\tau=0,0.5,1.0$ for the two configurations: (d-f) Case
VII - $\mathscr{C}$ and (g-i) Case VIII - $\mathscr{C}_{1}\mathscr{R}_{2}$.

###### (i.) Shock (Case IX):

When the right state in region 1 lies on the fast path in region 2 (three-
phase region) extended to region 1 (ice and gas region) referred to as
extended fast path, it results in only a single moving shock $\mathscr{S}$
(see Figure 10a, blue line). Note that the extended fast paths are not the
constant porosity contours in Region 1 (Figures 1 _c_ ,_g_). Mathematically,
the left and right states are connected by the relation

$\displaystyle\mathcal{C}_{r}$
$\displaystyle=\mathcal{H}_{r}+\mathcal{C}_{l}-\mathcal{H}_{l}.$ (73)

The shock speed $\Lambda_{\mathscr{S}}$ can then be evaluated using Rankine-
Hugoniot condition (52). The final solution of this case is provided in
Equation (53). This case corresponds to the part where meltwater percolates
and precipitates directly into ice, reducing the porosity equal to the left
state thereby extending the frozen fringe (Figures 10 _d_ -_f_). Physically,
the water seeps from the top layer to the bottom layer due to gravity and a
part of it precipitates due to heat loss to the surrounding cold firn. Since
the flux in the right state is zero, the evaluation depends on the difference
of either dimensionless composition or enthalpy. Figures 10 _a_ -_c_ (blue
line) shows an example of the solution where temperate firn with 45% porosity
($\varphi_{l}$) and 25% LWC ($\phi_{w,l}$) initially lies on a cold layer of
porosity 50% at T$=-15.84^{\circ}$C (Figure 10d). These initial conditions
correspond to a left state, $\textbf{u}_{l}=(0.8,0.25)^{T}$ and a right state
$\textbf{u}_{r}=(0.5,-0.05)^{T}$ (Figures 10a-c, blue line). A refreezing
shock $\mathscr{S}$ results that moves downwards while warming up the
surrounding snow by partially refreezing (Figures 10 _d_ -_f_). This reduces
the porosity behind the wetting front from $50\%$ to $45\%$, equal to the same
value as the left state thereby extending the previously formed frozen fringe
in the top layer into the bottom layer.

###### (ii.) 1-Contact discontinuity and 2-Shock (Case X):

This case occurs when the left state cannot be directly connected to the right
state along the extended fast path as discussed in Case IX. So in this case,
the left state $\textbf{u}_{l}$ in region 2 is more porous and connected to an
intermediate state, $\textbf{u}_{i}$, along the slow path where the
dimensionless composition of the intermediate state, $\mathcal{C}_{i}$, is
less than unity, to keep the medium unsaturated (see Figure 10a, red line).
Similar to the only shock $\mathscr{S}$ case (Case IX), the intermediate state
lies on the extended fast path from region 2 to region 1. Therefore the
solution is a combination of a stationary contact discontinuity
$\mathscr{C}_{1}$ and a moving shock $\mathscr{S}_{2}$. The first
characteristic wave, contact discontinuity $\mathscr{C}_{1}$, satisfies
Equation (44) for left state $\textbf{u}_{l}$ and intermediate state
$\textbf{u}_{i}$. Since the second characteristic lies on the Hugoniot locus,
a shockwave $\mathscr{S}_{2}$ results, that satisfies the Hugoniot-jump
condition (39) for the intermediate state $\textbf{u}_{i}$ and right state
$\textbf{u}_{r}$. Combining all these equations gives a simple relation for
dimensionless enthalpy and composition at the intermediate state as

$\displaystyle\mathcal{H}_{i}=\mathcal{H}_{l}\left(\frac{1-\mathcal{C}_{l}+\mathcal{H}_{l}}{1-\mathcal{C}_{r}+\mathcal{H}_{r}}\right)^{\frac{m-n}{n}}\quad\text{and}\quad\mathcal{C}_{i}=\mathcal{H}_{i}+\mathcal{C}_{r}-\mathcal{H}_{r}.$
(74)

The dimensionless velocity of the shock wave $\mathscr{S}_{2}$ (refreezing
front) is

$\displaystyle\Lambda_{\mathscr{S}_{2}}$
$\displaystyle=\frac{f(\textbf{u}_{i})-f(\textbf{u}_{r})}{\mathcal{H}_{i}-\mathcal{H}_{r}}=\frac{\mathcal{H}_{i}^{n}(1-\mathcal{C}_{i}+\mathcal{H}_{i})^{m-n}}{\mathcal{H}_{i}-\mathcal{H}_{r}}\quad\textrm{(since
$f(\textbf{u}_{r})=0$)}.$ (75)

The full solution is given in Equation (59) with shock speed (75). Since
$\mathcal{H}_{r}<0$, $\mathcal{H}_{i}>0$ and $f(\textbf{u}_{r})=0$, the speed
of the refreezing front $\mathscr{S}_{2}$ is slower than temperate firn case
(Case V) due to refreezing. Figures 10 _a_ -_c_ (red line) shows an example of
such a solution for light rainfall on a multilayered firn with coarse to fine
transition. The left state corresponds to the wetter, temperate firn on top
with a porosity of $70\%$ and a liquid water content $\phi_{w,l}=0.1$. The
right state corresponds to a cold and dry firn layer with 35% porosity and a
temperature of T$=-19.49^{\circ}$C. These values for left and right states
correspond to states being $\textbf{u}_{l}=(0.4,0.1)^{T}$ and
$\textbf{u}_{r}=(0.65,-0.08)^{T}$, respectively. Here, a stationary contact
discontinuity $\mathscr{C}_{1}$ stays at the surface leading to a growing
intermediate state $\textbf{u}_{i}=(0.893,0.163)$ formed due to partial
refreezing of the rainwater which warms the surrounding firn (Figures 10 _g_
-_i_). The intermediate state with $27\%$ porosity and $16.3\%$ LWC expands
with time as the refreezing front $\mathscr{S}_{2}$ infiltrates further into
the right state. Note that the porosity in the intermediate state is smaller
than the right state which decreases further due to meltwater refreezing
leading to the formation of a fresh frozen fringe. This phenomenon has been
observed in the model by meyer2017continuum who used the dimensional form of
Equation (75) with distinct densities of ice and water phases to explain the
field data from humphrey2012thermal.

###### (iii.) 1-Backfilling shock, 2-Jump and 3-Shock (formation of a
saturated region, Case XI):

Figure 11: Formation of a fully-saturated region when right state lies in
region 1: (a) Construction of solution in hodograph plane and their
corresponding self-similar analytical solutions for (b) dimensionless
composition and (c) dimensionless enthalpy with dimensionless velocity $\eta$.
The result shown with dark blue line consists of a backfilling shock,
$\mathscr{S}^{*}_{1}$, a jump, $\mathscr{J}_{2}$, and another refreezing
shock, $\mathscr{S}_{3}$, along with two intermediate states
$\textbf{u}_{i_{1}}=(\mathcal{C}_{i_{1}},\mathcal{H}_{i_{1}})^{T}$ and
$\textbf{u}_{i_{2}}=(\mathcal{C}_{i_{2}},\mathcal{H}_{i_{2}})^{T}$. The
evolution of the volume fractions of the three phases in the resulting system
at dimensionless times (d) $\tau=0$, (e) $0.1$, (f) $0.2$.

Here the slow path originating from the left state and the extended fast path
emanating from the right state do not intersect where $\mathcal{C}<1$, making
it a different case than Case X. If the intermediate state lies on the
saturated region, then the hyperbolic nature of the solution breaks down,
similar to Case VI for temperate firn (see Table 2). So the state solution
(69) is the same as provided for the temperate firn (Case VI) consisting of a
backfilling shock $\mathscr{S}^{*}_{1}$ moving upwards with speed (61), a
stationary jump $\mathscr{J}_{2}$ at the location of initial jump $\zeta=0$
and a “refreezing” front $\mathscr{S}_{3}$ moving downwards with speed
provided in Equation (63). The flux in the saturated region $q_{s}$ is again
provided by Equation (64). Note that the flux at the right state
$f(\textbf{u}_{r})=0$ since it lies in region 1. Invoking the ansatz for a
constant shock speed ratio, similar to Case VI, provides an analytic relation
for a constant flux in the saturated region, $q_{s}$, same as Equation (68)
with both shocks moving in opposite directions at constant speeds. However
since the flux at the right state is zero, the relation for the ratio of shock
speeds, $\mathfrak{R}$, earlier provided by Equation (67), now simplifies to

$\displaystyle\mathfrak{R}(\textbf{u}_{l},\textbf{u}_{r})$
$\displaystyle=\frac{-b-\sqrt{b^{2}-4ac}}{2a}\quad\text{where}$ (76)
$\displaystyle a$
$\displaystyle=\left(\frac{1-\mathcal{C}_{l}}{1-\mathcal{C}_{r}}\right),$
$\displaystyle b$
$\displaystyle=-\left(\frac{1-\mathcal{C}_{l}}{1-\mathcal{C}_{r}}\right)+\left(\frac{\mathcal{H}_{l}}{1-\mathcal{C}_{l}+\mathcal{H}_{l}}\right)^{n}-1\quad\text{
and }$ $\displaystyle c$
$\displaystyle=1-\left(\frac{1-\mathcal{C}_{l}+\mathcal{H}_{l}}{1-\mathcal{C}_{r}+\mathcal{H}_{r}}\right)^{m}\left(\frac{\mathcal{H}_{l}}{1-\mathcal{C}_{l}+\mathcal{H}_{l}}\right)^{n}.$

Figure 11 shows an example of this case where the left state corresponds to
wet and temperate layer with 80% porosity ($\varphi_{l}$) and 65% LWC
($\phi_{w,l}$) lying on top of very cold and less porous (T$=-30^{\circ}$C,
$\varphi_{r}=50\%$) firn. This corresponds to the left and right states being
$\textbf{u}_{l}=(0.85,0.65)^{T}$ and $\textbf{u}_{r}=(0.5,-0.095)^{T}$. The
first and second intermediate states are $\textbf{u}_{i_{1}}=(1,0.8)^{T}$ and
$\textbf{u}_{i_{2}}=(1,0.405)^{T}$ respectively. The upper shock, also called
rising perched water table, is almost four times faster than the lower,
refreezing front (see Figures 11 _b_ -_f_). Below the initial jump at
$\zeta=0$ the porosity is initially reduced which is further reduced by
refreezing, leading to formation of a frozen fringe (Figures 11 _e_ -_f_). The
second intermediate state below jump $\mathscr{J}_{2}$ with a newly formed
frozen fringe is unable to accommodate the whole volumetric flux of water,
leading to the formation of a rising perched water table (see Figures 11 _e_
-_f_). Once the rising perched water table reaches the surface, it will lead
to ponding and can eventually form runoff.

###### (iv.) 1-Backfilling shock, 2-Jump, 3-Contact discontinuity (formation
of impermeable ice lens, Case XII):

This final case captures the formation of impermeable ice lens through
convection of meltwater and associated latent heat. Ice lens forms through
Case XI, if the extended fast path ($\mathcal{H}=\mathcal{C}+\mathfrak{C}$)
emanating from the right state reach $\mathcal{C}\geq 1$ at the solidus
($\mathcal{H}=0$) or the porosity of the second intermediate state reaches 0.
This case is entirely governed by the right state and ice lens will form if
and only if the right state satisfies

$1-\mathcal{C}_{r}+\mathcal{H}_{r}\leq 0.$ (77)

The exact location of the right state does not matter if it lies in the region
of ice lens formation (grey region in Figure 12 _a_). In other words, an
impermeable ice lens via heat convection will form if and only if the cold
content of the firn exceeds the latent enthalpy of incoming meltwater. The
mathematical condition (77) in dimensional form has been given in
humphrey2021physical for distinct densities of ice and water. Furthermore,
this condition has been qualitatively discussed in colbeck1976analysis;
clark2017analytical. In this limit, the solution for Case XI breaks down due
to the formation of an impermeable ice lens. Figure 12 shows the region of ice
lens formation which is perpetuated by very low firn porosity or temperatures
in the right state. The final solution for this case is a backfilling shock,
$\mathscr{S}^{*}_{1}$ to the intermediate saturated region $i_{1}$. The first
saturated region is connected to the second intermediate state $i_{2}$ through
the jump $\mathscr{J}_{2}$. The second intermediate state $i_{2}$ is the
infinitesimally thin ice lens (see Figures 12 _e_ -_f_) which blocks the
further meltwater. The solution in this case can be written as

$\displaystyle\textbf{u}_{l}\xrightarrow{\mathscr{S}^{*}_{1}}\textbf{u}_{i_{1}}\mathrel{\mathchoice{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar
0\relax\symAMSa
04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar
0\relax\symAMSa
04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptscriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar
0\relax\symAMSa
04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar
0\relax\symAMSa
04B}{}}}\textbf{u}_{i_{2}}\xrightarrow{\mathscr{C}_{3}}{\textbf{u}_{r}}$

where the state solution is quantified as

u
$\displaystyle=\begin{cases}\textbf{u}_{l},\quad&\zeta<\Lambda_{\mathscr{S}^{*}_{1}}\\\
\textbf{u}_{i_{1}}=\begin{bmatrix}1\\\
1-\mathcal{C}_{l}+\mathcal{H}_{l}\end{bmatrix},\quad&\Lambda_{\mathscr{S}^{*}_{1}}<\zeta/\tau<0\\\
\textbf{u}_{i_{2}}=\begin{bmatrix}1\\\
0\end{bmatrix},\quad&0<\zeta<\textrm{d}\zeta\\\
\textbf{u}_{r},\quad&\zeta>\textrm{d}\zeta\end{cases},$ (78)

where $\textrm{d}\zeta$ is the infinitesimally small thickness of the ice
lens.

In this case, only a single shock moves upwards similar to the filling of a
bucket. The flux in the saturated region is simply zero from the mass balance
at the ice lens, i.e., $q_{s}=f(\textbf{u}_{i_{1}})=f(\textbf{u}_{i_{2}})=0$.
The speed of the rising perched water table (first shock
$\mathscr{S}^{*}_{1}$) is then

$\displaystyle\Lambda_{\mathscr{S}^{*}_{1}}(\textbf{u}_{l})=\frac{\textrm{d}\zeta_{U}}{\textrm{d}\tau}=\frac{f(\textbf{u}_{l})}{\mathcal{C}_{l}-1}=\frac{\mathcal{H}^{n}(1-\mathcal{C}_{l}+\mathcal{H}_{l})^{(m-n)}}{\mathcal{C}_{l}-1}$
(79)

which depends only on the left-state variables.

Figure 12: Impermeable ice lens formation: (a) Construction of solution in
hodograph plane and their corresponding self-similar analytical solutions for
(b) dimensionless composition and (c) dimensionless enthalpy with
dimensionless velocity $\eta$. The result shown with dark blue line consists
of a backfilling shock, $\mathscr{S}^{*}_{1}$, a jump, $\mathscr{J}_{2}$ and a
contact discontinuity, $\mathscr{C}_{3}$ along with two intermediate states
$\textbf{u}_{i_{1}}=(\mathcal{C}_{i_{1}},\mathcal{H}_{i_{1}})^{T}$ and
$\textbf{u}_{i_{2}}=(\mathcal{C}_{i_{2}},\mathcal{H}_{i_{2}})^{T}$. The second
intermediate state $\textbf{u}_{i_{2}}$ corresponds to the impermeable ice
lens. The grey region corresponds to $1-\mathcal{C}+\mathcal{H}\leq 0$ where
the right state resides to cause impermeable ice lens formation. The evolution
of the volume fractions of the three phases in the resulting system at
dimensionless times (d) $\tau=0$, (e) $0.02$, (f) $0.04$.

It can be deduced that either a cold firn (less $\mathcal{H}$) or a nearly
non-porous firn with $\mathcal{C}$ close to unity will induce the formation of
an impermeable ice lens via refreezing caused by heat advection. For example,
Figure 12 shows such a solution when wet and temperate layer of 70% porosity
and 65% LWC on top of a cold layer with 5% porosity at T$=-30^{\circ}$C. These
conditions correspond to $\textbf{u}_{l}=(0.95,0.65)^{T}$ and
$\textbf{u}_{r}=(0.95,-0.180)^{T}$ (Figures 12a-c). Here the intermediate
states lie at $\textbf{u}_{i_{1}}=(1,0.7)^{T}$ and
$\textbf{u}_{i_{2}}=(1,0)^{T}$ respectively. An impermeable ice lens forms as
a result (see Figures 12 _e_ -_f_) that blocks the flow of meltwater
downwards. Since there is no meltwater percolating in the lower layer, the
backfilling occurs very rapidly.

Table 2: Summary of all the analytical solutions presented in this paper along with related works which have either studied or observed the corresponding scenario. Case | Region |  State properties | Solution |  Physical relevance |  Related work(s)
---|---|---|---|---|---
$\textbf{u}_{l}$ | $\textbf{u}_{r}$ | Summary |  Equation(s) |  Construction figure
I | 2 | 2 |  $f(\textbf{u}_{l})=f(\textbf{u}_{r})$$\Leftrightarrow$$\textbf{u}_{l}$, $\textbf{u}_{r}$ connected by slow path | $\textbf{u}_{l}\xrightarrow{\mathscr{C}}\textbf{u}_{r}$ |  44, 46, 47 |  6 (green line) |  Steady meltwater flux inside a temperate firn with a jump in porosity. |  $-$
II | 2 | 2 |  $f(\textbf{u}_{l})<f(\textbf{u}_{r})$; $\textbf{u}_{l}$, $\textbf{u}_{r}$ connected by fast path; $\mathcal{C}_{l}<\mathcal{C}_{r}$ | $\textbf{u}_{l}\xrightarrow{\mathscr{R}}\textbf{u}_{r}$ |  45, 51 |  6 (blue line) |  A sudden decrease in meltwater flux inside a temperate firn with constant porosity. |  colbeck1976analysis; singh1997kinematic; clark2017analytical
III | 2 | 2 |  $f(\textbf{u}_{l})>f(\textbf{u}_{r})$; $\textbf{u}_{l}$, $\textbf{u}_{r}$ connected by fast path; $\mathcal{C}_{l}>\mathcal{C}_{r}$ | $\textbf{u}_{l}\xrightarrow{\mathscr{S}}\textbf{u}_{r}$ |  45, 52, 53 |  6 (red line) |  A sudden increase in meltwater flux inside a temperate firn with constant porosity. |  colbeck1971one; colbeck1972theory; singh1997kinematic; samimi2021time
IV | 2 | 2 |  $f(\textbf{u}_{l})<f(\textbf{u}_{r})$; $\textbf{u}_{l}$, $\textbf{u}_{r}$ cannot be connected by either paths; $\mathcal{C}_{i}<1$ | $\textbf{u}_{l}\xrightarrow{\mathscr{C}_{1}}\textbf{u}_{i}\xrightarrow{\mathscr{R}_{2}}\textbf{u}_{r}$ |  55, 57 |  7 (blue line) |  A sudden decrease in meltwater flux inside a temperate firn with a step change in porosity. |  $-$
V | 2 | 2 |  $f(\textbf{u}_{l})>f(\textbf{u}_{r})$; $\textbf{u}_{l}$, $\textbf{u}_{r}$ cannot be connected by either paths; $\mathcal{C}_{i}<1$ | $\textbf{u}_{l}\xrightarrow{\mathscr{C}_{1}}\textbf{u}_{i}\xrightarrow{\mathscr{S}_{2}}\textbf{u}_{r}$ |  58, 59, 52 |  7 (red line) |  A sudden increase in meltwater flux inside a temperate firn with a step change in porosity. |  $-$
VI | 2 | 2 |  $f(\textbf{u}_{l})>f(\textbf{u}_{r})$; $\textbf{u}_{l}$, $\textbf{u}_{r}$ cannot be connected by the combination of two paths in Region 2; slow path emanating from $\textbf{u}_{l}$ does not intersect with fast path emanating from $\textbf{u}_{r}$ where $\mathcal{C}<1$ | $\textbf{u}_{l}\xrightarrow{\mathscr{S}^{*}_{1}}\textbf{u}_{i_{1}}\mathrel{\mathchoice{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptscriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}}\textbf{u}_{i_{2}}\xrightarrow{\mathscr{S}_{3}}\textbf{u}_{r}$ |  61, 63, 64, 65, 67, 68, 69 |  8 |  A sudden increase in meltwater flux inside a temperate firn with a step decline in porosity with depth leading to formation of a rising perched water table. |  shadab2022analysis
VII | 1 | 1 |  $f(\textbf{u}_{l})=f(\textbf{u}_{r})=0$ | $\textbf{u}_{l}\xrightarrow{\mathscr{C}}\textbf{u}_{r}$ |  70 |  9 (blue line) |  A cold firn with a step change in porosity with no meltwater flux. |  $-$
VIII | 1 | 2 |  $-$ | $\textbf{u}_{l}\xrightarrow{\mathscr{C}_{1}}\textbf{u}_{i}\xrightarrow{\mathscr{R}_{2}}\textbf{u}_{r}$ |  57, 71, 72 |  9 (red line) |  A cold firn with no meltwater flux overlying a wet, temperate firn. |  $-$
IX | 2 | 1 |  $\textbf{u}_{l}$, $\textbf{u}_{r}$ connected by extended fast path | $\textbf{u}_{l}\xrightarrow{\mathscr{S}}\textbf{u}_{r}$ |  52, 53, 73 |  10 (blue line) |  A sudden increase in meltwater flux into a cold firn with constant porosity leading to the formation of frozen fringe. |  colbeck1976analysis; clark2017analytical; meyer2017continuum
X | 2 | 1 |  $\textbf{u}_{l}$, $\textbf{u}_{r}$ cannot be connected by extended fast path, $\mathcal{C}_{i}<1$ | $\textbf{u}_{l}\xrightarrow{\mathscr{C}_{1}}\textbf{u}_{i}\xrightarrow{\mathscr{S}_{2}}\textbf{u}_{r}$ |  59, 74, 75 |  10 (red line) |  A sudden increase in meltwater flux into a cold firn with a step change porosity leading to the formation of frozen fringe. |  colbeck1976analysis; clark2017analytical; meyer2017continuum
XI | 2 | 1 |  $\textbf{u}_{l}$, $\textbf{u}_{r}$ cannot be connected by extended fast path; slow path from $\textbf{u}_{l}$ does not intersect with extended fast path from $\textbf{u}_{r}$ where $\mathcal{C}<1$ | $\textbf{u}_{l}\xrightarrow{\mathscr{S}^{*}_{1}}\textbf{u}_{i_{1}}\mathrel{\mathchoice{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptscriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}}\textbf{u}_{i_{2}}\xrightarrow{\mathscr{S}_{3}}\textbf{u}_{r}$ |  61, 63, 64, 65, 66, 76, 68, 69 |  11 |  A sudden increase in meltwater flux into a cold firn with a step decrease porosity leading to the formation of frozen fringe and a rising perched water table. |  humphrey2021physical
XII | 2 | 1 |  $\textbf{u}_{r}$ lies in impermeable ice lens region satisfying $1-\mathcal{C}_{r}-\mathcal{H}_{r}\leq 0$ | $\textbf{u}_{l}\xrightarrow{\mathscr{S}^{*}_{1}}\textbf{u}_{i_{1}}\mathrel{\mathchoice{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptscriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}{\mathrel{}\mathrel{\mathop{\dabar@\dabar@}\limits^{\hbox{\set@color$\scriptstyle\,\mathscr{J}_{2}$}}}\mathrel{\mathchar 0\relax\symAMSa 04B}{}}}\textbf{u}_{i_{2}}\xrightarrow{\mathscr{C}_{3}}{\textbf{u}_{r}}$ |  78, 79 |  12 |  A sudden increase in meltwater flux into a cold firn with a step decrease porosity leading to the formation an impermeable ice lens and a rising perched water table. |  humphrey2021physical

##  4 Infiltration into a multilayered firn - Validation with numerical
simulations

Here we demonstrate the application of this theory to a realistic multilayered
firn leading to the formation of a perched firn aquifer. This final test shows
the infiltration into multilayered firn after a melt event combining two cases
proposed in Section 3 that are summarized in Table 2 (see Figure 13 _a_). This
problem summates the commonly studied wetting front propagation in a temperate
region (Case III) with the wetting front in a cold region and the formation of
a perched water table (Case XI) that is not studied well in the literature.
This problem demonstrates a delay in meltwater ponding at the surface due to a
decay in both firn porosity and temperature with depth. The firn is 70%
porous, dry and at 0∘C, above a dimensionless depth of $\zeta=1$ (Figures 13
_a_). At time $\tau=0$, the meltwater is generated at the surface ($\zeta=0$),
which keeps the liquid water content (denoted by $LWC$ or $\phi_{w}$) to 40%
making the firn very wet at the surface, $\zeta=0$. The analytic solution of
this problem can be constructed in the hodograph plane in two parts (Figure 13
_b_).

First, a wetting front $\mathscr{S}$ initially propagates downwards with a
constant dimensionless speed $\Lambda_{\mathscr{S}}=0.28$, as given in Case
III (Figures 13 _c_ -_e_ , red dashed line). The dimensionless time when the
initial front reaches the depth of transition $\zeta_{t}=1$ to form a
saturated region is
$\tau_{s}=\frac{\zeta_{t}}{\Lambda_{\mathscr{S}}}=1/0.28=3.57$, as shown in
Case XI. Afterwards the saturated region forms due to large meltwater flux
compared to the hydraulic conductivity of the second intermediate region
formed due to refreezing below the jump $\zeta>1$, making the medium pressure
driven instead and enforcing $q_{s}(\tau)\leq K_{i_{2}}$. This saturated
region expands in both directions as a perched water table (upper shock shown
by green dashed line in Figure 13) rising to the surface and the wetting front
${\mathscr{S}_{3}}$ that percolates in the cold region (lower shock shown with
blue dashed line) while refreezing a part of meltwater to warm the surrounding
snow. Thus, it shows a reduction in porosity from 30% to 21.2% for $\zeta>1$
behind the wetting front $\mathscr{S}_{3}$ at $\tau>\tau_{s}$ (Figure 13 _c_).
The perched water table ${\mathscr{S}^{*}_{1}}$ rises upwards with constant
dimensionless velocity $\Lambda_{\mathscr{S}^{*}_{1}}=-0.268$ until ponding
occurs. Meanwhile, the wetting front keeps percolating in the cold firn with
reduced velocity $\Lambda_{\mathscr{S}_{3}}=0.105$. Lastly, the ponding starts
at a time when the rising perched water table reaches at the surface so the
dimensionless ponding time can be calculated theoretically as
$\tau_{p}=\zeta_{t}/\Lambda_{\mathscr{S}}+(-\zeta_{t})/\Lambda_{\mathscr{S}^{*}_{1}}=1/0.28+(-1)/(-0.268)=7.30$.
All of these dimensionless shock speeds and times are computed analytically
and the resulting locations are graphed with dashed lines in Figures 13 _c_
-_e_. The theoretical results show an excellent comparison with the numerical
solutions shown in contour plots. The numerical solutions are obtained by
solving the governing model (32) performed in the absence of capillary effects
in between water and gas phases as well as heat diffusion along with the
treatment in the saturated region given by shadab2022hyperbolic. Further, the
densities of the water and the ice phases are assumed to be the same. The
computational domain $\zeta\in[0,2]$ is divided uniformly into 400 cells. The
boundary condition at the top surface ($\zeta=0$) is prescribed as the “Top
condition” in Figure 13 _a_ whereas the bottom boundary condition is not
required. More details about the numerical model are given in
shadab2023mechanism. This problem constitutes a very specific benchmark test
for firn hydrology simulators that are able to simulate variably saturated
flows. It shows how vertical heterogeneity in firn may lead to formation of
perched aquifers that can cause ponding at a later stage. It also illustrates
that an impermeable layer is not required to cause meltwater perching and
subsequently, ponding.

Figure 13: Infiltration into a multilayered firn with porosity and temperature
decay with depth: (a) Schematic diagram showing all of the layers (b)
Construction of solution in hodograph plane. The contours showing evolution of
the firn (c) porosity $\varphi$, (d) liquid water content $LWC$ or volume
fraction of water $\phi_{w}$ and temperature $T$ evaluated by the numerical
simulator. Here all dashed lines show analytic solutions computed from the
proposed theory. The thin, grey dashed lines show theoretically calculated
dimensionless times of saturation $\tau_{s}$ and ponding $\tau_{p}$. The
theoretical evolution of the initial wetting front $\mathscr{S}$ (red dashed
line) is computed from Case III whereas the dynamics of saturated region after
wetting front $\mathscr{S}$ reaches $\zeta=1$ shown by blue and green dashed
lines is computed by Case XI. Here $\delta$ and $t_{c}=\delta/K_{h}$ are
characteristic times with former being calculated from their definition 21.
For example, the characteristic depth is $\delta=5$ m and for
$K_{h}=5.45\times 10^{-4}$ m/s (see Table 1), the characteristic time comes
out to be $t_{c}=2.53$ hours. The dimensionless wetting front speeds can be
redimensionalized by multiplying with saturated (and no-matrix) hydraulic
conductivity $K_{h}$.

##  5 Conclusions

This work introduces a unified kinematic wave theory for meltwater
infiltration into firn that helps construct analytic solutions in a hodograph
plane. The theory neglects heat diffusion and capillary forces while assuming
a constant density for water and ice phases. We provide a suite of 12 cases of
melt infiltration into firn (except Case VII), inspired by nature, and
construct their analytic solution while connecting most of the cases given in
the literature. The previous works were predominantly limited to unsaturated
wetting fronts in cold and temperate firn. The proposed theory considers the
cases that have not been studied previously such as the formation of a perched
water table (Case VI and Case XI) and the region of formation of an
impermeable ice lens (Case XII). These simple cases can be combined to study a
more realistic problem, as was demonstrated by the problem of infiltration
into multilayered firn. There are several consequences of this work. First,
one can interpret the physics of the meltwater infiltration into firn without
running expensive numerical simulations. It will better constrain the process
of firn densification and subsequent parameterization for transforming the
altimetry data into gravimetry data. Second, these analytic solutions can thus
help develop better, cost-effective physics-based firn hydrological models
which can then be integrated with ice-sheet and Earth system models. Further,
these problems can serve as a benchmark for next generation of wet firn
hydrological models which currently show significant deviation due to a lack
of benchmark problems. Lastly, from advanced remote sensing and field
observation techniques, the proposed cases will help deduce the firn
properties such as heterogeneity and percolation parameters such as (relative)
permeability exponent. This comprehensive framework can significantly enhance
our understanding of wet firn hydrology, a component that has been poorly
understood, ultimately aiding in constraining its contribution to surface mass
balance loss from glaciers and, consequently, sea-level rise.

##  6 Acknowledgements and Funding

Support for this work was provided through NASA under Emerging World Grant
number NASA $18-$EW$18\\_2-0027$ and the University of Texas Institute for
Geophysics under the Blue Sky Student Fellowship. The authors would like to
acknowledge Dr. Surendra Adhikari and Dr. Andreas Colliander for the
discussions on meltwater infiltration at the Dye-2 site in Southwest
Greenland.

##  Appendix A Flux Gradient and Lemma related to eigendecomposition

The partial derivatives of the flux function given in equation (30) give:

$\displaystyle f_{,\mathcal{H}}$ $\displaystyle=\frac{\partial
f_{\mathcal{C}}}{\partial\mathcal{H}}=\begin{cases}0,&\mathcal{H}\leq 0\\\
\mathcal{H}^{n-1}(1-\mathcal{C}+\mathcal{H})^{m-n-1}(n(1-\mathcal{C})+m\mathcal{H}),&0<\mathcal{H}<\mathcal{C}\end{cases}$
(80) $\displaystyle f_{,\mathcal{C}}$ $\displaystyle=\frac{\partial
f_{\mathcal{H}}}{\partial\mathcal{C}}=\begin{cases}0,&\mathcal{H}\leq 0\\\
-(m-n)\mathcal{H}^{n}(1-\mathcal{C}+\mathcal{H})^{m-n-1},&0<\mathcal{H}<\mathcal{C}\end{cases}$
(81)

######  Lemma A.1.

Prove that $\lambda_{2}$ is non-negative and increases monotonically in the
direction of integral curves corresponding to the second eigenvector
$\textbf{r}_{2}$ (fast paths) in three phase region,
$0<\mathcal{H}<\mathcal{C}$.

###### Proof.

Since $\mathcal{H}>0$, $-\mathcal{C}<\mathcal{H}-\mathcal{C}$, and
$0\leq\mathcal{C}<1$,

$\displaystyle\lambda_{2}$
$\displaystyle=n\mathcal{H}^{n-1}(1-\mathcal{C}+\mathcal{H})^{m-n}>0$ (82)

Along the fast path, $\mathcal{C}=\mathcal{H}+\mathfrak{C}$,

$\displaystyle\lambda_{2}$
$\displaystyle=n\mathcal{H}^{n-1}(1-\mathcal{C}+\mathcal{H})^{m-n}$ (83)
$\displaystyle=n\mathcal{H}^{n-1}(1-(\mathcal{H}+\mathfrak{C})+\mathcal{H})^{m-n}$
(84) $\displaystyle=n\mathcal{H}^{n-1}(1+\mathfrak{C})^{m-n}$ (85)

Now, taking the derivative with respect to $\mathcal{H}$, we get

$\displaystyle\frac{\textrm{d}{\lambda}_{2}}{\textrm{d}\mathcal{H}}$
$\displaystyle=n(n-1)\mathcal{H}^{n-2}(1+\mathfrak{C})^{m-n}$ (86)

As $\mathcal{H}>0$ and $\mathfrak{C}=\mathcal{C}-\mathcal{H}<1$ in three phase
region,

$\displaystyle\frac{\textrm{d}{\lambda}_{2}}{\textrm{d}\mathcal{H}}>0,\quad\textrm{
for }n>1$ (87)

Therefore, $\lambda_{2}$ is non-negative and increases monotonically in the
direction of $\textbf{r}_{2}$ when $n>1$. ∎

######  Lemma A.2.

Prove that fast paths are parallel to constant porosity $\varphi$ contours in
three phase region, $0<\mathcal{H}\leq\mathcal{C}$.

###### Proof.

In the three-phase region, the porosity
$\varphi=\phi_{w}+\phi_{g}=\mathcal{H}+(1-\mathcal{C})$. For a constant
porosity $\textrm{d}\varphi=0$, therefore

$\displaystyle\textrm{d}\varphi$
$\displaystyle=\textrm{d}(\mathcal{H}+(1-\mathcal{C}))$ (88) $\displaystyle 0$
$\displaystyle=\textrm{d}\mathcal{H}-\textrm{d}\mathcal{C}$ (89)
$\displaystyle 1$
$\displaystyle=\frac{\textrm{d}\mathcal{H}}{\textrm{d}\mathcal{C}}$ (90)
$\displaystyle\mathcal{C}$ $\displaystyle=\mathcal{H}+\mathfrak{C}$ (91)

where $\mathfrak{C}$ again is an integration constant. Therefore, the constant
porosity $\varphi$ contours are the integral curves corresponding to second
eigenvector $\textbf{r}_{2}$ (fast paths) in three-phase region. ∎