2.1. Conserved quantities

We model firn as a three-phase system comprising liquid water (w), ice (i) and gas (g). These three phases are composed of two components, namely water (H₂O) and air (∼Nitrogen gas, N₂). The water component partitions into the liquid and ice, while the air component is confined to the gas phase. The first conserved variable is the water ( $\mathrm{H}_2\mathrm{O}$) composition (kg/m³), C, defined as the total mass of water component per unit representative elemental volume (REV) as

(1)\begin{align} {C = \rho_i \phi_i + \rho_w \phi_w,} \end{align}

where ρ_α is the density (kg/m³) and ϕ_α refers to volume fraction of the phase $\alpha \in \{w,i,g \}$. The formulation assumes that the water vapor component is negligible in the gas phase and that ice and water phases are pure. Assuming the same density for ice and water, $\rho_i\approx \rho_w = \rho$, and the volume fraction constraint, $\phi_i+\phi_w+\phi_g=1$, Eqn (1) further simplifies to

(2)\begin{align} C &= \rho(\phi_i +\phi_w )=\rho (1-\phi_g). \end{align}

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 (Jordan, Reference Jordan1991; Alexiades and Solomon, Reference Alexiades and Solomon1993; Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012; Carnahan and others, Reference Carnahan, Wolfenbarger, Jordan and Hesse2021). Therefore, the second conserved variable is chosen to be the enthalpy of the system (J/m³), H, which is defined as

(3)\begin{align} H &:= \rho_i \phi_i h_i(T) + \rho_w \phi_w h_w(T) + \rho_g \phi_g h_g(T), \end{align}

where h_α is the specific enthalpy (J/kg) of the phase α, 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_α of each phase $\alpha \in \{w,i,g\}$ can then be defined as

(4)\begin{align} h_i(T) &= \begin{cases}c_{p,i}(T-T_m),& T \lt T_m \quad (\textrm{or } H \lt 0) \\0, & T \geqslant T_m \quad(\textrm{or } H \geqslant 0) \end{cases}, \end{align}

(5)\begin{align} h_w(T) &= \begin{cases}0,& T \lt T_m \quad(\textrm{or } H \lt 0) \\c_{p,w}(T-T_m)+L, & T \geqslant T_m \quad(\textrm{or } H\geqslant 0)\end{cases}, \end{align}

(6)\begin{align} h_g(T) &= c_{p,g}(T-T_m). \end{align}

Here $c_{p,\alpha}$ is the specific heat capacity at constant pressure (J/kg·K) for phase α, 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} \lt c_{p,i}$ or $c_{p,w}$ (see Table 1). Hence, we make the simplification that the gas phase contribution to the total enthalpy of the system is negligible. After substituting Eqns (4) and (5) into Eqn (3), H can be ultimately formulated as

(7)\begin{align} H &= \begin{cases}\rho c_{p,i} \phi_i (T-T_m), & {T \lt T_m} \quad (\textrm{or } H \leqslant 0) \\ \rho \phi_w L, &{T= T_m} \quad (\textrm{or } 0 \lt H \lt CL) \\ \rho \phi_w \left(c_{p,w} (T-T_m) + L \right), &{T \gt T_m} \quad (\textrm{or } H \geqslant C L) \end{cases}. \end{align}

Table 1. A summary of simplified thermodynamic properties as well as flow properties of water in porous ice used in present work

Here the maximum enthalpy limit for the three-phase region is the product of composition, C, and the latent heat of fusion, L because it is the enthalpy of the fully molten system at the melting point (see Fig. 1a). The boundaries of the three-phase region are not included in the region $0 \lt H \lt CL$ because it strictly refers to the three-phase region. From the above formulation, we classify three regions, where region 1 ( $H \leqslant 0$) is comprised of ice and gas, region 2 ( $0 \lt H \lt CL$) contains all three phases and region 3 ( $H \geqslant CL$) corresponds to a no-matrix state consisting of only water and gas phases.

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. H = 0 or $\mathcal{H}=0$ and H = CL or $\mathcal{H}=\mathcal{C}$. The three regions are labeled in panels a and e.

The temperature and volume fractions of water, ice and gas phases can be evaluated from composition, C, and enthalpy, H, as shown in Fig. 1a–d, respectively. As shown in Eqn (2), the volume fraction of gas, ϕ_g, only depends on the composition, C, as shown in Fig. 1d. Next we formulate the governing equations for this model corresponding to the two conserved variables.

2.2. Transport model

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

(8)\begin{align} \frac{\partial C}{\partial t} + \nabla \cdot (\textbf{q} \rho) &= 0 \end{align}

(9)\begin{align} \frac{\partial H}{\partial t} + \nabla \cdot (\textbf{q} \rho h_w - \overline{\kappa} \nabla T) &= 0, \end{align}

where q is the volumetric flux of water phase (m³/m $^2\cdot$s) relative to ice phase. The effective thermal conductivity of the mixture $\overline{\kappa}$, typically a weighted average of the thermal conductivities of the phases, κ_α (W/m·K).

2.3. Constitutive relations

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

(10)\begin{align} \textbf{q} = - \frac{\mathrm{k}(\varphi) k_{rw}(s)}{\mu} (\nabla p - \rho \textbf{g}) \end{align}

where $\mathrm{k}$ is the absolute permeability (m²) which is a function of porosity φ $(\varphi=\phi_w+\phi_g =1-\phi_i)$, the ratio of void volume to the bulk volume, p is water pressure (Pa), µ is the viscosity of water (Pa·s) and g is the acceleration due to gravity vector (m²/s). The relative permeability for multiphase flow, k_rw, display complex hysteresis (Blunt, Reference Blunt2017), but here we only consider the simplest case with power law dependence. The relative permeability of the water phase, k_rw, is a function of the water saturation, s, which is the ratio of water phase volume to void volume, $s= \phi_w/1-\phi_i$. We assume that the water phase becomes immobile below a certain residual water saturation, $s_{w r}$. 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} \lt s \lt 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 has been determined to be $s_{w r}\sim0.07$ from lysimeter (Colbeck, 1976) and calorimeter (Coléou and Lesaffre, Reference Coléou and Lesaffre1998) techniques. However, as the ice is water-wet with a near-zero contact angle at the ice–water–air interface (Knight, Reference Knight1971), the residual water saturation during saturation rise (imbibition) is zero due to hysteresis in the relative permeability–capillary pressure curve (Carlson, Reference Carlson1981; Blunt, Reference Blunt2017). The phenomenon of hysteresis has largely 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 unsaturated regions (Colbeck, Reference Colbeck1972), such that the spatial variations in the difference between the water and air pressure (capillary pressure) are negligible at the problem length scales. See Smith (Reference Smith1983), Shadab and Hesse (Reference Shadab and Hesse2022, Reference Shadab and Hesse2024) for a more detailed discussion on neglecting the capillary pressure term in the 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 (Colbeck, Reference Colbeck1972; Shadab and Hesse, Reference Shadab and Hesse2022). Plugging it in Eqn (10) eliminates the diffusive, pressure term. The volumetric flux of water, q, finally takes the gravity-driven form

