# Numerical strategies for representing Richards' equation and its couplings in snowpack models

Kévin Fourteau, Julien Brondex, Clément Cancès, Marie Dumont

► To cite this version:

Kévin Fourteau, Julien Brondex, Clément Cancès, Marie Dumont. Numerical strategies for representing Richards' equation and its couplings in snowpack models. 2025. hal-05313594

HAL Id: hal-05313594

https://hal.science/hal-05313594v1

Preprint submitted on 14 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Numerical strategies for representing Richards' equation and its couplings in snowpack models

Kévin Fourteau<fnref n="1" />, Julien Brondex<fnref n="1" />, Clément Cancès<fnref n="2" />, and Marie Dumont<fnref n="1" />

**Correspondence:** Marie Dumont (marie.dumont@meteo.fr)

**Abstract.** The physical processes of heat conduction, liquid water percolation, and phase changes govern the transfer of mass and energy in snow. They are therefore at the heart of any physics-based snowpack model. In the last decade, the use of Richards' equation has been proposed to better represent liquid water percolation in snow. While this approach allows the explicit representation of capillary effects, it can also be problematic as it usually presents a large increase in numerical complexity and cost. This notably arises from the problem of applying a water retention curve in a fully-dry medium such as snow, leading to a divergence and degeneracy in Richards' equation. Moreover, the difficulty of representing both dry and wet snow in a single framework hinders the concomitant solving of heat conduction, phase changes, and liquid percolation. Rather, current models employ a sequential approach, which can be subject to non-physical overshoots. To treat these problems, we propose the use of a regularized water retention curve (WRC), that can be applied to dry snow. Combined with a variable switch technique, this opens the possibility of solving the energy and mass budgets in a fully consistent and tightly coupled manner, whether the snowpack contains dry and/or wet regions. To assess the behavior of the proposed scheme, we compare it to other implementations based on loose-coupling between processes and on the state-of-the-art strategies in snowpack models. Results show that the use of a regularized WRC with a variable switch greatly improve the robustness of the numerical implementation, consistently allowing timesteps greater or equal to 900 s, which results in faster and cheaper simulations. Furthermore, the possibility of solving the physical process in a fully-coupled and concomitant manner results in a slightly reduced error level compared to implementations based on the traditional sequential treatment. However, we did not observe any numerical oscillations and/or divergence sometimes associated with a sequential treatment. This indicates that a sequential treatment remains a potentially interesting tradeoff, favoring computational cost for a small decrease in precision.

## 1 Introduction

The wetting of snowpacks is a crucial stage in their evolution over time, with direct impacts on their environment. Notably, the intensity and timing of melt runoff have strong implications for the water availability in hydrological basins (Hock et al., 2019; Barnhart et al., 2020). Likely, the quantification of wet snow avalanche hazard relies on a precise determination of the liquid water content in snowpacks and of the depth of the percolation front (Baggi and Schweizer, 2009; Eckert et al., 2024). Liquid water percolation, also plays a crucial role in the formation of ice lenses and crusts (Wever et al., 2016; Quéno et al.,

<footnote n="1">¹Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Études de la Neige, Grenoble, France</footnote>

<footnote n="2">²Univ. Lille, CNRS, Inria, UMR 8524 - Laboratoire Paul Painlevé, 59000, Lille, France</footnote>
<figure><img src="image_3.png" /></figure>

2020), with direct implications for the local ecosystems (Tyler, 2010; Hansen et al., 2010). In this context, numerical snowpack models play a key role to assess the melting and wetting of snowpacks in various locations and under various conditions.

Historically, liquid water percolation in snowpack models has been implemented through the use of a so-called bucket-scheme, where snow layers are expected to retain liquid water until a certain threshold, after which all liquid water is instantaneously transferred downward (Bartelt and Lehning, 2002; Vionnet et al., 2012; Sauter et al., 2020). While this implementation is numerically efficient, it cannot capture certain effects, such as capillary barriers or capillary rise. To foster a more accurate representation of liquid water percolation, it has been proposed in the last decade to model it based on Richards' equation (Richards, 1931; Wever et al., 2014; D'Amboise et al., 2017). This more advanced description has notably been shown to better capture the timing associated with the wetting of the snowpack (Wever et al., 2015), although it does not address the challenge of simulating preferential flow, which is crucial to fully capture the complexity of liquid water percolation in snowpacks. (Wever et al., 2016; Moure et al., 2023).