(11)\begin{align} \textbf{q} = \frac{\mathrm{k}(\varphi) k_{rw}(s)}{\mu} \rho \textbf{g}. \end{align}

The absolute permeability of ice (m²), $\mathrm{k}$, and the relative permeability of water, k_rw, are assumed to be power laws (Kozeny, Reference Kozeny1927; Carman, 1937; Brooks and Corey, Reference Brooks and Corey1964; Bear, Reference Bear2013; Meyer and Hewitt, Reference Meyer and Hewitt2017) defined as

(12)\begin{align} \mathrm{k}(\varphi) &=\mathrm{k}_0 \varphi^m = \mathrm{k}_0 (1-\phi_i)^m , \end{align}

(13)\begin{align} k_{rw}(s) &= k^0_{rw} s^n = k^0_{rw} \left( \frac{\phi_w}{1-\phi_i}\right)^n , \end{align}

where $\mathrm{k}_0$ is a model constant, which can be considered as an absolute permeability (m²) when there is 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 Eqns (12) and (13) in Eqn (11) finally gives

(14)\begin{align} & \textbf{q}(\phi_i,\phi_w)\nonumber\\ &\quad = \begin{cases} \textbf{0}, &H\leqslant 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}} & 0 \lt H \lt CL \,\&\, {C \lt \rho}\\ \quad\! = K_h (1-\phi_i)^m \left( \frac{\phi_w}{1-\phi_i}\right)^n \hat{\textbf{g}}, & \end{cases} \end{align}

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. The symbol $K_h = \frac{\mathrm{k}_0 k^0_{rw}}{\mu}\rho g$ 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 (10–12) are only valid for a porous medium. Therefore, we will restrict our analysis to the porous media where φ < 1 and the flow is laminar, thus, Darcy’s law is applicable (Tek, Reference Tek1957).

2.4. Scaling

We nondimensionalize the model to make it scale-independent and find dominant terms governing the physics of the problem. The model is scaled using dimensionless variables for composition, $\mathcal{C}$, enthalpy, $\mathcal{H}$, temperature, $\mathcal{T}$, depth, ζ and time, τ, which are definedFootnote ¹ as

(15)\begin{equation} \mathcal{C} = \frac{C}{\rho}, \mathcal{H} = \frac{H}{\rho L}, \ \mathcal{T}=\frac{T-T_m}{T_m},\ \zeta = \frac{z}{\delta}, \ \text{and }\ \tau = \frac{t K_h}{\delta}. \end{equation}

Here the spatial coordinates (e.g. the depth coordinate z) are nondimensionalized by length scale of heterogeneity or the REV scale of the problem, δ. 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 Eqns (2) and (7), respectively, thus transform into the dimensionless forms

(16)\begin{align} \mathcal{C} &= \phi_i+\phi_w= 1-\phi_g , \end{align}

(17)\begin{align} \mathcal{H} &= \begin{cases} \mathcal{C} \, \mathrm{Ste} \, \mathfrak{c}_{p,r} \mathcal{T} , & \mathcal{H} \leqslant 0 \\ \phi_w, &0 \lt \mathcal{H} \lt \mathcal{C} \\ \mathcal{C} (\mathrm{Ste} \, \mathcal{T} + 1), & \mathcal{H} \geqslant \mathcal{C} \end{cases}, \end{align}

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 (17) and dimensional-specific enthalpy of water phase (5), the dimensionless temperature, $\mathcal{T}$, and dimensionless specific enthalpy of water phase, $\mathfrak{h}_w=h_w/L$, can be derived as

(18)\begin{align} \mathcal{T}(\mathcal{C},\mathcal{H}) &= \begin{cases} \frac{\mathcal{H}}{\mathcal{C} \, \textrm{Ste} \,\mathfrak{c}_{p,r}}, &{\mathcal{H} \leqslant 0} \\ 0, &0 \lt \mathcal{H} \lt \mathcal{C}\\ \frac{1}{\textrm{Ste}}\left(\frac{\mathcal{H}}{\mathcal{C}}-1 \right), &{\mathcal{H} \geqslant \mathcal{C}} \end{cases}, \qquad \textrm{and } \quad\nonumber\\ \mathfrak{h}_w (\mathcal{C},\mathcal{H}) &= \begin{cases} 0,& \mathcal{H} \leqslant 0 \\ 1 , & 0 \lt \mathcal{H} \lt \mathcal{C}\\ \frac{\mathcal{H}}{\mathcal{C}}, & \mathcal{H} \geqslant \mathcal{C} \\\end{cases}. \end{align}

Subsequently, the volume fractions of the phases, ϕ_α, and the porosity of the medium, φ, can be rewritten as functions of $\mathcal{C}$ and $\mathcal{H}$ as

(19)\begin{align} \phi_w(\mathcal{C},\mathcal{H}) &= \begin{cases}0,& \mathcal{H} \leqslant 0 \\ \mathcal{H} , & 0 \lt \mathcal{H} \lt \mathcal{C}\\ \mathcal{C}, & \mathcal{H} \geqslant \mathcal{C} \\\end{cases}, \quad\nonumber\\ \phi_i(\mathcal{C},\mathcal{H}) &= \begin{cases}\mathcal{C},& \mathcal{H} \leqslant 0 \\ \mathcal{C}-\mathcal{H} , & 0 \lt \mathcal{H} \lt \mathcal{C}\\ 0, & \mathcal{H} \geqslant \mathcal{C} \\\end{cases}, \end{align}

(20)\begin{align} \phi_g(\mathcal{C}) &= 1 - \mathcal{C} \quad \quad \textrm{and} \quad \quad\nonumber\\ \varphi(\mathcal{C},\mathcal{H}) &= \begin{cases}1-\mathcal{C},& \mathcal{H} \leqslant 0 \\ 1-\mathcal{C}+\mathcal{H} , & 0 \lt \mathcal{H} \lt \mathcal{C}\\ 1, & \mathcal{H} \geqslant \mathcal{C} \\\end{cases}. \end{align}

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 Fig. 1e–h, respectively. As shown in Eqn (2), the volume fraction of gas, ϕ_g, only depends on the dimensionless composition as illustrated in Fig. 1h.

The composition and enthalpy transport equations (8 and 9) thus take the dimensionless form

(21)\begin{align} \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, \end{align}

(22)\begin{align} \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}} - {\frac{\overline{\kappa}}{\kappa_w}} \frac{\textrm{Ste}}{Pe_{\mathcal{H}}}\nabla \mathcal{T} \right) = 0. \end{align}

Here κ_w is the thermal conductivity of the water phase. 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$, where $\alpha_T={\kappa_w}/{\rho c_{p,w}}$ (m²/s) is the thermal diffusivity of water. Moreover, the divergence and gradient operators are now scaled with inverse of characteristic depth, $1/\delta$.

While both conduction and advection can be important heat transport processes in firn (Shadab and others, Reference Shadab, Adhikari, Rutishauser, Grima and Hesse2024a), conduction does not affect the steady propagation of wetting and drying fronts analyzed here. This can be shown using thermodynamic and fluid flow parameters from Table 1. The value of Peclet number, $Pe_{\mathcal{H}}$, comes out to be $4\cdot 10^{3}$ for δ = 1 m. The ratio $\overline{\kappa}/\kappa_w \leqslant 1$ as its maximum value of unity is achieved when the REV only contains liquid water. The Stefan number Ste is a constant of value 3.43. Therefore, the value of $\frac{\overline{\kappa}}{\kappa_w}\frac{\text{Ste}}{Pe_{\mathcal{H}}}$ indicates about three orders of magnitude higher heat advection compared to heat conduction for gravity-driven infiltration in firn. Therefore, we can neglect the second-order heat conduction term in Eqn (22), which is also the necessary condition in the three-phase region ( $0 \lt \mathcal{H} \lt \mathcal{C}$ or $\mathcal{T}=0$) as $\nabla \mathcal{T} = 0$. Assuming local thermodynamic equilibrium, in the limit $Pe_{\mathcal{H}}\to \infty$, the system of dimensionless governing equations (21 and 22) then reduces to quasi-linear system of coupled hyperbolic equations,

(23)\begin{align} \frac{\partial \textbf{u}}{\partial \tau} + \nabla \cdot {\boldsymbol{f}}(\textbf{u})= \textbf{0} \end{align}

where $\textbf{u} = [\mathcal{C},\mathcal{H} ]^T$ is the vector of dimensionless conserved variables and ${\boldsymbol{f}}(\textbf{u}) = [ {\boldsymbol{f}}_{\mathcal{C}},{\boldsymbol{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

(24)\begin{align} {\boldsymbol{f}}:={\boldsymbol{f}}_{\mathcal{C}} = {\boldsymbol{f}}_{\mathcal{H}} = \begin{cases} \textbf{0},& \mathcal{H} \leqslant 0, \\ (1-\mathcal{C}+\mathcal{H})^m \left( \frac{\mathcal{H}}{1-\mathcal{C}+\mathcal{H}}\right)^n \hat{\textbf{g}}, & 0 \lt \mathcal{H} \lt \mathcal{C}.\end{cases} \end{align}

The above analysis shows the distinct system behaviors in the different regions based on fluxes (see Fig. 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 Eqn (24), it can be observed that the fluxes of enthalpy and composition in the system of governing equation (23) are identical in regions 1 and 2 defined by the symbol ${\boldsymbol{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 the form of latent heat (and specific heat) along with it.

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} \leqslant 0$) as well as region 3 comprising of three-phase region ( $0 \lt \mathcal{H} \lt \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}\geqslant \mathcal{C}$) is not considered in the present work.

In the next section, we will consider a simple problem of melt transport across a discontinuity to utilize the MOC for solving the system of hyperbolic PDEs (Lighthill and Whitham, Reference Lighthill and Whitham1955; LeVeque, Reference LeVeque1992) given in Eqn (23).

3.1. General structure of reaction fronts

Consider the following one-dimensional initial value problem with two constant states, known as a Riemann problem. See LeVeque Reference LeVeque1992 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 (ζ). In that case, the dimensionless flux vector (24) reduces to

(25)\begin{align} f: = \begin{cases} 0,& \mathcal{H} \leqslant 0 \\ (1-\mathcal{C}+\mathcal{H})^m \left( \frac{\mathcal{H}}{1-\mathcal{C}+\mathcal{H}}\right)^n, & 0 \lt \mathcal{H} \lt \mathcal{C}\end{cases} \end{align}

The system of dimensionless composition and enthalpy conservation equations (23) can be written in one-dimensional depth coordinates, ζ, as

(26)\begin{align} \textbf{u}_{\tau} + \textbf{f}(\textbf{u})_{\zeta} = \textbf{0}, \quad \tau \in \mathbb{R}^{+}, \quad \zeta \in \mathbb{R}, \end{align}

with initial conditions

(27)\begin{align} \textbf{u}=\begin{cases} \textbf{u}_l, &\zeta \lt 0 \\ \textbf{u}_r, &\zeta \gt 0 \end{cases} \end{align}

where the flux vector in ζ direction is $\textbf{f}(\textbf{u})=[f,f]^T$ and the subscripts τ and ζ refer to the partial derivatives with respect to the dimensionless time and depth, respectively. The subscripts l and r refer to the state of the system on the left and right sides of a discontinuity. During melt infiltration into firn the left (right) state corresponds to the top (bottom) layer around a discontinuity. An example of an initial discontinuity with left state, u_l, and right state, u_r, is shown in Fig. 3b and c. The flow is toward the direction of gravity, assumed to be in $+ \zeta$-direction. The solution to the Riemann problem for well-behaved systems of two coupled nonlinear PDEs is characterized by the formation of an intermediate state, u_i, bounded by two waves $\mathscr{W}_1$ and $\mathscr{W}_2$ (LeVeque, Reference LeVeque1992). This solution structure, observed in Fig. 4, can be represented as

(28)\begin{align} \textbf{u}_l\xrightarrow{\mathscr{W}_1}\textbf{u}_i \xrightarrow{\mathscr{W}_2}\textbf{u}_r \end{align}

Figure 4. Solution of the Riemann problem introduced in Figure 3. (a) Evolution of dimensionless composition, $\mathcal{C}$, in space, ζ, and time, τ. (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, u_i, corresponds to a state between the fronts. The system (26) can be recast into a quasilinear form using chain rule as

(29)\begin{align} {\textbf{u}}_\tau + \nabla_{\textbf{u}} \textbf{f}(\textbf{u}) \, \textbf{u}_\zeta = \textbf{0} \end{align}

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

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

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, and which are given in Appendix A. The system of advection equations (29) 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 step change allows the introduction of the similarity variable

(31)\begin{equation} \eta = \frac{\zeta}{\tau} \end{equation}

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

(32)\begin{align} (\textbf{A} - \lambda_p \textbf{I})~\textbf{r}_p = \textbf{0}, \quad p\in \{1,2 \}, \end{align}

where the flux gradient is $\textbf{A}=\nabla_{\textbf{u}} \textbf{f}(\textbf{u})$ and the eigenvector is $\textbf{r}_p=\mathrm{d} \textbf{u}/\mathrm{d} \eta$ corresponding to the eigenvalue λ_p. Here the eigenvalues λ ₁ and λ ₂ are the characteristic propagation speeds of the waves $\mathscr{W}_1$ and $\mathscr{W}_2$, respectively. The associated eigenvectors, $\textbf{r}_p=\mathrm{d} \textbf{u}/\mathrm{d} \eta$, give the pathways through the $\mathcal{C}-\mathcal{H}$ plane, also referred to as the hodograph plane, that satisfy the conservation equations (see Fig. 5b, e.g.).

Constant solutions of the conservation laws (29) satisfy (32) trivially because $\mathrm{d} \textbf{u}/\mathrm{d} \eta=0$. Solutions of the conservation laws (29) that vary continuously must instead satisfy the eigenvalue problem (32). Finally, discontinuities in the solution of (29) must satisfy the Rankine–Hugoniot (R–H) jump condition (LeVeque, Reference LeVeque1992), which is derived from the discrete conservation of mass and enthalpy around the discontinuity, and given by

(33)\begin{align} \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}]}, \end{align}

where $\Lambda_{\mathscr{S}}$ is the shock speed, $[\,\cdot\,]$ 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 an intermediate state.

3.3. Construction of the solution in the $\mathcal{C}-\mathcal{H} $ hodograph plane

The self-similar solutions are constructed by identifying directions in the $\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 directions is determined by the nonlinear algebraic system of equations arising from the R–H jump condition (33), described by shock fronts. First, we consider a system 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 (Shadab and Hesse, Reference Shadab and Hesse2022). Although there can be at most one moving wave for simple cases where the medium remains unsaturated ( $C \lt \rho$, $\mathcal{C} \lt 1$ or $\phi_g \gt 0$), there can be two moving waves when a fully saturated region appears, i.e. $\mathcal{C}=1$. Below we enumerate 12 distinct self-similar solutions to the Riemann problem that describe a variety of firn hydrological processes and summarize them in Table 2.

Table 2. Summary of all analytic solutions presented in this paper along with related works that either studied or observed the corresponding scenario.

3.3.1. Region 2 only (three-phase region)

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

(34)\begin{align} \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}. \end{align}

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, λ ₂, and (b) propagation velocity of second front with respect to the melt given by $\phi_w\lambda_2$ in the $\mathcal{C}-\mathcal{H}$ hodograph plane for m = 3 and n = 2. The slow path r₁ and fast path r₂ 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 λ ₂ 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 Fig. 5a for m = 3 and n = 2. Therefore, all solutions governed by hyperbolic PDEs (26) 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 Fig. 5b. As expected, the propagation velocity $\phi_w \lambda_2$ decreases toward the lower boundary of the three-phase region ( $\mathcal{H}=0$). These eigenvalues yield two corresponding, linearly independent eigenvectors in the $\mathcal{C}-\mathcal{H}$ hodograph plane given by