In this framework, the simple bucket-scheme is replaced by a partial differential equation of liquid water mass budget with capillary effects. However, Richards' equation can be notoriously difficult to solve numerically, due to its non-linearities and potential degeneracy (Farthing and Ogden, 2017). Notably, it is usually associated with an adaptive time stepping strategy (e.g. Celia et al., 1990; Forsyth et al., 1995; Wever et al., 2014; D'Amboise et al., 2017), as the use of iterative methods to solve non-linear equations can fail at large timesteps. Therefore, the implementation of Richards' equation in snowpack models can imply a significant increase in numerical cost, hindering its broad use, in particular for simulations over large areas and long periods. For instance, the current implementation of Richards' equation in the detailed snowpack model Crocus (D'Amboise et al., 2017) is known to display signs of numerical instabilities and to require small timesteps of the order of 30 seconds or less (D'Amboise et al., 2017, and M. Lafaysse, personal communication, 2024), while the model is usually run with a 900 s timestep. This drastically hinder the routine use of Richards' equation in Crocus, as the two-order of magnitude increase in numerical cost can be considered as a too expensive trade-off. Moreover, due to the difficulty of representing both dry and wet snow in a unified framework for Richards' equation, current implementations in snowpack models assume that a small amount of liquid water is always present and need to modify the capillary behavior of snow in consequence (Wever et al., 2014; D'Amboise et al., 2017). This lack of unified treatment of dry and wet snow hinders the concomitant numerical solving of liquid water percolation with other physical processes (such as phase change or heat conduction). Rather, Richards' equation is solved sequentially (Wever et al., 2014; D'Amboise et al., 2017).

The main goal of this article is to investigate how the specific numerical implementation of liquid water percolation in snow affects the behavior of snowpack models, with a focus on their robustness and numerical cost. For this, we rely on various techniques proposed in applied mathematics, where the efficient solving of Richards' equation has been an active research subject (e.g. Forsyth et al., 1995; Sadegh Zadeh, 2011; Farthing and Ogden, 2017; Bassetto et al., 2020). Specifically, we explore whether the use of regularized capillary laws and the concomitant solving of Richards' equation with other physical processes are beneficial for the numerical behavior of snowpack models. The article is organized as follows. Section 2 presents a consistent system of equations describing energy and mass conservation in snowpacks that applies naturally to both dry and wet snow. Section 3 presents toy-models based on different numerical implementations of the heat and mass budget equations
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

and on the standard implementations used in state-of-the art detailed snowpack models as well. Simple test cases, representing three distinct situations of liquid water infiltration in snowpacks are presented in Section 4. Finally, the performance of the toy-models on these test cases are discussed in Section 5.

## 2 Deriving a consistent system of equations for the energy and mass budgets in snowpacks

While a real snowpack is a complex 3D structure composed of intertwined ice, air, and liquid water, such complexity is challenging to model. Rather, snowpack models rely on a macroscopic and 1D framework (e.g. Jordan, 1991; Bartelt and Lehning, 2002; Vionnet et al., 2012; Sauter et al., 2020). Here, macroscopic means that snow is not treated as an actual multiphasic medium, but rather as a homogeneous material (Torquato, 2002). In this framework, snow is characterized by macroscopic (sometimes referred to as effective) material properties, such as its thermal capacity (Calonne et al., 2014) or its hydraulic conductivity (Calonne et al., 2012).

The core of any snowpack model thus requires to determine and solve the equations governing the evolution of the macroscopic energy and mass contents of snow. This can be done from the first principles of energy and mass conservation, complemented by material laws characterizing the effective material properties. The goal of this section is to present a set of equations that govern the physical evolution of a snowpack and that apply both in dry and wet snow. Note that while this article assumes a 1D framework, as usually done in snowpack models, it could easily be transposed to a 2D or 3D configuration.

### 2.1 Snow internal energy conservation

As a first equation governing the evolution of a snowpack, we consider the energy conservation of snow, understood here as the combination of the ice matrix and the air within (and excluding potential liquid water). In the case of dry snow, and neglecting the influence of water vapor for simplicity, the temporal evolution of the snow energy is given by a classical conservation equation, i.e.:

$$ \partial_t h_s + \nabla \cdot F_{\text{cond}} = Q_{\text{abs}} \qquad (1) $$

where $h_s$ is the energy content of snow (expressed in J m⁻³), $F_{\text{cond}}$ the heat conduction flux (in W m⁻²), and $Q_{\text{abs}}$ a volumetric energy source (in W m⁻³), for instance the energy source due to shortwave absorption within the snowpack (van Dalum et al., 2019; Picard and Libois, 2024). Moreover, we classically assume that the heat conduction flux $F_{\text{cond}}$ follows Fourier's law:

$$ F_{\text{cond}} = -\lambda \nabla T \qquad (2) $$

where $\lambda$ is the thermal conductivity of snow (in W K⁻¹ m⁻¹) and $T$ its temperature (K). These equations, and the value of the snow thermal conductivity in terms of the microstructure, can for instance be derived from homogenization methods (e.g.
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

Boutin et al., 2010; Calonne et al., 2011, 2014; Bouvet et al., 2024). Moreover, the temperature and the energy content are related through the thermal capacity of snow:

$$h_s = \int_{T_0}^{T} c_p dT = c_p (T - T_0) \qquad (3)$$

where $T_0$ is a reference temperature, and $c_p$ is the volumetric thermal capacity of snow (in J m⁻³ K⁻¹; Tubini et al., 2021). Note that we have assumed for simplicity that $c_p$ does not depend on temperature, as regularly done in snowpack models (e.g. in the Crocus model; Vionnet et al., 2012). Also, we made the choice to take the fusion temperature of pure water as the reference temperature. This implies that the energy of snow is strictly negative below its fusion point, and vanishes at the fusion temperature. Using homogenization methods, the volumetric thermal capacity of snow can be shown (e.g. Calonne et al., 2014) to be

$$c_p = \phi_i c_i + (1 - \phi_i) c_a \approx \phi_i c_i \qquad (4)$$

where $\phi_i$ is the ice volumetric fraction, and $c_i$ and $c_a$ the volumetric thermal capacity of ice and air, respectively.

Specificities arise in the case of wet snow. Indeed, due to the presence of latent heat during the phase transition between ice and liquid water, the temperature of snow is blocked at the fusion temperature in case of ice-liquid water coexistence. The wet part of a snowpack thus becomes isothermal, and the internal thermal gradient vanishes, leading to a zero heat conduction flux. This phenomenon is already taken naturally into account in Eq 2, and thus does not require further development. A second point is that in the case of wet snow, input or removal of energy leads to phase change rather than temperature change. This phenomenon can be taken into account in Eq 1 by introducing an energy source term, related to the melting and refreezing of water (e.g. Tubini et al., 2021). The energy budget of ice matrix in snow is thus governed by:

$$\partial_t h_s - \nabla \cdot (\lambda \nabla T) = Q_{\text{abs}} + Q_{\text{freeze}} \qquad (5)$$

where $Q_{\text{freeze}}$ (in W m⁻³) is the rate of latent heat released/absorbed during the freezing/melting of water (taken positive for freezing).

## 2.2 Liquid water conservation

Since the heat budget and temperature of snow are directly related to the melting and freezing of water, it is necessary to also treat the liquid water budget. This can be done with the use of Richards' equation (Wever et al., 2014, 2015, 2016; D'Amboise et al., 2017). Note that in this paper we only consider "matric" water flow, and do not consider fast preferential flow (Vogel et al., 2000; Lewandowska et al., 2004; Wever et al., 2016). This limitation is discussed in Section 5.4. Under these conditions, the mass conservation of liquid water is given by:
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

$$ \rho_w \partial_t \theta_w + \rho_w \nabla \cdot F_w = -M_{\text{freeze}} \quad (6) $$

where $\theta_w$ is the volumetric liquid water content (LWC; expressed in m³ of water per m³ of snow), $\rho_w$ the density of liquid water (assumed constant; in kg m⁻³), $F_w$ is the liquid water flux (in m³ of water per m² of snow per s), and $M_{\text{freeze}}$ is the rate of freezing water (counted positive in the case of freezing and expressed in kg of water per m³ of snow per s). The rate of freezing water $M_{\text{freeze}}$ is directly related to the rate of energy release/absorption $Q_{\text{freeze}}$ of Eq 5 through $Q_{\text{freeze}} = L_{\text{fus}}M_{\text{freeze}}$, with $L_{\text{fus}}$ the specific enthalpy of fusion of water (in J kg⁻¹).

The liquid water flux is assumed to follow the Darcy-Buckingham law (Sposito, 1978):

$$ F_w = -K^{\text{sat}} k_r (\nabla \psi + \cos(\gamma)) \quad (7) $$

where $K^{\text{sat}}$ is the saturated hydraulic conductivity (in m s⁻¹), which does not depend on the LWC but only the properties of the ice matrix, $k_r$ the relative hydraulic conductivity, that depends on the LWC, $\psi$ the liquid matric potential (in m and that can be negative due to capillary forces), and $\gamma$ the slope angle. In the same way that the snow energy and temperature are related through Eq 3, the LWC and the matric potential are related through the use of a water retention curve (WRC). This WRC can be expressed with the use of, for instance, the van Genuchten or Brook models (Brooks and Corey, 1964; van Genuchten, 1980; Lenhard et al., 1989). In this article, and consistently with previous snow-related work (Wever et al., 2014; D'Amboise et al., 2017), we use a van Genuchten model for the WRC in snow:

$$ \theta_w = \begin{cases} \theta_r + (\theta_s - \theta_r)(1 + |\alpha\psi|^n)^{\frac{1-n}{n}} & \text{if } \psi < 0 \\ \theta_s & \text{otherwise} \end{cases} \quad (8) $$

where $n$ and $\alpha$ (in m⁻¹) are parameters characterizing the shape of the WRC and which depends on the snow microstructure (Yamaguchi et al., 2012), and $\theta_r$ and $\theta_s$ are referred to as the residual and saturation LWC, respectively. Physically, $\theta_s$ corresponds to the maximum amount of liquid water a snow sample can hold (thus typically assumed to correspond to the total porosity), and $\theta_r$ to the minimum amount of water that can be reached with draining through capillary and gravity forces.

To solve Richards' equation, it is as well necessary to provide an expression for the dependence of the relative hydraulic conductivity $k_r$ to the liquid water content. Following Mualem (1976), in complement of the van Genuchten WRC, the hydraulic conductivity is usually taken as

$$ k_r = \sqrt{S} \left( \frac{\int_0^S \frac{dS}{\psi}}{\int_0^1 \frac{dS}{\psi}} \right)^2 = \sqrt{S} \left( 1 - \left( 1 - S^{\frac{n}{n-1}} \right)^{\frac{n-1}{n}} \right)^2 \quad (9) $$
<figure><img src="image_10.png" /></figure>

where $S = \frac{\theta_w - \theta_r}{\theta_s - \theta_r}$ is referred to as the saturation degree. It can be noted that the relative hydraulic conductivity is an increasing function of $\theta_w$, that vanishes at the residual LWC $\theta_r$ and reaches unity value at water saturation.

However, Richard's equation presents specific problems in the case of water-saturated and fully-dry materials. In the case of a water-saturated materials, the liquid water content reaches a maximum value, while the matric potential can keep increasing. In other words, in the saturated case, the equation becomes degenerate and can no longer be expressed in terms of liquid water content. This is a classical difficulty associated with Richard's equation, usually circumvented by using $\psi$ as the primary variable to describe the material (Celia et al., 1990; Farthing and Ogden, 2017). It can also be treated using a variable switch technique, effectively changing from $\theta_w$ to $\psi$ as a primary variable depending on the degree of saturation of the material (Diersch and Perrochet, 1999; Krabbenhøft, 2007; Bassetto et al., 2022). This latter technique is explored below in the article. The problem of a fully-dry material is more specific to snow, and requires a dedicated treatment.

### 2.2.1 The problem of Richards's equation for dry snow

Indeed, in their usual forms, the WRCs used as material laws in Richards' equation do not apply in dry snow. In Eq. 8, the matric potential $\psi$ diverges when the LWC approaches its residual value $\theta_r$. In material solely driven by capillary and gravity flow, this implies that $\theta_w$ cannot drop below this residual value, i.e. that liquid water is always present. To the contrary, snow can become fully dry when energy is removed and the residual liquid water freezes. In such a case, the LWC can fall below the residual liquid water content and even vanish. Then, the matric potential $\psi$ becomes undefined. In snowpack models, this issue of dry snow has been circumvented in previous implementations of Richards' equation by keeping a small liquid water content, even in the case of snow below its fusion temperature (Wever et al., 2014; D'Amboise et al., 2017). This however hinders a consistent treatment of coupled heat and liquid water budget, as snow below the fusion is considered dry while solving the heat budget but wet while solving the liquid water budget. Furthermore, this technique requires the residual LWC to be artificially modified, in order to remain strictly below the liquid water content at all times. Thus, snow containing very little liquid water will tend to percolate, even though the unmodified WRC would rather imply that the liquid water should be held still under capillary forces.

This issue of the disappearance of a phase and of the divergence of the capillary forces is not only encountered in snow modeling. It is for instance present in underground nuclear-waster storage, where phases can appear and disappear (e.g. Bourgeat et al., 2009). To the best of our knowledge, two types of solutions to this problem have been proposed. First, the gradient of matric potential can be rewritten in terms of the gradient of liquid water content using the chain-rule (similarly to Bourgeat et al., 2009). The liquid water flux is thus taken as $-K^{\text{sat}}k_r(\nabla\psi + \cos(\gamma)) = -K^{\text{sat}}k_r(\partial_{\theta_w}\psi\nabla\theta_w + \cos(\gamma))$. While the derivative $\partial_{\theta_w}\psi$ diverges near the residual point, the product $-K^{\text{sat}}k_r\partial_{\theta_w}\psi$ vanishes near the residual point (as $k_r$ rapidly goes to zero). This product can thus be extended at and below the residual point with a value of zero, allowing the computation of the liquid water flux to be defined at and below the residual point. However, using this technique for the computation of the liquid water
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

flux presents one main issue. The use of the liquid water content gradient as the driving force for the liquid water flux implies that the state of equilibrium is characterized by a uniform liquid water content rather than uniform water (matric + gravitational) potential. While this is not a problem in the case of a homogeneous medium, it will negatively impact the representation of capillary barriers in a stratified medium, unless a specific treatment is implemented at the interfaces between the different strata (Amaziane et al., 2012). This is a major issue for stratified snowpack, as capillary barriers are crucial features that need to be captured (Wever et al., 2016). The second technique proposed in the literature to circumvent the singularity of the WRC near the residual point is to regularize it, so that it does not diverge (Beaude et al., 2019; Moure et al., 2023). The matric potential simply assumes a finite value near and below the residual point.

For our article, we rely on the second technique, namely the regularization of the WRC. Concretely, this is done by limiting the divergence of the retention curve up to a critical liquid water content $\theta_{lim}$ and thus a corresponding limit matric potential $\psi_{lim}$. Below that point, the liquid water content is allowed to further decrease (down to 0 in the case of fully dry and cold snow) while the matric potential remains at the value $\psi_{lim}$. This was done by choosing a unique limit saturation degree $S_{lim} = \frac{\theta_{lim}-\theta_r}{\theta_s-\theta_r}$, applied to all snow, below which the WRC reaches its plateau. Some regularized WRC are depicted in Fig. 1 for three different snow samples. The value of $S_{lim}$ needs to be taken small (typically $10^{-10}$), in order not to restrict too much the sharp increase of the matric potential at low LWCs. With this modification of the WRC, $\theta_{lim}$ assumes a role similar to that of the residual value: it corresponds to the point that cannot be further dried with water flow alone. To be consistent with this idea, we propose to also modify the expression of the relative hydraulic conductivity such that it vanishes at $\theta_{lim}$ rather than $\theta_r$. Specifically, we take

$$k_r = \begin{cases} \sqrt{S} \left( \frac{\int_{S_{lim}}^{S} \frac{dS}{\psi}}{\int_{S_{lim}}^{S_s} \frac{dS}{\psi}} \right)^2 = \sqrt{S} \left( 1 - \left( \frac{1-S^{\frac{n-1}{n-1}}}{1-S^{\frac{n-1}{n-1}}} \right)^{\frac{n-1}{n}} \right)^2 & \text{if } S > S_{lim} \\ 0 & \text{otherwise} \end{cases} \quad (10)$$

which vanishes when $\theta_w = \theta_{lim}$. As $\theta_{lim}$ corresponds to the point where liquid cannot flow anymore, we refer to as it as the retention point, to distinguish it from the residual point $\theta_r$. Note that while we have constructed the regularization as a plateau, other choices could be made. For instance, the WRC could be regularized through a linear function (Moure et al., 2023). With the help of regularized WRCs, Richard's equation now applies to describe both dry and wet snow, and can naturally handle regime changes between them.

## 2.3 Ice budget

Besides the energy and liquid water budgets, the process of water melting and freezing needs to be accounted for in the ice budget as well, and this in order to properly close the mass budget of the snowpack. There are two potential ways to account for melting refreezing event in the ice budget while respecting mass balance: it can be either viewed as a decrease or increase in density at constant volume, or as a loss or gain of volume at constant density. To the best of our knowledge, the choice in snowpack models is usually to treat refreezing as a density increase and the melting as a thinning of the snowpack, with no
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 1.** Examples of the regularized water retention curves used in this work, for three different snow density and surface specific area (SSA).</figcaption></figure>

direct impact on the density (Bartelt and Lehning, 2002; Vionnet et al., 2012). The justification behind this choice follows the observation that melting snow is usually of a high density, and thus that the melting process should not act as a de-densification mechanism. Yet, as phase changes occur directly within the snow microstructure, at the surface of the porous ice matrix, we rather believe that both melting and refreezing rather impact the snow density, without direct impact on the thickness. With this idea, the relatively high-density nature of melting snow could be attributed to the low viscosity of wet snow (Vionnet et al., 2012) that easily compacts under mechanical stress. Therefore, for the rest of the paper we assume that both type of phase changes affect the snow density, and not directly the snow layer thicknesses.

210

215

The ice mass budget within the snowpack is thus given by

$$ \rho_i \partial_t \phi_i = M_{\text{freeze}} \qquad (11) $$

where $\phi_i$ is the ice volume fraction and $\rho_i$ the ice density (in kg m$^{-3}$).

## 2.4 Resulting equations

Replacing Eq 5 by its sum with $\rho_w L_{\text{fus}}$ Eq 6, we can rewrite the system of equation as:
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

$$220 \quad \left\{ \begin{array}{l} \partial_t h - \nabla \cdot (\lambda \nabla T + \rho_w L_{\text{fus}} \nabla \cdot (K^{\text{sat}} k_r (\nabla \psi + \cos(\gamma)))) - Q_{\text{abs}} = 0 \\ \rho_w \partial_t \theta_w - \rho_w \nabla \cdot (K^{\text{sat}} k_r (\nabla \psi + \cos(\gamma))) + M_{\text{freeze}} = 0 \\ \partial_t m_i - M_{\text{freeze}} = 0 \end{array} \right. \quad (12)$$

where $h = h_s + \rho_w L_{\text{fus}} \theta_w = c_p (T - T_0) + \rho_w L_{\text{fus}} \theta_w$ is the sum of energy carried by the snow and the liquid water within. The first equation corresponds to the total energy budget. It includes the processes of heat conduction and advection of energy with liquid water as fluxes, as well as the volumetric source of energy. The release/absorption of latent heat is not included as a source term in this equation. Indeed, such phase change does not change the total energy budget, as latent heat is converted into snow sensible heat, both included in $h$.

In its current form, this system of equation includes 6 unknown, namely $h$, $T$, $\psi$, $\theta_w$, $M_{\text{freeze}}$ and $\phi_i$ ($K^{\text{sat}}$, $k_r$, and $\lambda$ being material properties that can be deduced from these unknowns). To this system of equation can be added 2 constitutive laws: Eq. 3 relating the snow temperature to its sensible heat and Eq. 8 relating the matric potential to the liquid water content.

The system of equations is thus not yet closed, as it remains to express how $M_{\text{freeze}}$ relates to the other variables. To the best of our knowledge, there is no broadly accepted theoretical or experimental work providing the rate of melting or freezing of water in snow (Moure et al., 2023). As most snowpack models (Jordan, 1991; Bartelt and Lehning, 2002; Vionnet et al., 2012; Sauter et al., 2020), we assume that liquid water and the snow are in local thermodynamic equilibrium. Specifically, it is assumed that snow below the fusion temperature and liquid water cannot coexist. This assumption is consistent with the slow matric flow considered in this article. In this framework, $h$ and $\theta_w$ become directly related: in the case where $h$ is below the fusion value, $\theta_w$ vanishes; in the case where $h$ is above the fusion value, $\theta_w$ is proportional to the excess of energy. In this case, the liquid water content and the temperature can be expressed as a function of the total energy:

$$\theta_w = \begin{cases} 0 & \text{if } h < 0 \\ \frac{h}{\rho_w L_{\text{fus}}} & \text{otherwise} \end{cases} \quad (13)$$

and

$$240 \quad T = \begin{cases} T_0 + \frac{h}{c_p} & \text{if } h < 0 \\ T_0 & \text{otherwise} \end{cases} \quad (14)$$

The rate of freezing/melting can then be derived, as a diagnostic, from the closure of Richard's equation. It physically corresponds to the amount of frozen and melted water required to maintain the local thermodynamic equilibrium. It can also be eliminated from the system of equations by summing Eqs 6 and 11, simply resulting in:
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

$$245 \quad \left\{ \begin{array}{l} \partial_t h - \nabla \cdot (\lambda \nabla T + \rho_w L_{\text{fus}} \nabla \cdot (K^{\text{sat}} k_r (\nabla \psi + \cos(\gamma)))) - Q_{\text{abs}} = 0 \\ \rho_w \partial_t \theta_w + \rho_i \partial_t \phi_i - \rho_w \nabla \cdot (K^{\text{sat}} k_r (\nabla \psi + \cos(\gamma))) = 0 \end{array} \right. \quad (15)$$

with the second equation now corresponding to the total mass budget.

Finally, we note that the recent article of Moure et al. (2023) proposes to relax the assumption of local thermal equilibrium and to introduce a finite rate of phase change, derived from the upscaling of the Frenkel-Wilson equation. This implies that the ice and liquid water temperatures in a snow sample are in general different and can be above or below the fusion point. This modeling framework, composed of four partial differential equations, is briefly presented in Appendix A. However, as observed in the Appendix, the timescale of relaxation towards local thermodynamical equilibrium appears to be much smaller than the timescale of matric water movement and heat diffusion, supporting the standard assumption of local equilibrium in snowpack models.

## 2.5 Choosing primary variables in the presence of regime changes

A difficulty with the system presented in Eq 15 (complemented with its material laws) is the presence of 4 distinct regimes, where parts of the equations behaves differently and sometimes degenerate. The first regime corresponds to dry snow, when the temperature is below the freezing point. Here, the thermodynamical state of snow can be characterized by either its temperature or its energy content (the two being related by the thermal capacity). To the contrary, the LWC and the matric potential cannot be used to describe the state of snow, as they assume constant, degenerate, values. In the second regime, the snow is above its fusion point, but the LWC remains below its retention value. Here, the snow can be characterized by its energy content or its LWC, but not by its temperature nor its matric potential. The third regime corresponds to the classical unsaturated Richards' case, where the LWC is above the retention point. Here, the snow can be characterized by its energy content, its LWC, or its matric potential, but not by the temperature. The fourth regime corresponds to the water-saturated regime. Here, the snow can only be characterized by its matric potential. This regime corresponds to the saturated regime of Richards' equation. A summary of the different regimes and the variables that can be used to characterize them is given in the Table 1. As seen, there is no natural variable that can be used to describe the thermodynamical state of snow over the 4 different regimes. While the energy content could be used to characterize snow in regimes 1, 2, and 3, it degenerates in regime 4. Similarly, while the matric potential is a good candidate to describe regimes 3 and 4 (as usually done when numerically solving Richards' equation; Celia et al., 1990; Farthing and Ogden, 2017), it fails in regime 1 and 2.

To circumvent this issue, we rely on a variable switch technique by parametrization (Brenner and Cancès, 2017), through the use of a fictitious variable meant to behave as the energy content in the first three regimes and as the water matric potential
<figure><img src="image_19.png" /><figcaption>**Table 1.** Summary of the various regimes encountered in snow. For each regime, the relevant variables that can be used to characterize snow are given, as well as the degenerated variables that cannot be used to characterize snow. In all cases, the ice volume fraction $\phi_i$ is also needed as a second variable.</figcaption></figure>

<table><thead><tr><th></th><th>Dry snow</th><th>Wet non-flowing snow</th><th>Wet unsaturated snow</th><th>Wet saturated snow</th></tr></thead><tbody><tr><th>Relevant variables</th><td>h, T</td><td>h, θ<sub>w</sub></td><td>h, θ<sub>w</sub>, ψ</td><td>ψ</td></tr><tr><th>Degenerated variables</th><td>θ<sub>w</sub>, ψ</td><td>T, ψ</td><td>T</td><td>T, h, θ<sub>w</sub></td></tr></tbody></table>

275 in the last one (Bassetto et al., 2021, 2022). For this, we introduce a new variable $\chi$ (expressed in J m$^{-3}$) used to parameterize the {$h, \psi$} graph. Specifically, we choose $\chi$ such that:

$$h(\chi, \phi_i) = \begin{cases} \chi & \text{if } \chi < \theta_s(\phi_i)\rho_w L_{\text{fus}} \\ \theta_s(\phi_i)\rho_w L_{\text{fus}} & \text{otherwise} \end{cases} \qquad (16)$$

$$\psi(\chi, \phi_i) = \begin{cases} \psi(\theta_w(\chi, \phi_i)) & \text{if } \chi < \theta_s(\phi_i)\rho_w L_{\text{fus}} \\ \frac{\chi - \theta_s(\phi_i)\rho_w L_{\text{fus}}}{\beta} & \text{otherwise} \end{cases} \qquad (17)$$

where $\beta$ is an arbitrary constant, introduced to respect the unit homogeneity of Eq 17, and where the dependence of $\phi_i$ as also been made explicit. Thanks to this variable $\chi$, we are now able to characterize the energy content and matric potential of snow in all regimes, and from it to derive the values of all other relevant quantities for snow ($\theta_w, T, k_r$, etc).

Once parameterized with $\chi$ our system of equation becomes:

$$\begin{cases} \partial_t h(\chi, \phi_i) - \nabla \cdot (\lambda \nabla T(\chi, \phi_i) + \rho_w L_{\text{fus}} \nabla \cdot (K^{\text{sat}} k_r (\nabla \psi(\chi, \phi_i) + \cos(\gamma))) - Q_{\text{abs}} = 0 \\ \rho_w \partial_t \theta_w(\chi, \phi_i) + \rho_i \partial_t \phi_i - \rho_w \nabla \cdot (K^{\text{sat}} k_r (\nabla \psi(\chi, \phi_i) + \cos(\gamma))) = 0 \end{cases} \qquad (18)$$

285 that can be solved searching for $\chi$ and $\phi_i$. We have explicitly shown the dependency of $h, T, \psi$, and $\theta_w$ on both $\chi$ and $\phi_i$, but one should keep in mind that the effective properties $\lambda, K^{\text{sat}}$, and $k_r$ also depends on the value of $\chi$ and $\phi_i$. With this set of equations, we now have a consistent description of the energy and mass budgets in snowpacks, that naturally applies to both dry and wet snow.

## 3 Numerical implementations

290 The goal of this section is to present and explore different numerical implementations of the energy and mass budgets in snow, in order to investigate the impacts of the implementation on the results and robustness of the models. We decided to focus this
<figure><img src="image_20.png" /></figure>

comparison on (i) the use of a regularized WRC with variable switching for matric flow, and (ii) the loosely or tightly-coupled nature of the implementation, where processes can either be solved sequentially (i.e. with operator-splitting) or concomitantly (i.e. with tight-coupling) (Steefel and McQuarrie, 1996; Keyes et al., 2013). By varying the regularization of the WRC and the degree of operator-splitting, we have implemented five snowpack toy-models. Below, we start by briefly presenting the common features shared by all implementations, and then discuss the concrete differences between them. These differences between the numerical implementations are summarized in Table 2.

## 3.1 A common core between all models

In order to be directly comparable, the models share portions of their numerical implementations. This includes compaction under mechanical stress (assuming a linear viscosity between stress and deformation and which is solved after the thermodynamics), the constitutive laws defining the material properties of snow, the spatial and temporal discretization, and the method for solving the resulting non-linear systems of equations.

### 3.1.1 Constitutive laws

To be fully closed, the equations of energy and mass conservation need to be complemented with constitutive laws, prescribing the material properties of the snow material. We assume the thermal conductivity $\lambda$ of snow to follow Calonne et al. (2011) and the saturated hydraulic conductivity $K^{sat}$ to follow Calonne et al. (2012). The WRC (Eq. 8) and the relative hydraulic conductivity (Eq. 10) require defining the parameter $\alpha$ and $n$. For this, we follow Yamaguchi et al. (2012). The residual LWC $\theta_r$ is taken as 0.02 (Yamaguchi et al., 2010; Wever et al., 2014; D'Amboise et al., 2017) and the saturated LWC $\theta_s$ as $1 - \phi_i$ (i.e. the total porosity). Finally, for the compaction of the snowpack under mechanical settling, we use the linear compactive viscosity of Vionnet et al. (2012), including its dependence to the LWC. The precise expressions for the constitutive laws are given in Appendix B.

### 3.1.2 Finite Volume discretization

To be numerically solved, the system of Eqs. 18 needs to be discretized in time and space. For the time discretization, we use a simple backward Euler time-stepping. This choice is motivated by the overall stability of the method, mitigating the apparition of overshoots and oscillations in the case of large-time steps (Butcher, 2008). For the spatial discretization, we use a finite volume scheme. Briefly, this numerical method relies on dividing the snowpack into $N$ cells and performing energy and mass balances on each individual cell. In this framework, each cell is described by its average $\chi$ and $\phi_i$ values. The evolution of these average values is obtained by integrating the system of equations over the different cells, and applying the divergence theorem to transform the divergence operator into fluxes at the cells' boundaries.
<figure><img src="image_21.png" /></figure>

To be mathematically closed, the heat and liquid water fluxes at the interface between two adjacent cells need to be reconstructed from the cells' averages. For the heat fluxes, this is done by computing the temperature gradient from one cell center to the other, and defining an interfacial thermal conductivity. As classically done in finite volume schemes, this interfacial thermal conductivity $\lambda^{k+1/2}$ (between cells $k$ and $k+1$) is taken as some average of the thermal conductivities of the adjacent cells. In our case, we take it as the weighted harmonic average. This choice ensures that the heat flux vanishes when the thermal conductivity of one of the two cells vanishes (Kadioglu et al., 2008).

For the computation of the liquid water, we discretize as well the matric potential gradient using the average values in the cell and the center-to-center distance. For the computation of the interfacial hydraulic conductivity, we use the so-called upstream mobility formulation (Bassetto et al., 2021). In this framework, the hydraulic conductivity of the interface is split into the product of the saturated conductivity $K_{sat}^{k+1/2}$ and the relative conductivity $k_r^{k+1/2}$. As with the thermal conductivity, the interfacial saturated conductivity is taken as the harmonic average of the saturated conductivities of the neighboring cells. This ensures that the liquid water flux vanishes when on the cell is impermeable. The interfacial relative conductivity $k_r^{k+1/2}$ between cell $k$ and $k+1$ is taken as the upstream value, i.e.:

$$k_r^{k+1/2} = \begin{cases} k_r^k & \text{if } \psi^k > \psi^{k+1} + \cos(\gamma)d^{k+1/2} \\ k_r^{k+1} & \text{otherwise} \end{cases} \qquad (19)$$

where $d^{k+1/2}$ is the distance between the neighboring cell centers and $\gamma d^{k+1/2}$ thus is the vertical distance between the cell centers. With this upstream choice, the numerical scheme becomes monotonic (in the sense of Eq. (3.18) of Bassetto et al., 2021). This property is notably beneficial for the convergence of the Newton (or other iterative) method, when solving the resulting non-linear system of equations.

### 3.1.3 Truncated Newton method with adaptive time-stepping

The systems of equations to be numerically solved in the different models are non-linear and require an iterative scheme. For that, we rely on a Newton method, with a stopping criterion when the iterated estimations is close enough to the solution. In models 4 and 5, this criterion is complemented with a criterion on mass conservation, as these numerical scheme are not naturally mass-conservative. While the Newton algorithm provides a relatively fast convergence when the starting estimation is close to the solution, it is not a globally convergent algorithm. In other words, it is possible that the algorithm diverges or get stuck in cycles. To improve its convergence performance, we implemented two strategies. First, as usually done with Richards' equation, we use an adaptive time-step. In the case where the algorithm fails to converge after a given number of iterations, the algorithm is rewind to the start of the time-step and its value halved. By default, we set the maximum number of iterations to 25. In the rest of the article, we thus make the difference between the so-called "default timestep", corresponding to the timestep that is ideally used if no convergence problems are encountered, and the "adaptive timestep", that is effectively used to solve the nonlinear equations involving matric flow. Secondly, we use the so-called truncation method when a transition from one regime to another occurs (Bassetto et al., 2020). Indeed, the transitions between regimes are corner points, characterized
<figure><img src="image_22.png" /></figure>

by discontinuities in the derivatives. Using the derivative computed on one side of a corner point to derive the evolution of the estimate in the other side is therefore problematic and can lead to overshoots, impeding the convergence towards the solution. To avoid this problem, each time a transition from one regime to another occurs in a cell of the snowpack, $\chi$ is set back in the vicinity of the corner point. In practice, $\chi$ is placed at a small $\epsilon$ value before or after the corner point, in order to fall in the good regime. We set this value of $\epsilon$ to be $10^{-5}$ J m$^{-3}$. Finally, we note that for solving Richards' equation the SNOWPACK (Wever et al., 2014) and Crocus (D'Amboise et al., 2017) models decided to follow a modified Picard method (Celia et al., 1990), where the contribution of the change in hydraulic conductivity in the Jacobian is not taken into account. As a test, we also run simulations using the modified Picard rather than the Newton method for models 4 and 5. This did not modify the results discussed later.

## 3.2 Differences between the five different numerical implementations

## 3.3 Model 1: Fully coupled with a regularized WRC

For the first model, we use the regularized WRC and solve the energy and mass budgets in a tightly-coupled manner, that is to say that the system of Eqs 18 is solved concomitantly for both $\chi$ and $\phi_i$. In this version, there is therefore no degree of operator-splitting when solving the thermodynamical state of the snowpack.

## 3.4 Model 2: Partially coupled with a regularized WRC

While the numerical implementation of model 1 appears to be the most consistent, it results in a large non-linear system that can be numerically expensive to solve. As a compromise, some degree of operator-splitting (i.e. sequentiality) between processes can be introduced. For model 2, we thus used the regularized WRC with a decoupling of the computation of the energy and liquid water budgets from that of the ice budget. The motivation behind it is that the timescale for the evolution of the density is usually longer (of the order of the day) than that of the energy and liquid water evolution (of the order of the hour or less). Concretely, in Eq. 18 we thus first solve the evolution of $\chi$ (assuming a fixed $\phi_i$). Then, we find the evolution of $\phi_i$ by closing the mass budget.

## 3.5 Model 3: Full operator-splitting with a regularized WRC

Keeping with the idea of using operator-splitting to reduce the numerical cost and complexity of the model, the third implementation is based on a decoupling of the energy, liquid water, and ice budgets solving. Specifically, we first evolve $\chi$ under the process of heat conduction and phase changes only. We then, evolve $\chi$ again, this time under the process of matric liquid water flow. Finally, $\phi_i$ is updated by re-applying phase transition and closing the mass budget.

## 3.6 Model 4 and 5: Full operator-splitting without a regularized WRC

The last two models are meant to emulate the techniques used so far in snow models to handle the degeneracy of the retention curve in dry snow. Specifically, model 4 is based on the SNOWPACK implementation (Wever et al., 2014) and model 5 on
<figure><img src="image_23.png" /><figcaption>**Table 2.** Description of different implemented toy-models based on (i) the regularization of the WRC, (ii) whether the primary variable is a switch (behaving as *h* in unsaturated snow and *ψ* in saturated snow) or *ψ* for solving the matric flow process, and (iii) the degree of couplings between thermodynamics processes.</figcaption></figure>

<table><caption>**Table 3.** Description of three different test cases considered in this article, with their durations and external forcings.</caption><thead><tr><th></th><th>Model 1</th><th>Model 2</th><th>Model 3</th><th>Model 4 and 5</th></tr></thead><tbody><tr><th>WRC</th><td>Regul.</td><td>Regul.</td><td>Regul.</td><td>Non Regul.</td></tr><tr><th>Primary variable for liquid flow</th><td>Switch</td><td>Switch</td><td>Switch</td><td>ψ</td></tr><tr><th>Degree of coupling</th><td>- Conduction, phase change, and liquid flow coupled.</td><td>- Conduction and phase change coupled.<br/>- Liquid flow decoupled</td><td>- Conduction, phase change, and liquid flow decoupled</td><td>- Conduction, phase change, and liquid flow decoupled</td></tr></tbody></table>

<table><thead><tr><th></th><th>Test case 1<br/>Melting snowpack</th><th>Test case 2<br/>Low-intensity rain</th><th>Test case 3<br/>High-intensity rain</th></tr></thead><tbody><tr><th>Duration</th><td>6 days</td><td>1 day</td><td>2 day</td></tr><tr><th>Shortwave</th><td>840 W m<sup>-2</sup> peak</td><td>420 W m<sup>-2</sup> peak</td><td>420 W m<sup>-2</sup> peak</td></tr><tr><th>Longwave and Turbulent fluxes</th><td>200 to 310 W m<sup>-2</sup></td><td>310 W m<sup>-2</sup></td><td>310 W m<sup>-2</sup></td></tr><tr><th>Rainfall rate</th><td>No rain</td><td>constant 1 mm h<sup>-1</sup></td><td>15 mm h<sup>-1</sup> for 4 h</td></tr></tbody></table>

the Crocus implementation (D’Amboise et al., 2017). The specificities of the implementation are given in Appendix C, but we recall that both techniques are based on the idea of introducing a small amount of liquid water in cold snow, shifting the residual liquid water content value below it, and using *ψ* as the primary variable. Moreover, these models rely on a decoupled resolution of the thermodynamical processes.

385

390

## 4 Numerical tests

The goal of this section is to build simple synthetic examples, meant to represent different situations in which liquid water percolation might occur in snowpacks, and to investigate the differences in behavior between the various implementations.
<figure><img src="image_24.png" /></figure>

## 4.1 Test case 1: Melting snowpack

The first test case is meant to represent the situation of melting snowpack, releasing liquid water near the surface that then percolates downward. For the initialization of the snowpack, we rely on a simulation performed with the Crocus snowpack (Vionnet et al., 2012). This yields a realistic stratigraphy, with thinner cells near the surface in order to better capture surface effects. The initial state of the simulation was chosen to obtain a dry snowpack near its peak snow water equivalent. This defines the ice volume fraction, temperature, and specific surface area at the start of our simulations. Note that the specific surface area is needed as it plays a role in the hydraulic properties of the snowpack (Calonne et al., 2012; Yamaguchi et al., 2012). However, as we did not implement metamorphism laws in our toy model, the specific surface area is simply kept constant during the simulation.

For the top boundary conditions, we impose an incoming diurnal energy flux at the surface of the snowpack (oscillating between 310 W m⁻² at noon and 200 W m⁻² at night), encapsulating the effect of incoming longwave radiations and turbulent heat fluxes. We also impose diurnal shortwave radiations (peaking at 840 W m⁻² at noon and vanishing at night), that penetrate within the snowpack (with an albedo of 0.7 and an e-folding depth of 5 cm; Libois et al., 2013). No liquid, nor solid, precipitation is considered in this test case. As the bottom boundary condition, we impose a constant heat flux of 10 W m⁻², emulating the heat flux received from a warm ground below, and a free-drainage condition for the liquid water flux. For each model, simulations are run for 6 days, so that water percolates to the bottom of snowpack. Default timesteps vary between 112 and 7200 s (we recall that because of the adaptive timestep strategy, the actual timestep used for some process might be shorter than the prescribed default value).

## 4.2 Test case 2: Water infiltration under low-intensity rain

This second test case is meant to represent the slow infiltration of rain under a long but low-intensity event. The initialization is the same as in case 1, based on the output of the same Crocus simulation. For the top boundary condition, we impose a constant surface energy flux of 310 W m⁻², a shortwave radiation flux peaking at 420 W m⁻² at noon (lower than in the first test case to represent cloudiness), and a constant rain flux of 1 mm h⁻¹ lasting the whole simulation. For the bottom boundary condition, we use the same conditions as in the test case 1 (constant heat flux and free-drainage). Simulations were performed with each model for 1 day, in order for shallow infiltration to occur. Default timesteps range from 112 to 7200 s.

## 4.3 Test case 3: Water infiltration under high-intensity rain

The third test case represents the case of a short but high-intensity rain event. The initialization and boundary conditions are taken as in the test cases 2, except for the incoming rain flux that is now of 15 mm h⁻¹ and lasts for 4 h. This results in an abrupt input of liquid water in the snowpack, rapidly percolating downward. Simulations were run for about 2 days, long enough for deep percolation to occur. Default timesteps range from 112 to 7200 s.
<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 2.** Results of the reference simulations for the three investigated test cases: a) surface melting with a deep infiltration of surface melt, b) a low-intensity rain with limited infiltration, c) a high-intensity rain event with a fast and deep water infiltration.</figcaption></figure>

For all three test cases, reference simulations were performed with the fully coupled model (model 1) and a timestep of 1 s.
The resulting LWC profiles over time are displayed in Fig. 2.

# 5 Results and discussion

## 5.1 Robustness and numerical cost

The first difference in behavior between the different implementations is their respective numerical cost, i.e. the computation time required to perform a given simulation. In the presence of an adaptive timestep, required to accommodate the strong non-linearity of the liquid matric flow and of the WRC, we observed that the numerical cost of the models is, at first order, controlled by the ability of the models to stick to large adaptive timesteps. Figure 3 displays the adaptive timesteps required by the models in the three test cases presented above. While the models using a regularized retention curve with a variable switch (models 1, 2, and 3) are able to maintain a timestep of around 1 h for most of the simulations, models 4 and 5 requires the adaptive timestep to regularly drop, sometimes below 10 s. In the end, models 4 and 5 require about one order of magnitude more steps to perform the same simulation, resulting in significantly increases numerical cost and computation time. For instance, in the case of the test case 3 with a default timestep of 3600 s, models 4 and 5 require a computation time of about 500 s, while models 1, 2, and 3 require about 130 s, 25 s, and 45 s, respectively. Note that these numbers should only be interpreted qualitatively to assess numerical complexity, as the precise computation time is also influenced by factors that we did not control for (such as the level of code optimization or memory-caching). Therefore, the use of a regularized WRC with a variable switch appears as a very favorable practice to increase the robustness of the numerical scheme and to decrease its numerical cost by accommodating larger timestesp. Our understanding for this behavior is that the combination of regularized WRC and not using the matric potential as the primary variable near the retention point hinders the presence of abrupt, and nonlinear, variations of the numerical solution, thus making the problem easier to solve. To the contrary, using the matric potential as the primary variable is known to be inefficient for low LWC materials and to require timesteps that are much smaller than the
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 3.** Adaptive timesteps required to simulate matric flow in each test case (columns) and for each model implementation (rows). The default value for the timestep is set to 3600 s. In each plot, the marker marks the total number of timesteps required and the average timestep size.</figcaption></figure>

actual timescale of liquid water transport (Forsyth et al., 1995; Sadegh Zadeh, 2011).

Focusing on models 1, 2, and 3 in Fig. 3, it appears that these three implementations require very similar timesteps, apart from a notable drop for model 3 in the first and second test case. Their differences in numerical cost should be controlled by the complexity of the assembly and solving of the matrix-form of the equations. According to our numerical test cases, model 2 and 3 have about the same numerical cost, while our implementation of model 1 requires a longer computational time. This is to be expected, as model 1 requires the largest matrix to be assembled and solved ($2N \times 2N$), while the other model handle smaller matrices ($N \times N$). However, we also note that a large portion of the up-cost of model 1 comes from the computation of derivatives with respect to $\phi_i$ for the assembly of the Jacobian. In our implementation, computing these derivatives comes with a high cost, as we rely on generic automatic differentiation without any attempt at optimization. If these derivatives with respect to $\phi_i$ are ignored (thus transforming the scheme from a true Newton to a form of modified Picard method), the numerical solving of the coupled system of Eqs. 18 can be efficiently cheapen. For instance, in the aforementioned test case 3 with a default timestep of 3600 s, the computational time can be lowered from 130 to about 50 s, i.e. in the same range as models 2 or 3.
<figure><img src="image_30.png" /></figure>

## 5.2 Error-level and timestep convergence

Besides their numerical cost, the different implementations also yield different levels of error when compared to the reference simulations. While the error level generally tends to decrease with shorter timesteps, not all models perform the same for decreasing timesteps. Here, we thus investigate the degree error of each model, for default timesteps ranging from 7200 to 112 s. For that, we compute the Root Mean Square Difference (RMSD) between the models' outputs and the reference simulations. Figure 4 displays the RMSDs in LWC and temperature for the three test cases as a function of the default timestep. Note that the RMSD is computed over the whole snowpack, while errors tend to be spatially located near the percolation front. Therefore, the RMSD value is smaller, typically of an order of magnitude, than the errors occurring at the percolation front. As expected, the overall trend is a decrease of the RMSDs with small timestep. However, while the unregularized models 4 and 5 perform similarly as the other models with large default timesteps, they do not show a clear decrease of error at smaller timesteps. Our understanding is that the strategy of modifying the residual LWC based at the start of each timestep impacts the physics and changes the solution to which the models converge with small timesteps. This results in a plateau or even an increase of the associated RMSDs in Fig. 4. Concerning the regularized models (1, 2, and 3), our results suggest that model 1 consistently yields the lowest RMSDs, for timesteps of 900 s and less. The comparison of the errors of models 2 and 3 depends on the specific test case and the variable of interest, with no clear advantage of one over the other. We note that under some circumstances (e.g. test case 2 and timestep of 900 s), one model can appear better in a metric (e.g. LWC RMSD) while it is outperformed in another metric (e.g. temperature RMSD).

## 5.3 Is full-coupling beneficial ?

As seen above, the use of a tightly-coupled resolution in numerical models usually provides a slightly better precision and a better robustness (Keyes et al., 2013; Brondex et al., 2023), at the expense of a more complex numerical resolution. Indeed, by ensuring physical consistency between the different variables, tightly-coupled resolutions are known to prevent overshooting, which could develop into numerical oscillations or even divergence (Brondex et al., 2023). One of the motivation behind the investigation of the impact of operator-splitting on the numerical resolution of the energy and mass budgets with liquid water matric flow was thus to examine whether a tightly-coupled resolution was warranted in this case. To our surprise, the problem of heat conduction, liquid matric flow, and phase changes, appears to be quite robust to operator splitting. Indeed, even in the case of a quite large timestep of 4 hr, we did not observe any clear sign of numerical instabilities or very large overshoots in the models using operator-splitting. Our understanding is that this overall insensibility to operator-splitting comes from the nature of the involved processes, which do not play in the same regimes. Indeed, the process of heat conduction only occurs in dry snow, the process of matric flow in wet snow, and the process of phase change at the wet-dry interface. The possibility of interaction between processes thus remains limited and the large physical inconsistencies, that can be problematic with operator-splitting, did not emerge in most situations. Therefore, the standard practice of operator-splitting seems to be well-justified, at least when representing the interaction of matric flow with heat conduction and phase changes in snowpack models.
<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /><figcaption>**Figure 4.** Root Mean Squared Difference in LWC and temperature between the models and a reference simulation as a function of the default timestep size.</figcaption></figure>

## 5.4 Inclusion of other physical mechanisms

In this article, we considered the coupled mechanisms of heat conduction, matric liquid water percolation, and liquid-solid phase changes. While we focused on these mechanisms, as they can be found in all snowpack models (with more or less degrees of complexity), others thermodynamical mechanisms could be included as well. This is for instance the case of water vapor diffusion and vapor phase changes. To the best of our understanding, vapor-related mechanisms could be easily added to the systems of the thermodynamical equations governing snowpacks. For that, one would need to (i) introduce a new partial differential equation governing the mass balance of water vapor, (ii) introduce a water vapor component in the energy balance equation, and (iii) close the system by prescribing the phase change term for the gas phase (or assuming water vapor thermodynamical equilibrium with the other phases) and parametrizing the diffusivity of vapor. This could be done based on works of Calonne et al. (2014), Jafari et al. (2020), Brondex et al. (2023), and Bouvet et al. (2024).

Similarly, an important mechanism relating to the transport of liquid water in snowpacks is the process of preferential (i.e. non-matric) flow (Schneebeli, 1995; Yamaguchi et al., 2018). While the exact mechanisms, and therefore description, of preferential flow in snowpacks remain unclear at this point, the presence of fast flowing, and out-of-equilibrium with the rest of the snow layer, liquid water appears as a prerequisite for the formation of internal ice layers (Fierz et al., 2009). As percolation chimneys are likely to play a key role for preferential flow, the use of a so-called dual-porosity/dual-permeability model (e.g.
<figure><img src="image_33.png" /></figure>

Vogel et al., 2000; Lewandowska et al., 2004) has notably been proposed to take into account preferential flow in the SNOW-PACK model (Wever et al., 2016; Quéno et al., 2020). As preferential flow share many similarities with matric flow (being essentially a faster version), the governing equations are subjected to the same degenerate behaviors in the case of a dry or a fully-saturated medium. Therefore, the techniques explored in this article, namely the regularization of the WRC to handle dry media with the use of a switch variable, could also benefit the numerical implementation of preferential flow. Moreover, as preferential flow is a fast and out-of-equilibrium process, it can introduce stiffness in the equations governing snowpacks (Fazio, 2001), potentially resulting in overshoots and oscillations when using operator splitting with a large timestep (Brondex et al., 2023; Fourteau et al., 2024). As discussed above, the derivation of a fully-consistent system of equations, that applies from dry to water-saturated snow, enables a tightly-coupled numerical resolution, which could be key to handle stiff and fast systems of equations (Keyes et al., 2013).

# 6 Conclusions

This article focuses on the numerical implementation of matric water flow and its interaction with other physical processes in snowpack models using Richards' equation. While the use of Richards' equation can improve the representation of water flow, its complexity and numerical cost can sometimes hinder widespread adoption. To overcome this issue, we explored a new numerical implementation of Richard's equation. We started by recalling the governing equations of energy and mass conservation in snowpacks, that govern the evolution of the snow material both in dry and wet zones. However, for these equations to be directly applicable for both dry and wet snow, it is necessary to provide an expression of the liquid water flow that also applies in dry snow (where simply no flow is expected). For this, we regularized the water retention curve, in order to avoid the divergence of the matric potential as the snow material dries. Another issue with the concomitant and consistent representation of a whole snowpack is that dry and wet snow are not described by the same thermodynamical variables, as dry snow is described in terms of temperature/energy while wet snow is rather described in terms of liquid water content/matric potential. To seamlessly handle this change in the physical variables characterizing the snow material, we introduced a variable switch technique implemented thanks to the introduction of a fictitious variable building on the idea introduced in Brenner and Cancès (2017). This defines a single thermodynamical variable to describe all states of snow, from dry to water-saturated. Eventually, we thus obtained a consistent system of two equations (energy and mass conservation) governing the evolution of snowpacks and which allows a unified treatment of dry and wet snow.

To compare the behavior and performance of this new description, we have implemented five toy snowpack models. This includes three models, all using a regularized WRC with a variable switch but relying on various degrees of operator splitting during the resolution, as well as two models based on the treatment of Richards' equation in the SNOWPACK and Crocus snow models (i.e. using a non-regularized retention curve). Based on three test cases, representing various situations of snowpack humidification, we observe that the use of a regularized WRC with variable switch significantly increases the robustness of the numerical implementation, which can run with timesteps of 15 min and above, while the non-regularized models require smaller
<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

timesteps to handle liquid water percolation. On the other hand, the possibility of a fully tightly-coupled implementation seems to have a limited impact compared to an operator-splitting implementation, both in terms of stability and timestep sensitivity. While this article focused on matric liquid flow, we believe that the methodology put forward, namely the use of regularized laws, variable switch, and physically consistent description of the various thermodynamical states of snow, could also be useful for the description of preferential flow. Indeed, as preferential flow is usually modeled as a faster version of matric flox, the ability of regularized models to nonetheless handle relatively large timesteps could prove crucial to limit their numerical cost.

## Appendix A: Ice-Water thermodynamics equilibrium in snow

In this section, we briefly discuss the relaxation of the thermodynamics equilibrium assumption between the ice and water phases in snow. This appendix is largely based on the recent work of Moure et al. (2023). If we assume that the ice and liquid water are out-of-equilibrium, it is necessary to distinguish between the ice temperature and the liquid water temperature. It is thus no-longer possible to write a single energy-equation, as in Eq. 15. Similarly, we also need to distinguish between the solid and liquid mass budgets. The out-of-equilibrium ice and liquid water system is now described by $T_i$, $T_w$, $m_i$, and $m_w$ (the ice and liquid temperatures and the ice and liquid mass content, expressed in kg m⁻³) and the four equations provided by Moure et al. (2023):

$$ 
\left\{ 
\begin{aligned}
\partial_t(c_i m_i T_i) + \nabla \cdot J_{u,i} &= W_{\text{SSA}} \frac{\lambda_i}{r_i} (T_{\text{int}} - T_i) \\
\partial_t(c_w m_w T_w) + \nabla \cdot (J_w c_w T_w) + \nabla \cdot J_{u,w} &= W_{\text{SSA}} \frac{\lambda_w}{r_w} (T_{\text{int}} - T_w) \\
\partial_t m_i &= -\rho_i R_m W_{\text{SSA}} (T_{\text{int}} - T_0) \\
\partial_t m_w + \nabla \cdot J_w &= \rho_i R_m W_{\text{SSA}} (T_{\text{int}} - T_0)
\end{aligned}
\right. 
\qquad (A1) $$

where $J_{u,i}$, $J_{u,w}$, and $J_w$ are heat conduction fluxes in the ice, in the water, and the water flux, respectively, $R_m = \frac{c_w}{\beta L_{\text{fus}}}$ (with $\beta$ the kinetic attachment coefficient for ice growth from liquid water and $c_w$ the thermal capacity of liquid water), $W_{\text{SSA}}$ is the wet surface area, $T_{\text{int}}$ the temperature ice-liquid water interface, $\lambda_i$ and $\lambda_w$ are the ice and liquid water thermal conductivities, and $r_i$ and $r_w$ two microstructural parameters estimated to be 0.06 $d_i$ and 0.081 $d_i$, respectively (with $d_i$ the diameter of the grain in the snow microstructure; Moure et al., 2023). The temperature $T_{\text{int}}$ is given as a linear combination of the ice, water, and melting temperature:

$$ T_{\text{int}} = \frac{\frac{c_w \rho_w}{\beta} T_0 + \frac{\lambda_i}{r_i} T_i + \frac{\lambda_w}{r_w} T_w}{\frac{c_w \rho_w}{\beta} + \frac{\lambda_i}{r_i} + \frac{\lambda_w}{r_w}} = b_0 T_0 + b_i T_i + b_w T_w. \qquad (A2) $$

In this system of equations, the temperature of the ice-liquid interface governs the chemical equilibrium. Melting occurs if $T_{\text{int}} > T_0$, freezing occurs if $T_{\text{int}} < T_0$, and local equilibrium is achieved when $T_{\text{int}} = T_0$.
<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /><figcaption>**Figure A1.** Evolution of ice, liquid water, and interface temperatures towards equilibrium, starting from a non-equilibrium situation.</figcaption></figure>

A difficulty in this equation is relating $W_{SSA}$ to the actual specific surface area. Here, we follow Moure et al. (2023) and assume that $W_{SSA}$ scales linearly with the saturation degree $S$ (i.e. liquid water filling 50% of the porosity would wet 50% of the interfaces). However, we note that as the liquid water is the wetting phase in snow, it preferentially covers the ice surface, and the interface rapidly gets wet when liquid water is present. Therefore, the assumption that the wet SSA scales linearly with the saturation degree is likely an underestimation, which slows down the relaxation towards thermodynamics equilibrium in snow.

To estimate the timescale of relaxation, we have implemented a simple model simulating relaxation towards local equilibrium (ignoring spatial fluxes for simplicity in the system of Eqs. (A1)). Figure A1 shows the evolution of the temperatures of the ice, the liquid water, and the interface, with $T_i = 273.15 = T_0$ K, $T_w = 274.15$ K, $m_i = 275$ kg m⁻³, $m_w = 100$ kg m⁻³, and a SSA of $2$ m² kg⁻¹ as initial conditions. This illustrates that starting from a situation of non-equilibrium, thermodynamical equilibrium is achieved within about 20 ms. As this timescale of relaxation is much shorter than the timescale of evolution of the snowpack, the assumption of thermodynamical equilibrium between the ice and liquid water appears justified.

## Appendix B: Constitutive laws

This appendix presents the constitutive laws chosen for the models. For the thermal conductivity (expressed in W K⁻¹ m⁻¹) we have

$$ \lambda = 2.5 \times 10^{-6} \phi_i \rho_i - 1.23 \times 10^{-4} \phi_i \rho_i + 0.024 \qquad (B1) $$

and for the saturated hydraulic conductivity (expressed in m s⁻¹)

$$ K^{\text{sat}} = \frac{g \rho_w}{\mu_w} 3 r_{\text{es}}^2 \exp(-0.0130 \phi_i \rho_i) \qquad (B2) $$
<figure><img src="image_38.png" /></figure>

where $g = 9.81 \text{ m s}^{-2}$ is the gravity acceleration, $\rho_w$ is the density of liquid water and $\mu_w = 1.79 \times 10^{-3} \text{ Pa s}$ is the dynamic viscosity of water at $0^{\circ}\text{C}$, and $r_{es} = 3/(\text{SSA}\rho_i\phi_i)$ is the equivalent optical radius. For the $\alpha$ and $n$ parameters, defining the van Genuchten WRC, we have

$$585 \qquad \alpha = 4.4 \times 10^6 \left( \frac{\phi_i \rho_i}{2r_{es}} \right)^{-0.98} \qquad (B3)$$

and

$$n = 1 + 2.7 \times 10^{-3} \left( \frac{\phi_i \rho_i}{2r_{es}} \right)^{0.61}. \qquad (B4)$$

Finally, for the compactive viscosity $\eta$ (expressed in Pa s) we have

$$\eta = \frac{3.05 \times 10^{-7}}{1 + 60\theta_w} \frac{\phi_i \rho_i}{250} \exp(0.1(T_0 - T) + 0.023\phi_i \rho_i). \qquad (B5)$$

## Appendix C: Details on the implementation of models without WRC regularization

In this appendix, we briefly present the implementations of models 4 and 5, which are based on the treatment of Richards' equation in the SNOWPACK (Wever et al., 2014) and Crocus (D'Amboise et al., 2017) models, respectively. Note, that even though our implementations are based on SNOWPACK and Crocus, they cannot be considered as strict copies. Rather, our implementations mimic how the treatment of Richards' equation in SNOWPACK and Crocus could be translated into other models.

### C1 Model 4

The implementation of this model is based on the SNOWPACK implementation (Wever et al., 2014) and its publicly available source code (master branch, f89d8c17 commit; checked on 16/12/2024). As explained in the main part of the manuscript, we follow here a sequential scheme, where the process of heat conduction, surface energy fluxes, and phase changes are first solved (and liquid percolation treated in a second step). This yields a LWC field, potentially with wet zones that flows and dry zones where the WRC is undefined. This issue is circumvented in two steps: (i) a small amount of liquid water is introduced, to avoid fully-dry media and (ii) the residual point of the WRC is lowered-down such that the LWC lies on a defined part of the WRC. Concretely, to solve Richards' equation, at the beginning of each adaptive timestep:

1. We compute the residual LWC of the snow material $\theta_r = \max(0, \min(\theta_w - \theta_{acc}/10, \max(0, \min(0.02, \max(\theta_r^{\text{prev}}, 0.75\theta_w))))$, with $\theta_w$ the LWC at the start of the timestep, $\theta_r^{\text{prev}}$ the residual LWC used in the previous timestep, and $\theta_{acc}$ the error aimed for on $\theta_w$ during the non-linear iteration. This step defines a residual LWC that is positive, equals the default value of 0.02 when
<figure><img src="image_39.png" /></figure>

the snow LWC is above 0.027, and tends to 75% of the LWC when it is below 0.027. This corresponds to the lines 395 to 397 in the vanGenuchten.cc source file of SNOWPACK.

2. We compute the matric potential of "dry" layers. For this, we compute the matric potential of layers if they were just above their residual points. We then compute the minimum (maximum in absolute value) of these matric potentials, and define it as the matric potential that "dry" layers should have. This corresponds to the lines 967 to 969 in the ReSolver1d.cc source file of SNOWPACK.

3. We compute the minimum required LWCs that need to be present in the snow layers to reach the "dry" matric potential. This corresponds to line 976 to 999 in the ReSolver1d.cc source file of SNOWPACK.

4. If needed, we increase the LWC of layers that fall below their minimum required LWCs (computed in step 3 above), and recompute the residual LWC as in step 1. In order to ensure mass conservation, this is done by taking mass from the ice. This corresponds to the lines 1017 to 1024 in the ReSolver1d.cc source file of SNOWPACK.

After this step, we obtain consistent fields of LWC and WRC for the whole snowpack. Richards' equation (liquid water conservation under matric flow) is solved using the matric potential as the primary potential. For consistency with the other toy-models, we use a Newton method for the non-linear iterations. As with the other models, the stopping criterion is based on the error when trying to close the PDE, with an additional criterion on mass conservation (as using the matric potential as the primary variable affects mass conservation Celia et al., 1990).

## C2 Model 5

Model 5 is based on the Crocus implementation of Richards' equation (D'Amboise et al., 2017) and its source code (damboise_dev branch, 67bda59d commit; checked on 16/12/2024). It follows the overall same strategy as model 4. The main difference relates to the definition of the residual LWC. Concretely, at the beginning of each (adaptive) timestep:

1. We compute the residual points of the snow material as $\theta_r = \min(0.02, 0.75\theta_w)$, with $\theta_w$ the LWC at the start of the timestep. This results in a residual LWC that equals the standard value of 0.02 in well-wet material ($\theta_w > 0.027$), and that can drop to zero in the case of fully-dry material. This corresponds to the lines 3417 in the snowcro.F90 source file of Crocus.

2. We compute the matric potential of "dry" layers. This corresponds to the minimum matric potential of the layers when they are at the so-called "pre-wetting" level ($\theta_w = 10^{-5}$; D'Amboise et al., 2017). This corresponds to the lines 3408 to 3438 in the snowcro.F90 source file of Crocus.

3. If the LWC of a cell is below the "pre-wetting" level, the LWC is increased so that the matric potential reach the "dry" matric potential. This is done without ensuring mass conservation.

4. The residual points are re-computed as in step 1.

After this step, liquid water matric flow is solved using the same strategy as in model 4.
<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

*Code and data availability.* The implementations of the five toy models have been published as Fourteau (2025), available at
https://doi.org/10.5281/zenodo.14753491.

*Author contributions.* The research was designed by KF, JB, and MD. The mathematical formulation was derived by KF, JB, CC, and MF. The code was developed by KF with help from CC. The manuscript was written by KF with the help of all co-authors.

*Competing interests.* The authors declare having no competing interests.

*Acknowledgements.* We acknowledge the SNOWPACK and Crocus developers, notably for making their source code easily available. We thank Laurent Oxarango, Nander Wever and Michael Lombardo for the fruitful discussions on liquid water transport in snow. Kévin Fourteau's current position and Julien Brondex's past position were funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (IVORI; grant no. 949516). This work was supported by the S-NOW project funded by the Institut des Mathématiques pour la Planète Terre.
<figure><img src="image_42.png" /></figure>

References

Amaziane, B., El Ossmani, M., and Jurak, M.: Numerical simulation of gas migration through engineered and geological barriers for a deep repository for radioactive waste, Computing and visualization in science, 15, 3–20, https://doi.org/10.1007/s00791-013-0196-1, 2012.

655 Baggi, S. and Schweizer, J.: Characteristics of wet-snow avalanche activity: 20 years of observations from a high alpine valley (Dischma, Switzerland), Natural Hazards, 50, 97–108, https://doi.org/10.1007/s11069-008-9322-7, 2009.

Barnhart, T. B., Tague, C. L., and Molotch, N. P.: The Counteracting Effects of Snowmelt Rate and Timing on Runoff, Water Resources Research, 56, e2019WR026634, https://doi.org/10.1029/2019WR026634, 2020.

660 Bartelt, P. and Lehning, M.: A physical SNOWPACK model for the Swiss avalanche warning: Part I: numerical model, Cold Regions Sci. Tech., 35, 123–145, https://doi.org/10.1016/S0165-232X(02)00074-5, 2002.

Bassetto, S., Cancès, C., Enchéry, G., and Tran, Q. H.: Robust Newton Solver Based on Variable Switch for a Finite Volume Discretization of Richards Equation, in: Finite Volumes for Complex Applications IX - Methods, Theoretical Aspects, Examples, edited by Klöfkorn, R., Keilegavlen, E., Radu, F. A., and Fuhrmann, J., vol. 253, pp. 385–393, Springer International Publishing, https://doi.org/10.1007/978-3-030-43651-3_35, 2020.

665 Bassetto, S., Cancès, C., Enchéry, G., and Tran, Q.-H.: Upstream mobility finite volumes for the Richards equation in heterogenous domains, ESAIM: Mathematical Modelling and Numerical Analysis, 55, 2101–2139, 2021.

Bassetto, S., Cancès, C., Enchéry, G., and Tran, Q.-H.: On several numerical strategies to solve Richards’ equation in heterogeneous media with finite volumes, Computational Geosci., 26, 1297–1322, https://doi.org/10.1007/s10596-022-10150-w, 2022.

670
Beaude, L., Brenner, K., Lopez, S., Masson, R., and Smai, F.: Non-isothermal compositional liquid gas Darcy flow: formulation, soil-atmosphere boundary condition and application to high-energy geothermal simulations, Computational Geosciences, 23, 443–470, https://doi.org/10.1007/s10596-018-9794-9, 2019.

Bourgeat, A., Jurak, M., and Smaï, F.: Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository, Computational Geosciences, 13, 29–42, https://doi.org/10.1007/s10596-008-9102-1, 2009.

675
Boutin, C., Auriault, J.-L., and Geindreau, C.: Homogenization of coupled phenomena in heterogenous media, vol. 149, John Wiley & Sons, 2010.

Bouvet, L., Calonne, N., Flin, F., and Geindreau, C.: Multiscale modeling of heat and mass transfer in dry snow: influence of the condensation coefficient and comparison with experiments, The Cryosphere, 18, 4285–4313, https://doi.org/10.5194/tc-18-4285-2024, 2024.

680
Brenner, K. and Cancès, C.: Improving Newton’s Method Performance by Parametrization: The Case of the Richards Equation, SIAM Journal on Numerical Analysis, 55, 1760–1785, https://doi.org/10.1137/16M1083414, 2017.

Brondex, J., Fourteau, K., Dumont, M., Hagenmuller, P., Calonne, N., Tuzet, F., and Löwe, H.: A finite-element framework to explore the numerical solution of the coupled problem of heat conduction, water vapor diffusion, and settlement in dry snow (IvoriFEM v0.1.0), Geoscientific Model Development, 16, 7075–7106, https://doi.org/10.5194/gmd-16-7075-2023, 2023.

685
Brooks, R. and Corey, A. T.: Hydraulic properties of porous media, Hydrology Paper No 3, http://hdl.handle.net/10217/61288, last accessed: 22/10/2024, 1964.

Butcher, J.: Numerical methods for ordinary differential equations, "John Wiley & Sons, Ltd, The Atrium, Southern Gate, Chichester, West Sussex PO19 8SQ, England", 2008.
<figure><img src="image_43.png" /></figure>

Calonne, N., Flin, F., Morin, S., Lesaffre, B., du Roscoat, S. R., and Geindreau, C.: Numerical and experimental investigations of the effective thermal conductivity of snow, Geophysical Research Letters, 38, https://doi.org/10.1029/2011GL049234, 2011.

690 Calonne, N., Geindreau, C., Flin, F., Morin, S., Lesaffre, B., Rolland du Roscoat, S., and Charrier, P.: 3-D image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy, The Cryosphere, 6, 939–951, https://doi.org/10.5194/tc-6-939-2012, 2012.

Calonne, N., Geindreau, C., and Flin, F.: Macroscopic Modeling for Heat and Water Vapor Transfer in Dry Snow by Homogenization, The Journal of Physical Chemistry B, 118, 13393–13403, https://doi.org/10.1021/jp5052535, 2014.

695 Celia, M. A., Bouloutas, E. T., and Zarba, R. L.: A general mass-conservative numerical solution for the unsaturated flow equation, Water Resources Research, 26, 1483–1496, https://doi.org/10.1029/WR026i007p01483, 1990.

D'Amboise, C. J. L., Müller, K., Oxarango, L., Morin, S., and Schuler, T. V.: Implementation of a physically based water percolation routine in the Crocus/SURFEX (V7.3) snowpack model, Geoscientific Model Development, 10, 3547–3566, https://doi.org/10.5194/gmd-10-3547-2017, 2017.

700 Diersch, H.-J. and Perrochet, P.: On the primary variable switching technique for simulating unsaturated-saturated flows, Advances in Water Resources, 23, 271–301, https://doi.org/10.1016/S0309-1708(98)00057-8, 1999.

Eckert, N., Corona, C., Giacona, F., Gaume, J., Mayer, S., van Herwijnen, A., Hagenmuller, P., and Stoffel, M.: Climate change impacts on snow avalanche activity and related risks, Nature Reviews Earth & Environment, pp. 1–21, https://doi.org/10.1038/s43017-024-00540-2, 2024.

705 Farthing, M. W. and Ogden, F. L.: Numerical Solution of Richards’ Equation: A Review of Advances and Challenges, Soil Science Society of America Journal, 81, 1257–1269, https://doi.org/10.2136/sssaj2017.02.0058, 2017.

Fazio, R.: Stiffness in numerical initial-value problems: A and L-stability of numerical methods, Int. J. Math. Educ. Sci. Tech., 32, 752–760, https://doi.org/10.1080/002073901753124619, 2001.

710 Fierz, C., Armstrong, R. L., Durand, Y., Etchevers, P., Greene, E., McClung, D. M., Nishimura, K., Satyawali, P. K., and Sokratov, S. A.: The International Classification for Seasonal Snow on the Ground, UNESCO-IHP, Paris, last accessed: 20/10/2021, 2009.

Forsyth, P., Wu, Y., and Pruess, K.: Robust numerical methods for saturated-unsaturated flow with dry initial conditions in heterogeneous media, Advances in Water Resources, 18, 25–38, https://doi.org/10.1016/0309-1708(95)00020-J, 1995.

Fourteau, K.: Supplementary Material to "Numerical strategies for representing Richards' equation and its couplings in snowpack models", https://doi.org/10.5281/zenodo.14753491, [code, data set], 2025.

715 Fourteau, K., Brondex, J., Brun, F., and Dumont, M.: A novel numerical implementation for the surface energy budget of melting snowpacks and glaciers, Geoscientific Model Development, 17, 1903–1929, https://doi.org/10.5194/gmd-17-1903-2024, 2024.

Hansen, B. B., Aanes, R., and Sæther, B.-E.: Feeding-crater selection by high-arctic reindeer facing ice-blocked pastures, Canadian Journal of Zoology, 88, 170–177, https://doi.org/10.1139/Z09-130, 2010.

720 Hock, R., Rasul, G., Adler, C., Cáceres, B., Gruber, S., Hirabayashi, Y., Jackson, M., Kääb, A., Kang, S., Kutuzov, S., Milner, A., Molau, U., Morin, S., Orlove, B., and Steltzer, H.: High Mountain Areas, in: IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, edited by Po'rtner, H.-O., Roberts, D., Masson-Delmotte, V., Zhai, P., Tignor, M., Poloczanska, E., Mintenbeck, K., Alegri'a, A., Nicolai, M., Okem, A., Petzold, J., Rama, B., and Weyer, N., pp. 131–202, Cambridge University Press, Cambridge, UK and New York, NY, USA, https://doi.org/10.1017/9781009157964.004, 2019.
<figure><img src="image_44.png" /></figure>

Jafari, M., Gouttevin, I., Couttet, M., Wever, N., Michel, A., Sharma, V., Rossmann, L., Maass, N., Nicolaus, M., and Lehning, M.: The Impact of Diffusive Water Vapor Transport on Snow Profiles in Deep and Shallow Snow Covers and on Sea Ice, Frontiers in Earth Science, 8, https://doi.org/10.3389/feart.2020.00249, 2020.

725

Jordan, R. E.: A one-dimensional temperature model for a snow cover: Technical documentation for SNTHERM. 89, http://hdl.handle.net/11681/11677, last accessed: 01/09/2023, 1991.

Kadioglu, S. Y., Nourgaliev, R. R., and Mousseau, V. A.: A Comparative Study of the Harmonic and Arithmetic Averaging of Diffusion Coefficients for Non-linear Heat Conduction Problems, Tech. rep., Idaho National Laboratory, Idaho Falls, Idaho 83415, https://doi.org/10.2172/928087, 2008.

730

Keyes, D. E., McInnes, L. C., Woodward, C., Gropp, W., Myra, E., Pernice, M., Bell, J., Brown, J., Clo, A., Connors, J., Constantinescu, E., Estep, D., Evans, K., Farhat, C., Hakim, A., Hammond, G., Hansen, G., Hill, J., Isaac, T., Jiao, X., Jordan, K., Kaushik, D., Kaxiras, E., Koniges, A., Lee, K., Lott, A., Lu, Q., Magerlein, J., Maxwell, R., McCourt, M., Mehl, M., Pawlowski, R., Randles, A. P., Reynolds, D., Rivière, B., Rüde, U., Scheibe, T., Shadid, J., Sheehan, B., Shephard, M., Siegel, A., Smith, B., Tang, X., Wilson, C., and Wohlmuth, B.: Multiphysics simulations: Challenges and opportunities, The International Journal of High Performance Computing Applications, 27, 4–83, https://doi.org/10.1177/1094342012468181, 2013.

735

Krabbenhøft, K.: An alternative to primary variable switching in saturated-unsaturated flow computations, Advances in Water Resources, 30, 483–492, https://doi.org/https://doi.org/10.1016/j.advwatres.2006.04.009, 2007.

Lenhard, R. J., Parker, J. C., and Mishra, S.: On the Correspondence between Brooks-Corey and van Genuchten Models, Journal of Irrigation and Drainage Engineering, 115, 744–751, https://doi.org/10.1061/(ASCE)0733-9437(1989)115:4(744), 1989.

740

Lewandowska, J., Szymkiewicz, A., Burzyński, K., and Vauclin, M.: Modeling of unsaturated water flow in double-porosity soils by the homogenization approach, Advances in Water Resources, 27, 283–296, https://doi.org/https://doi.org/10.1016/j.advwatres.2003.12.004, 2004.

Libois, Q., Picard, G., France, J. L., Arnaud, L., Dumont, M., Carmagnola, C. M., and King, M. D.: Influence of grain shape on light penetration in snow, The Cryosphere, 7, 1803–1818, https://doi.org/10.5194/tc-7-1803-2013, 2013.

745

Moure, A., Jones, N., Pawlak, J., Meyer, C., and Fu, X.: A Thermodynamic Nonequilibrium Model for Preferential Infiltration and Refreezing of Melt in Snow, Water Resources Research, 59, e2022WR034035, https://doi.org/https://doi.org/10.1029/2022WR034035, 2023.

Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resources Research, 12, 513–522, https://doi.org/https://doi.org/10.1029/WR012i003p00513, 1976.

750

Picard, G. and Libois, Q.: Simulation of snow albedo and solar irradiance profile with the two-stream radiative transfer in snow (TARTES) v2.0 model, EGUSphere, 2024, 1–42, https://doi.org/10.5194/egusphere-2024-1176, 2024.

Quéno, L., Fierz, C., van Herwijnen, A., Longridge, D., and Wever, N.: Deep ice layer formation in an alpine snowpack: monitoring and modeling, The Cryosphere, 14, 3449–3464, https://doi.org/10.5194/tc-14-3449-2020, 2020.

Richards, L. A.: Capillary conduction of liquids through porous mediums, Physics, 1, 318–333, https://doi.org/10.1063/1.1745010, 1931.

755

Sadegh Zadeh, K.: A mass-conservative switching algorithm for modeling fluid flow in variably saturated porous media, Journal of Computational Physics, 230, 664–679, https://doi.org/https://doi.org/10.1016/j.jcp.2010.10.011, 2011.

Sauter, T., Arndt, A., and Schneider, C.: COSIPY v1.3 – an open-source coupled snowpack and ice surface energy and mass balance model, Geosci. Model Dev., 13, 5645–5662, https://doi.org/10.5194/gmd-13-5645-2020, 2020.

Schneebeli, M.: Development and stability of preferential flow paths in a layered snowpack, IAHS Publications-Series of Proceedings and Reports-Intern Assoc Hydrological Sciences, 228, 89–96, 1995.

760
<figure><img src="image_45.png" /></figure>

Sposito, G.: The statistical mechanical theory of water transport through unsaturated soil: 2. Derivation of the Buckingham-Darcy Flux Law, Water Resources Research, 14, 479–484, https://doi.org/https://doi.org/10.1029/WR014i003p00479, 1978.

Steefel, C. I. and McQuarrie, K.: Approaches to modeling of reactive transport in porous media, Reviews in mineralogy, 34, 83–130, https://doi.org/https://doi.org/10.1515/9781501509797-005, 1996.

Torquato, S.: Random Heterogeneous Materials, Springer Science+Business Media New York, 2002.

Tubini, N., Gruber, S., and Rigon, R.: A method for solving heat transfer with phase change in ice or soil that allows for large time steps while guaranteeing energy conservation, The Cryosphere, 15, 2541–2568, https://doi.org/10.5194/tc-15-2541-2021, 2021.

Tyler, N. J. C.: Climate, snow, ice, crashes, and declines in populations of reindeer and caribou (Rangifer tarandus L.), Ecological Monographs, 80, 197–219, https://doi.org/https://doi.org/10.1890/09-1070.1, 2010.

van Dalum, C. T., van de Berg, W. J., Libois, Q., Picard, G., and van den Broeke, M. R.: A module to convert spectral to narrowband snow albedo for use in climate models: SNOWBAL v1.2, Geoscientific Model Development, 12, 5157–5175, https://doi.org/10.5194/gmd-12-5157-2019, 2019.

van Genuchten, M. T.: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Science Society of America Journal, 44, 892–898, https://doi.org/https://doi.org/10.2136/sssaj1980.03615995004400050002x, 1980.

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geoscientific Model Development, 5, 773–791, https://doi.org/10.5194/gmd-5-773-2012, 2012.

Vogel, T., Gerke, H., Zhang, R., and Van Genuchten, M.: Modeling flow and transport in a two-dimensional dual-permeability system with spatially variable hydraulic properties, Journal of Hydrology, 238, 78–89, https://doi.org/https://doi.org/10.1016/S0022-1694(00)00327-9, 2000.

Wever, N., Fierz, C., Mitterer, C., Hirashima, H., and Lehning, M.: Solving Richards Equation for snow improves snowpack meltwater runoff estimations in detailed multi-layer snowpack model, The Cryosphere, 8, 257–274, https://doi.org/10.5194/tc-8-257-2014, 2014.

Wever, N., Schmid, L., Heilig, A., Eisen, O., Fierz, C., and Lehning, M.: Verification of the multi-layer SNOWPACK model with different water transport schemes, The Cryosphere, 9, 2271–2293, https://doi.org/10.5194/tc-9-2271-2015, 2015.

Wever, N., Würzer, S., Fierz, C., and Lehning, M.: Simulating ice layer formation under the presence of preferential flow in layered snow-packs, The Cryosphere, 10, 2731–2744, https://doi.org/10.5194/tc-10-2731-2016, 2016.

Yamaguchi, S., Katsushima, T., Sato, A., and Kumakura, T.: Water retention curve of snow with different grain sizes, Cold Regions Science and Technology, 64, 87–93, https://doi.org/https://doi.org/10.1016/j.coldregions.2010.05.008, 2010.

Yamaguchi, S., Watanabe, K., Katsushima, T., Sato, A., and Kumakura, T.: Dependence of the water retention curve of snow on snow characteristics, Annals of Glaciology, 53, 6–12, https://doi.org/10.3189/2012AoG61A001, 2012.

Yamaguchi, S., Hirashima, H., and Ishii, Y.: Year-to-year changes in preferential flow development in a seasonal snowpack and their dependence on snowpack conditions, Cold Regions Science and Technology, 149, 95–105, https://doi.org/https://doi.org/10.1016/j.coldregions.2018.02.009, 2018.