(35)\begin{align} \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}. \end{align}

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

(36)\begin{align} \frac{\mathrm{d} \textbf{u}}{\mathrm{d} \eta} &= \frac{1}{\nabla_{\textbf{u}} \lambda_p \cdot \textbf{r}_p} \textbf{r}_p \end{align}

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

(37)\begin{align} \textbf{u}(\textbf{u}_0,\eta) &= \textbf{u}_0 + \int_{\lambda_p(\textbf{u}_0)}^\eta \frac{1}{\nabla_{\textbf{u}} \lambda_p \cdot \textbf{r}_p} \textbf{r}_p ~\mathrm{d} \eta'. \end{align}

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

(38)\begin{align} \mathcal{C}&=1+ \mathcal{H}+\mathfrak{C} \mathcal{H}^{-\frac{n}{m-n}}. \end{align}

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

(39)\begin{align} \mathcal{C} &= \mathcal{H} + \mathfrak{C} . \end{align}

The speed of second characteristic λ ₂ is nonnegative and increases monotonically in the direction of integral curves corresponding to the second eigenvector r₂ (fast paths) in three-phase region for n > 1 (see Lemma A.1 in Appendix A and Fig. 5b). Additionally, the slow path corresponds to constant flux contours as $\lambda_1=0$ and the fast path corresponds to constant porosity, φ, contours (Lemma A.2 in Appendix A).

Figure 6. The simple solutions of a 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 the hodograph plane and their corresponding self-similar analytical solutions for (b) dimensionless composition and (c) dimensionless enthalpy with dimensionless velocity η. 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.

Solutions in the three-phase region

We will now discuss the different analytical solutions within the three-phase region, tailored to glaciological scenarios. All scenarios are summarized in Table 2. The discussion below assumes that u₀ is the left state, u_l, and describes the set of permissible right states u_r. We begin with cases that lead to solutions with a single front before considering cases leading to two fronts and the formation of an intermediate state.

(a) Stationary linear reaction front (Case I)

This case resembles a steady meltwater flux into a temperate firn with a step reduction in porosity at shallow depth (Fig. 6d–f). The integral curves associated with λ ₁ and r₁, known as the first characteristic field ( $\lambda_1,\textbf{r}_1$), are constant flux lines in the $\mathcal{C}-\mathcal{H}$ hodograph plane. Any right state u_r, along the integral curve, connected to u₀ by a stationary discontinuity is a weak solution of Equation (29). 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 (38) (constant flux lines) satisfying

(40)\begin{align} \frac{\varphi_l}{\varphi_r}&= \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). \end{align}

The complete solution in this case takes the form

(41)\begin{align} \textbf{u} &= \begin{cases}\textbf{u}_l, \quad & \zeta \lt 0 \\ \textbf{u}_r, \quad & \zeta \gt 0\end{cases}. \end{align}

Fig. 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 ζ < 0 to 55.3% in ζ > 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 Fig. 6d–f at different times.

(b) Moving nonlinear reaction front

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

(42)\begin{align} \mathcal{C}_r&=\mathcal{C}_l + ( \mathcal{H}_r-\mathcal{H}_l). \end{align}

(i) Rarefaction wave (Case II)

This case resembles a sudden drop of meltwater influx into a temperate firn with constant porosity (Fig. 6g–i). If the characteristic speed λ ₂ varies smoothly from left to right state, any right state u along the fast path is connected to left state u₀ by a continuously varying saturation front (Fig. 6a, blue lines). The propagation velocity, λ ₂, 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 (Clark and others, Reference Clark, Nijssen and Luce2017). Rarefaction waves are a weak solution of Eqn (32) if the resultant profile of u is single-valued. This condition is satisfied if u lies on the branch of the integral curve r₂ emanating from u₀ in the direction of increasing λ ₂ (see Figs. 5 and 6a, blue lines). The analytical solution concerning the self-similar variable η for a rarefaction wave on an integral curve can be evaluated from Eqn (37), which comes out to be

(43)\begin{align} \mathcal{H}&= \sqrt[n-1]{\frac{\eta}{n (1-\mathcal{C}_l + \mathcal{H}_l)^{m-n}}}\quad \text{and} \end{align}

(44)\begin{align} \mathcal{C} &= \mathcal{C}_l -\mathcal{H}_l + \sqrt[n-1]{\frac{\eta}{n (1-\mathcal{C}_l + \mathcal{H}_l)^{m-n}}}. \end{align}

The final solution in this case takes the form

(45)\begin{align} \textbf{u} &= \begin{cases}\textbf{u}_l, \quad & \eta \lt \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) \lt \eta \lt \lambda_2(\textbf{u}_r) \\\textbf{u}_r, \quad & \eta \gt \lambda_2(\textbf{u}_r)\end{cases} \end{align}

where the speed of the second characteristic $\lambda_2(\cdot)$ is evaluated from Eqn (34). An example of this case is when meltwater flux instantly drops, leading to a smoothing drainage front. Fig. 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 Fig. 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 Eqn (45) are m = 3 and n = 2. This drying/rarefaction front has been studied by Clark and others Reference Clark, Nijssen and Luce(2017) and was observed in models studying the Dye-2 site in Greenland on 12 August 2016 after the meltwater flux ceased (Samimi and others, Reference Samimi, Marshall and MacFerrin2020; Vandecrux and others, Reference Vandecrux2020; Colliander and others, Reference Colliander2022; Shadab and others, Reference Shadab, Adhikari, Rutishauser, Grima and Hesse2024a).

(ii) Shock wave (Case III)

This case resembles a sudden increase of meltwater influx into a temperate firn with constant porosity (Fig. 6j–l). If the right state u lies on the opposite branch of the integral curve (shown by red line in Fig. 6a, e.g.), 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 u₀ by a discontinuous reaction front that propagates with a velocity, $\Lambda_{\mathscr{S}}(\textbf{u}_0,\textbf{u})$, which can be calculated from the R–H jump condition (33) 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 u₀ 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 r₂ (39) from the Hugoniot jump condition (33), since the flux of enthalpy and composition are the same. The dimensionless velocity of the shock from Eqn (33) is then

(46)\begin{align} \Lambda_{\mathscr{S}}(\textbf{u}_+,\textbf{u}_-) = \frac{\mathrm{d} \zeta}{\mathrm{d} \tau} = \frac{f(\textbf{u}_+)-f(\textbf{u}_-)}{\mathcal{H}_+-\mathcal{H}_{-}} = \frac{f(\textbf{u}_+)-f(\textbf{u}_-)}{\mathcal{C}_+-\mathcal{C}_{-}}, \end{align}

where the subscript − refers to the left (top) state and + is the right (bottom) state for this particular configuration where $f(\textbf{u}_l) \gt f(\textbf{u}_r)$. The final solution thus takes the form

(47)\begin{align} \textbf{u} &= \begin{cases}\textbf{u}_l, \quad & \zeta / \tau \lt \Lambda_{\mathscr{S}}(\textbf{u}_r,\textbf{u}_l) \\ \textbf{u}_r, \quad & \zeta / \tau \gt \Lambda_{\mathscr{S}}(\textbf{u}_r,\textbf{u}_l) \end{cases}. \end{align}

When meltwater flux instantly increases, it leads to a sharp wetting front propagating in the direction of gravity. Fig. 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 ϕ_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 Fig. 6a–c. The wetting front has been discussed (Colbeck, Reference Colbeck1972; Humphrey and others, Reference Humphrey, Harper and Pfeffer2012; Clark and others, Reference Clark, Nijssen and Luce2017; Meyer and Hewitt, Reference Meyer and Hewitt2017, e.g.) and observed in models (Vandecrux and others, Reference Vandecrux2020; Samimi and others, Reference Samimi, Marshall, Vandecrux and MacFerrin2021; Colliander and others, Reference Colliander2022, e.g.) and field observations including at the Dye-2 site in Greenland on 9 August 2016 when meltwater percolates in temperate firn (Heilig and others, Reference Heilig, Eisen, MacFerrin, Tedesco and Fettweis2018; Samimi and others, Reference Samimi, Marshall and MacFerrin2020).

(c) Two fronts with an intermediate state

The solution profile contains a single reaction front if u_l and 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, 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 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 are parallel to the upper boundary of the three-phase region ( $\mathcal{H}=\mathcal{C}$). The two possible intermediate states are given by the intersections of the integral curves emanating from u_l and u_r (see Fig. 7a, e.g.). Only one intersection yields a physically realistic single-value solution. The correct intersection is selected by requiring that the propagation speed increases monotonically from u_l and u_r. A single-valued solution is ensured if and only if u_l and u_i lie on the slow path first and then u_i is connected to u_r along the fast path. The reactive melt transport system considered here only allows two solutions for this case:

(48)\begin{align} \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, \end{align}

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 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 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 the hodograph plane and their corresponding self-similar analytical solutions for (b) dimensionless composition and (c) dimensionless enthalpy with dimensionless velocity η. 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)

This case resembles a sudden drop in meltwater flux inside temperate firn where porosity also reduces with depth (Fig. 7d–f). In this case, $f(\textbf{u}_l) \lt f(\textbf{u}_r)$ and the resulting first characteristic wave is a contact discontinuity, $\mathscr{C}_1$, which satisfies Eqn (38) for left state u_l and the intermediate state u_i. The second characteristic wave is a rarefaction which is governed by Eqns (39), (43) and (44) for intermediate state u_i and right state u_r. Combining all these equations results in a nonlinear algebraic equation to evaluate $\mathcal{H}_i$ that is given by

(49)\begin{align} & {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}\nonumber\\ &\quad =\left({1+ \mathcal{H}_l-\mathcal{C}_l}\right) \left(\frac{\mathcal{H}_i}{\mathcal{H}_l}\right)^{-\frac{n}{m-n}} \end{align}

which can be re-written in terms of porosities, φ, as

(50)\begin{align} {\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}&={\varphi_l} \left(\frac{\mathcal{H}_i}{\mathcal{H}_l}\right)^{-\frac{n}{m-n}}. \end{align}

Next, the composition at intermediate state, $\mathcal{C}_i$, can be computed from Eqns (39) and (34) which corresponds to the fast path and the speed of the second characteristic, λ ₂, respectively. The final solution in this case takes the form

(51)\begin{align} \textbf{u} &= \begin{cases}\textbf{u}_l, \quad & \zeta \lt 0 \\ \textbf{u}_i, \quad & 0 \lt \zeta / \tau \lt \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) \lt \eta \lt \lambda_2(\textbf{u}_r) \\\textbf{u}_r, \quad & \eta \gt \lambda_2(\textbf{u}_r)\end{cases}. \end{align}

Physically, the first contact discontinuity represents the constant flux of water entering the bottom, low porosity layer (Fig. 7e and f). The rarefaction shows the drainage of the wetter firn due to gravity. Fig. 7a 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\%$ 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 Fig. 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\%$. Fig. 7b and 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. Fig. 7d–f show 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 with a sudden rise in meltwater flux (Fig. 7g–i), and as such $f(\textbf{u}_l) \gt f(\textbf{u}_r)$. As a result, the first characteristic wave is a contact discontinuity, $\mathscr{C}_1$, that satisfies Eqn (38) for the left state u_l and intermediate state 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 (33) for the intermediate state and right state. Combining Eqns (38) and (39) gives a simple relation for dimensionless enthalpy at intermediate state, $\mathcal{H}_i$, to be

(52)\begin{align} \mathcal{H}_i&=\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}}. \end{align}

Then Eqn (39) for the intermediate state, u_i, and right state, u_r, provides the value of dimensionless composition at intermediate state, $\mathcal{C}_i$. The final solution in this case takes the form

(53)\begin{align} \textbf{u} &= \begin{cases}\textbf{u}_l, \quad & \zeta \lt 0 \\ \textbf{u}_i, \quad & 0 \lt \zeta / \tau \lt \Lambda_{\mathscr{S}}(\textbf{u}_i,\textbf{u}_r) \\ \textbf{u}_r, \quad & \zeta / \tau \gt \Lambda_{\mathscr{S}}(\textbf{u}_i,\textbf{u}_r)\end{cases}. \end{align}

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 (Fig. 7h and i). Fig. 7a 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 Fig. 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 Fig. 7g–i). Fig. 7b and 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)

This case resembles a temperate firn with a step decrease in porosity at depth receiving a sudden increase of meltwater influx (Fig. 8d–f), which results in ponding at the porosity contrast and a rising water table. 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} \lt 1$, as shown in Fig. 8a. In other words, the bottom layer is unable to accommodate the flux of meltwater from the top layer, as discussed for soils in Shadab and Hesse (Reference Shadab and Hesse2022). 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) PDEs (see Shadab and Hesse, Reference Shadab and Hesse2022 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 (24) anymore.

Nevertheless, we can still construct a full analytical solution to this problem using the extended kinematic wave approximation proposed by Shadab and Hesse (Reference Shadab and Hesse2022). In this case, the solution consists of three waves including a backfilling shock moving upwards (Fig. 8d–f), 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@@@

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

(54)\begin{align} \mathcal{C}_{i_1} = 1 \quad \text{and} \quad \mathcal{H}_{i_1}=1-\mathcal{C}_l+\mathcal{H}_l. \end{align}

This state is observed right next to the left state and can be considered as the rising perched water table in the region ζ < 0. The speed of this backfilling, upper shock in dimensionless form is again given by the R–H jump condition as

(55)\begin{align} \Lambda_{\mathscr{S}^*_1} = \frac{\mathrm{d} \zeta_U}{\mathrm{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}} \lt 0, \end{align}

where ζ_U is the location of the upper shock. The shock moves upwards due to choking as the numerator is positive, because flux $f(\textbf{u}_l)$ is larger 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.

Figure 8. Formation of a fully saturated region in temperate firn (Case VI): (a) Construction of solution in the hodograph plane and their corresponding self-similar analytical solutions for (b) dimensionless composition and (c) dimensionless enthalpy with dimensionless velocity η. 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) τ = 0, (e) 0.5, (f) 1.0.

Next, the first intermediate state i ₁ is connected to the second intermediate state i ₂ through a stationary jump $\mathscr{J}_2$ at the location of initial jump at ζ = 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 (Shadab and Hesse, Reference Shadab and Hesse2022). 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 ₂, are provided by the right state. The second intermediate state also lies at $\mathcal{C}=1$ on the fast path (39) (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,

(56)\begin{align} \mathcal{C}_{i_2} = 1 \quad \text{and} \quad \mathcal{H}_{i_2} = 1 - \mathcal{C}_r + \mathcal{H}_r . \end{align}

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

(57)\begin{align} \Lambda_{\mathscr{S}_3} = \frac{\mathrm{d} \zeta_L}{\mathrm{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} \gt 0 \end{align}

where ζ_L is the dimensionless location of the lower shock. Similar to a two-layer soil discussed in Shadab and Hesse (Reference Shadab and Hesse2022), 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

(58)\begin{align} q_s(\tau) =\frac{\zeta_U(\tau) - \zeta_L(\tau)}{\frac{\zeta_U(\tau)}{K_{i_1}}-\frac{\zeta_L(\tau)}{K_{i_2}}}, \end{align}

where $K_{i_1}$ and $K_{i_2}$ are the dimensionless saturated hydraulic conductivities at the first and second intermediate states given by

(59)\begin{align} 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 \end{align}

using Eqns (54) and (56). 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 Eqns. (57) and (55) along with the definition of flux in the saturated region (58) gives the location of the shocks and the flux in the saturated region. Similar to the case of two-layered soils in Shadab and Hesse (Reference Shadab and Hesse2022), to find the analytic value of $q_s(\tau)$, the ratio of shock speeds can be considered a constant as an ansatz given by

(60)\begin{align} \frac{\Lambda_{{\mathscr{S}^*_1}}}{\Lambda_{{\mathscr{S}_3}}} = \frac{\zeta_U(\tau)}{\zeta_L(\tau)} = \mathfrak{R}\quad\mbox{for} \quad 0 \lt \tau \leqslant \tau_p, \end{align}

where $\mathfrak{R}$ is the constant ratio of shock speeds which is negative and τ_p is the dimensionless time of ponding when the upward moving shock reaches 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 (60), shock speeds (55) and (57), and flux in the saturated region (58), $\mathfrak{R}$ comes out to be the solution of a quadratic equation

(61)\begin{align} \mathfrak{R}(\textbf{u}_l,\textbf{u}_r) &= \frac{-b - \sqrt{b^2 - 4ac}}{2a} \quad \text{where}\\ a &= \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\right.\nonumber\\ &\quad\left. \left(\frac{\mathcal{H}_r}{1-\mathcal{C}_r+\mathcal{H}_r} \right)^n \right], \nonumber \\ b &= - \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]\nonumber\\ &\quad + \left(\frac{\mathcal{H}_l}{1-\mathcal{C}_l+\mathcal{H}_l} \right)^n - 1 \quad \text{and } \nonumber \\ c &= 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 .\nonumber \end{align}

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, 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 the hodograph plane and their corresponding self-similar analytical solutions for (b) dimensionless composition and (c) dimensionless enthalpy with dimensionless velocity η. 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$.

Subsequently, substituting Eqns (59) and (60) in (58) gives the time-invariant dimensionless flux in the saturated region q_s during $0 \lt \tau \leqslant \tau_p$ as

(62)\begin{align} 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}. \end{align}

The time and space invariant flux q_s leads to the result using Eqns (55) and (57) that the shock speeds are constant for $0 \lt \tau \leqslant \tau_p$. Therefore, the shock locations before ponding vary linearly with time. Lastly, the full solution for this case can be summarized as

(63)\begin{align} \textbf{u} &= \begin{cases}\textbf{u}_l, \quad & \zeta \lt \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} \lt \zeta / \tau \lt 0\\ \textbf{u}_{i_2}=\begin{bmatrix}1\\1-\mathcal{C}_r+\mathcal{H}_r \end{bmatrix}, \quad & 0 \lt \zeta / \tau \lt \Lambda_{\mathscr{S}_3} \\ \textbf{u}_r, \quad & \zeta / \tau \gt \Lambda_{\mathscr{S}_3} \end{cases}. \end{align}

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. In terms of firn processes, the lower (low porosity) 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 (Fig. 8e and f). In this example, 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 (Fig. 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 (Fig. 8d–f). The speeds of both fronts $\mathscr{S}_1^*$ and $\mathscr{S}_3$ are constant before ponding occurs.

3.3.2. Region 1 only (ice and gas region), Case VII

While the previous cases consisted of temperate firn, this case resembles a cold firn with a step reduction in porosity at depth (e.g. porous firn on top of less porous firn or glacier ice), and no meltwater influx, resulting in a dry, static system (Fig. 9d–f). In region 1 ( $\mathcal{H}\leqslant 0$), since both fluxes are zero the system is not strictly hyperbolic and leads to a single wave. As the characteristic speed, λ_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 will not be any transport between the two states since the fluxes of both composition and enthalpy are zero on either side. In other words, u_l will be connected to u_r by a stationary contact discontinuity $\mathscr{C}$ as

\begin{align*} \textbf{u}_l\xrightarrow{\mathscr{C}}\textbf{u}_r. \end{align*}

The resulting solution thus takes the form

(64)\begin{align} \textbf{u} &= \begin{cases}\textbf{u}_l, \quad & \zeta \lt 0 \\ \textbf{u}_r, \quad & \zeta \gt 0\end{cases}. \end{align}

In this example, a temperate, 40% porous firn lies 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 Fig. 9a–c (blue lines) and d–f. However, in reality, the firn may compact due to overburden (Cuffey and Paterson, Reference Cuffey and Paterson2010) 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 a sharp transition at ζ = 0 due to the absence of heat conduction.

3.4. Solutions transitioning between regions

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

This case resembles a dry, porous cold firn layer on top of wet temperate firn (Fig. 9g–f), without additional meltwater influx. 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 lower boundary of the three-phase region ( $\mathcal{H}=0$) where the intermediate state lies (see Fig. 9a, e.g.), i.e.

(65)\begin{align} \mathcal{H}_i = 0. \end{align}

Moreover, the intermediate state, $\mathcal{C}_i$, lies on the fast path (constant porosity line) in region 2 satisfying Eqn (39) which gives

(66)\begin{align} \mathcal{C}_i &= \mathcal{C}_r - \mathcal{H}_r . \end{align}

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 Eqn (51).

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 Fig. 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$ (Fig. 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 Fig. 9b–c (red lines) and 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 corresponds to temperate and wet firn overlying initially cold firn (see Figs 10–12). 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} \lt 0$) which lies only in the right state with a sharp transition from the left or intermediate state as heat conduction is not considered in this model. There can be four scenarios 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, 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 the hodograph plane and their corresponding self-similar analytical solutions for (b) dimensionless composition and (c) dimensionless enthalpy with dimensionless velocity η. 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)

This case corresponds to a sudden increase in meltwater flux into temperate, wet firn overlying cold firn, resulting in meltwater percolation into the deeper layers and formation of a frozen fringe (Fig. 10d–f). 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 Fig. 10a, blue line). Note that the extended fast paths are not the constant porosity contours in region 1 (Fig. 1c and g). Mathematically, the left and right states are connected by the relation

(67)\begin{align} \mathcal{C}_r&=\mathcal{H}_r+ \mathcal{C}_l-\mathcal{H}_l . \end{align}

The shock speed $\Lambda_{\mathscr{S}}$ can then be evaluated using the Rankine-Hugoniot condition (46). The final solution of this case is provided in Eqn (47). 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. Fig. 10a–c (blue line) show an example of the solution where temperate firn with 45% porosity (φ_l) and 25% LWC ( $\phi_{w,l}$) initially lies on a cold layer of porosity 50% at T $=-15.84^\circ$C (Fig. 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$ (Fig. 10a–c, blue line). A refreezing shock $\mathscr{S}$ that results moves downwards while warming up the surrounding snow by partially refreezing (Fig. 10d–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 is similar to Case IX (wet firn overlying cold firn), but with a decreased porosity in the underlying cold firn leading to the formation of a frozen fringe (Fig. 10g–i). Here, 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 u_l in region 2 is more porous and connected to an intermediate state, 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 Fig. 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 Eqn (38) for left state u_l and intermediate state u_i. Since the second characteristic lies on the Hugoniot locus, a shockwave $\mathscr{S}_2$ results, that satisfies the Hugoniot-jump condition (33) for the intermediate state u_i and right state u_r. Combining all these equations gives a simple relation for dimensionless enthalpy and composition at the intermediate state as

(68)\begin{align} \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 . \end{align}

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

(69)\begin{align} \Lambda_{\mathscr{S}_2} = \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} \nonumber\\ &\qquad\quad \textrm{(since}\ f(\textbf{u}_r)=0\textrm{)}. \end{align}

The full solution is given in Eqn (53) with shock speed (69). Since $\mathcal{H}_r \lt 0$, $\mathcal{H}_i \gt 0$ and $f(\textbf{u}_r)=0$, the speed of the refreezing front $\mathscr{S}_2$ is slower than the temperate firn case (Case V) due to refreezing. Fig. 10a–c (red line) show 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 (Fig. 10g–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 Meyer and Hewitt Reference Meyer and Hewitt(2017) who used the dimensional form of Eqn (69) with distinct densities of ice and water phases to explain the field data from Humphrey and others (Reference Humphrey, Harper and Pfeffer2012).

Figure 11. Formation of a fully saturated region when right state lies in region 1: (a) Construction of solution in the hodograph plane and their corresponding self-similar analytical solutions for (b) dimensionless composition and (c) dimensionless enthalpy with dimensionless velocity η. 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) τ = 0, (e) 0.1, (f) 0.2.

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

This is similar to Case X, where wet and temperate firn layer lies on top of cold and less porous firn, and the system receives a sudden increase in meltwater flux (Fig. 11d–f). However, the temperature of the underlying cold firn is reduced, and the liquid water content in the above temperate layer increases which leads to the formation of a rising perched water table. The slow path originating from the left state and the extended fast path emanating from the right state do not intersect where $\mathcal{C} \lt 1$, making it different from 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 analytic solution (63) is the same as provided for the temperate firn (Case VI) consisting of a backfilling shock $\mathscr{S}^*_1$ moving upwards with speed (55), a stationary jump $\mathscr{J}_2$ at the location of initial jump ζ = 0 and a ‘refreezing’ front $\mathscr{S}_3$ moving downwards with speed provided in Eqn (57). The flux in the saturated region q_s is again provided by Eqn (58). Note that the flux at the right state is $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 Eqn (62) 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 Eqn (61), now simplifies to

(70)\begin{align} \mathfrak{R}(\textbf{u}_l,\textbf{u}_r) &= \frac{-b - \sqrt{b^2 - 4ac}}{2a} \quad \text{where}\\ a &= \left(\frac{1-\mathcal{C}_l}{1-\mathcal{C}_r} \right), \nonumber \\ b &= - \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 } \nonumber \\ c &= 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 .\nonumber \end{align}

Fig. 11 shows an example of this case where the left state corresponds to wet and temperate layer with 80% porosity (φ_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 a rising perched water table, is almost four times faster than the lower, refreezing front (see Fig. 11b–f). Below the initial jump at ζ = 0 the porosity is further reduced by refreezing, leading to the formation of a frozen fringe (Fig. 11e–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 Fig. 11e–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 and 3-Contact discontinuity (formation of an impermeable ice layer, Case XII)

This final case captures the formation of an impermeable ice layer through advection of meltwater and associated latent heat. An ice layer forms through the construction in Case XI, if the extended fast path ( $\mathcal{H}=\mathcal{C}+\mathfrak{C}$) emanating from the right state reaches $\mathcal{C}\geqslant 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 an ice layer will form if and only if the right state satisfies

(71)\begin{equation} 1-\mathcal{C}_r+\mathcal{H}_r\leqslant 0. \end{equation}

Figure 12. Impermeable ice layer formation: (a) Construction of solution in the hodograph plane and their corresponding self-similar analytical solutions for (b) dimensionless composition and (c) dimensionless enthalpy with dimensionless velocity η. 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 layer. The grey region corresponds to $1-\mathcal{C}+\mathcal{H}\leqslant 0$ where the right state resides to cause impermeable ice layer formation. The evolution of the volume fractions of the three phases in the resulting system at dimensionless times (d) τ = 0, (e) 0.02, (f) 0.04.

The exact location of the right state does not matter if it lies in the region of ice layer formation (grey region in Fig. 12a). In other words, an impermeable ice layer via heat advection will form if and only if the cold content of the firn exceeds the latent enthalpy of incoming meltwater. The mathematical condition (71) in dimensional form has been given in Humphrey and others (Reference Humphrey, Harper and Meierbachtol2021) for distinct densities of ice and water. In this limit, the solution for Case XI breaks down due to the formation of an impermeable ice layer. Fig. 12 shows the region of ice layer formation which requires either 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 ₁. The first saturated region is connected to the second intermediate state i ₂ through the jump $\mathscr{J}_2$. The second intermediate state i ₂ is the infinitesimally thin ice layer (see Fig. 12e–f) that blocks further meltwater percolation. The solution in this case can be written as

where the state solution is quantified as

(72)\begin{align} \textbf{u} &= \begin{cases}\textbf{u}_l, \quad & \zeta \lt \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} \lt \zeta / \tau \lt 0\\ \textbf{u}_{i_2}=\begin{bmatrix}1\\0 \end{bmatrix}, \quad & 0 \lt \zeta \lt \mathrm{d} \zeta \\ \textbf{u}_r, \quad & \zeta \gt \mathrm{d} \zeta \end{cases}, \end{align}

where $\mathrm{d} \zeta$ is the infinitesimally small thickness of the ice layer.

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 layer, 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

(73)\begin{align} \Lambda_{\mathscr{S}^*_1}(\textbf{u}_l) = \frac{\mathrm{d} \zeta_U}{\mathrm{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} \end{align}

which depends only on the left-state variables.

It can be deduced that either a cold firn (lower $\mathcal{H}$) or a nearly non-porous firn with $\mathcal{C}$ close to unity will induce the formation of an impermeable ice layer via refreezing caused by heat advection. For example, Fig. 12 shows such a solution when a wet and temperate layer of 70% porosity and 65% LWC lies 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$ (Fig. 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 layer forms as a result (see Fig. 12e–f) that blocks the flow of meltwater downwards. Since there is no meltwater percolating in the lower layer, the backfilling occurs very rapidly.

4.1. Meltwater infiltration into a multilayered firn—formation of a perched water table

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 Fig. 13a). 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 has not been studied in the firn 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 ζ = 1 (Fig. 13a). At time τ = 0, the meltwater is generated at the surface (ζ = 0), which keeps the liquid water content (denoted by LWC or ϕ_w) to 40% making the firn very wet at the surface, ζ = 0. The analytic solution of this problem can be constructed in the $\mathcal{C}-\mathcal{H} $ hodograph plane in two parts (Fig. 13b).

First, a wetting front $\mathscr{S}$ initially propagates downwards with a constant dimensionless speed $\Lambda_{\mathscr{S}}=0.28$, as given in Case III (Fig. 13c–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 by refreezing below the jump ζ > 1. This makes the flow pressure-driven instead of gravity-driven and enforces $q_s(\tau) \leqslant K_{i_2}$. This saturated region expands in both directions as a perched water table (upper shock shown by green dashed line in Fig. 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 firn. Thus, it shows a reduction in porosity from 30% to 21.2% for ζ > 1 behind the wetting front $\mathscr{S}_3$ at $\tau \gt \tau_s$ (Fig. 13c). 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 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 Fig. 13c–e.

These theoretical results show an excellent comparison with the numerical solutions shown in contour plots. The numerical solutions are obtained by solving the governing model (26) performed in the absence of capillary effects in between water and gas phases as well as heat conduction along with the treatment in the saturated region given by Shadab and Hesse (Reference Shadab and Hesse2024) and implementation of the enthalpy method in Shadab and others Reference Shadab, Adhikari, Rutishauser, Grima and Hesse(2024a). 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 (ζ = 0) is prescribed as the ‘Top condition’ in Fig. 13a whereas the bottom boundary condition is not required. 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 the 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.

4.2. Limitations of the present theory and further work

The unified kinematic wave theory makes necessary hydrologic and thermodynamic simplifications to make the system hyperbolic and amenable to solution using the method of characteristics. Some of these simplifications can be relaxed, but doing so will make the solutions more involved. On the hydrologic side, the theory neglects capillary forces in unsaturated firn, which may introduce significant errors at small spatial scales and grain sizes. Thus, it assumes that meltwater advection occurs faster as a wavefront than through diffusion via capillary suction. Capillary forces account for <10% of the total force, including gravity, when the meltwater flux is 10⁻⁸ m/s, although the percentage rapidly increases for smaller fluxes (Colbeck, Reference Colbeck1974b). The theory assumes that the density of water and ice are the same, which underestimates the porosity reduction due to refreezing when a wetting front percolates into a cold, dry firn. The kinematic model could easily be extended to account for different densities. The present model does not assume compaction of ice due to factors such as overburden, temperature, time and initial snow properties (e.g. crystal shape, packing geometry), which is present in both wet and dry firn (Bader, Reference Bader1954; Mellor, Reference Mellor1977; Herron and Langway, Reference Herron and Langway1980; Amory and others, Reference Amory2024). Lastly, the present theory assumes a one-dimensional flow in the direction of gravity with piecewise homogeneous firn and does not consider preferential flow which might enhance the speed and depth of meltwater percolation at higher melt fluxes (Wever and others, Reference Wever, Würzer, Fierz and Lehning2016; Vandecrux and others, Reference Vandecrux2020; Jones and others, Reference Jones, Moure and Fu2024).

On the thermodynamic side, the present theory assumes local thermodynamic equilibrium, which may not be a suitable assumption for rapid infiltration. Further, the theory neglects heat conduction, which may be dominant at smaller spatial scales and meltwater fluxes. Heat conduction leads to the formation of a conductive boundary layer ahead of the wetting front (Shadab and others, Reference Shadab, Adhikari, Rutishauser, Grima and Hesse2024a). The width of the thermal boundary layer is controlled by the Peclet number for enthalpy transport, as defined in Section 2.4. A high Peclet number typically leads to a thinner thermal boundary layer, indicating the dominance of heat advection over heat conduction.

However, the kinematic theory could be extended to incorporate several additional physical processes. It could be extended to account for hysteresis in relative permeabilities, which allows melt to enter dry snow without having to overcome a saturation threshold but to leave behind residual melt during drying. The method of characteristics can also be applied to study the kinematics of non-equilibrium processes, though the solutions are no longer self-similar. This could be used to analyze the percolation of melt that is not in thermal equilibrium with the ice or the coarsening of ice grains due to interaction with melt. The theory could be extended to add tracer transport that could be used to study changes in isotopic composition of the ice or dissolved gases in the melt and their effect on the climatic record in ice cores. Finally, the theory could be extended to incorporate salt and used to study the effect of colligative melting point depression in firn processes or sea-water intrusion into firn. As such, the kinematic theory is a useful tool to probe the physics of advective processes in firn and provide benchmark analytic solutions to a range of coupled non-linear processes.