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)

Julien Brondex, Kévin Fourteau, Marie Dumont, Pascal Hagenmuller, Neige

Calonne, François Tuzet, Henning Löwe

► To cite this version:

Julien Brondex, Kévin Fourteau, Marie Dumont, Pascal Hagenmuller, Neige Calonne, et al.. 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, 2023, 16 (23), pp.7075-7106. 10.5194/gmd-16-7075-2023. hal-04335220

HAL Id: hal-04335220

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

Submitted on 14 Dec 2023

**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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# 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)

Julien Brondex<fnref n="1" />, Kévin Fourteau<fnref n="1" />, Marie Dumont<fnref n="1" />, Pascal Hagenmuller<fnref n="1" />, Neige Calonne<fnref n="1" />, François Tuzet<fnref n="1" />, and Henning Löwe<fnref n="2" />

**Correspondence:** Julien Brondex (julien.brondex@univ-grenoble-alpes.fr)

Received: 17 May 2023 – Discussion started: 8 June 2023

Revised: 17 October 2023 – Accepted: 17 October 2023 – Published: 6 December 2023

**Abstract.** The poor treatment (or complete omission) of water vapor transport has been identified as a major limitation suffered by currently available snowpack models. As vapor and heat fluxes are closely intertwined, their mathematical representation amounts to a system of nonlinear and tightly coupled partial differential equations that are particularly challenging to solve numerically. The choice of the numerical scheme and the representation of couplings between processes are crucial to ensure an accurate and robust solution that guarantees mass and energy conservation while also allowing time steps in the order of 15 min. To explore the numerical treatments fulfilling these requirements, we have developed a highly modular finite-element program. The code is written in Python. Every step of the numerical formulation and solution is coded internally, except for the inversion of the linearized system of equations. We illustrate the capabilities of our approach to tackle the coupled problem of heat conduction, vapor diffusion, and settlement within a dry snowpack by running our model on several test cases proposed in recently published literature. We underline specific improvements regarding energy and mass conservation as well as time step requirements. In particular, we show that a fully coupled and fully implicit time-stepping approach enables accurate and stable solutions with little restriction on the time step.

## 1 Introduction

Over the last few decades, snow models of various complexity have been developed for a myriad of applications, including avalanche forecasting (Morin et al., 2020), water resources management (Magnusson et al., 2015), glacier mass balance assessment (e.g., van Pelt et al., 2012; Sauter et al., 2020), or projections of future climate evolution (Krinner et al., 2018). They range from single-layer snow schemes to detailed snowpack models that provide an explicit description of the vertical distribution of physical properties, such as the Crocus (Brun et al., 1989, 1992) and SNOWPACK (Bartelt and Lehning, 2002) models. However, even the most detailed snow models suffer from major weaknesses (Menard et al., 2021). Their inability to reproduce inverted density gradients as observed in Arctic snowpacks (Domine et al., 2016; Barrere et al., 2017), where strong temperature gradients induce a significant water vapor flux that redistributes the ice mass from the basal to upper layers, has been pinpointed as one of those (Domine et al., 2019). More generally, vapor transport is involved in many processes, such as the redistribution of water vapor isotopes in polar firn (Touzeau et al., 2018) or snow metamorphism (e.g., Sturm and Benson, 1997), with implications for snowpack stability (e.g., Pfeffer and Mrugala, 2002). To address the aforementioned limitations, efforts have been made in two directions: (i) the implementation of ad hoc water vapor transfer modules in existing models (e.g., Touzeau et al., 2018; Jafari et al., 2020) and (ii) the development of stand-alone modules.

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

<footnote n="2">²WSL Institute for Snow and Avalanche Research SLF, Flüelastr. 11, 7260 Davos, Switzerland</footnote>
els to explore the numerical treatment of the coupled heat and water vapor transfer problem (e.g., Simson et al., 2021; Schürholt et al., 2022b).

Heat conduction and water vapor diffusion are two-way coupled: (1) temperature gradients within the snowpack drive water vapor fluxes and (2) vapor fluxes transport latent heat that redistributes energy within the snowpack and feeds back on temperature gradients whenever sublimation/deposition occurs (Yosida et al., 1955; Sturm and Johnson, 1992; Albert and McGilvary, 1992; Pinzer et al., 2012). Furthermore, the microstructure evolves as water vapor deposits on or sublimates from the ice phase. This affects the effective thermal conductivity and effective vapor diffusion coefficient, which in turns feeds back on the energy and vapor fluxes (Yosida et al., 1955; Jaafar and Picot, 1970; Sturm and Johnson, 1992; Calonne et al., 2011; Riche and Schneebeli, 2013). Two mathematical descriptions of the macroscopic heat and water vapor budget in dry snow have been proposed. The first one was derived by Calonne et al. (2014) using the two-scale expansion method. The second one was introduced by Hansen and Foslien (2015) using mixture theory. Both models account for the interactions between energy and vapor fluxes through phase change. However, both models were derived on an invariant microstructure, thus neglecting the feedback of phase change on the distribution of the ice volume fraction, i.e., the fraction of the volume occupied by ice for a given volume of snow.

While the heat equation is at the core of any snowpack model and was therefore implemented at the very early stage of the decades-long development of these models (e.g., Brun et al., 1989; Bartelt and Lehning, 2002), efforts to include water vapor transfer are recent. Because it avoids in-depth modification of the code, which is cumbersome and prone to bug dissemination (Menard et al., 2021), first-order operator splitting is the most natural way to couple newly developed subprocesses to processes already incorporated in previous versions of a model. Basically, this method consists of solving all subprocesses of a coupled process sequentially. Variants of this method have been used by Touzeau et al. (2018) and Jafari et al. (2020) to implement water vapor diffusion in Crocus and SNOWPACK, respectively. However, this approach may require a shorter time step than the main time step of the model. In the present case, the deposition/sublimation rate is controlled by the magnitude of the oversaturation/undersaturation of vapor, which is highly temperature sensitive, and is modulated by the kinetics of the absorption/desorption of water molecules on the ice surfaces (Fourteau et al., 2021). It follows that solving the vapor mass balance and heat equations sequentially is prone to instabilities whenever the time step is too large relative to the considered kinetics. Thus, both Touzeau et al. (2018) and Jafari et al. (2020) have reported time step constraints of 1 s and 1 min, respectively. This is considerably shorter than the time step of 15 min used in Crocus and SNOWPACK for operational avalanche forecasting, for example. Further-

more, some important feedbacks were not accounted for: (i) Touzeau et al. (2018) did not include any phase-change-related latent heat effects in the heat equation and (ii) neither of the two studies considered the evolution of effective parameters (including viscosity) due to deposition/sublimation.

An alternative approach is to treat the heat conduction and water vapor diffusion as a single monolithic process. In two companion papers, Schürholt et al. (2022b) and Simson et al. (2021) developed their own stand-alone models to explore this path. More specifically, Schürholt et al. (2022b) used the Python-based FEniCS finite-element computing platform in order to solve both the models of Hansen and Foslien (2015) and of Calonne et al. (2014) while also accounting for deposition/sublimation feedback on the ice volume fraction field. However, because re-meshing was not supported in FEniCS at the time of the study, they did not include snow settlement. In contrast, Simson et al. (2021) proposed a numerical approach that combines a deforming mesh procedure based on a Lagrangian method to solve for settlement with a classical finite difference method applied on the deformed mesh to solve the model of Hansen and Foslien (2015) only. Both of these works constitute extremely valuable contributions to the understanding of nonlinear feedback. However, a number of limitations are left to be addressed: (1) requirements regarding the time step are still not clearly established and (2) a careful assessment of mass and energy conservation is also lacking.

Snowpack models are always said to be based on mass and energy conservation (e.g., Jordan, 1991; Bader and Weilennann, 1992; Brun et al., 1989, 1992; Bartelt and Lehning, 2002; Sauter et al., 2020). While this is usually true for their mathematical formulation, any numerical implementation might cause violations of mass and energy conservation. Commonly, these numerical details do not receive sufficient attention. This is reflected, for example, by the fact that none of the previous snow intercomparison projects (Slater et al., 2001; Etchevers et al., 2004; Essery et al., 2009; Krinner et al., 2018) contain an assessment of the accuracy of mass and energy conservation in the numerical implementation. This contrasts with dedicated numerical experiments conducted in other intercomparison projects – for example, in the climate modeling community (Irving et al., 2021). Numerical errors are frequently argued to be small. However, as snowpack model runs must cover seasons or centuries, even small numerical inconsistencies can cause drifts or lead to significant problems on these timescales. As we will show in this paper, achieving strict mass and energy conservation is particularly subtle in the highly coupled, nonlinear situations outlined above.

In this study, we aim to pinpoint the most suitable numerical treatments based on time step requirement, conservation of energy, and conservation of mass criteria. To this end, we unify previous model developments within a comprehensive, stand-alone, finite-element core written in Python (Brondex et al., 2023). Contrary to Schürholt et al. (2022b), each step
of the numerical formulation and solution is coded internally, except for the inversion of the linearized system of equations which relies on the standard NumPy linear-algebra library. In this way, we have complete control of every detail of the numerical recipe and can thus explore various solving strategies. We demonstrate the improvements within established benchmark scenarios and carve out the origin of errors in the numerical solution of the conservation laws. We also discuss the treatment of boundary conditions (BCs) on vapor.

The paper is organized as follows: in Sect. 2, we introduce the mathematical models derived by Calonne et al. (2014) and Hansen and Foslien (2015) and underline specific issues; in Sect. 3, we go through the details of our numerical implementation, specifying main differences compared to previous work; in Sect. 4, our model is tested on numerical benchmarks and appropriate numerical approaches are highlighted; and Sect. 5 summarizes our work and is an opening on the implications of our findings for future work.

## 2 Mathematical models

In this section, we present the mathematical models derived by Calonne et al. (2014) and Hansen and Foslien (2015) and point out relevant issues.

### 2.1 The Calonne model

Starting from the physical phenomena occurring at the pore scale – specifically (i) the heat conduction through air and ice, (ii) the water vapor diffusion in the pore space, and (iii) the sublimation and deposition of vapor at the ice-pore interface – Calonne et al. (2014) used the two-scale expansion method in order to derive a closed system of equations governing the heat and water vapor budget at the macroscopic scale. The main advantage of this approach is that the exact expression of the effective properties (such as the snow thermal conductivity) and of the source terms naturally arises as the macroscopic system of equations is derived. However, it must be stressed that Calonne et al. (2014) did not include any equation governing the evolution of the pore space at the microscale related to water molecules depositing on or sublimating from the ice-pore interface. As a consequence, this macroscopic model implicitly assumes an invariant microstructure, i.e., the ice volume fraction does not evolve over time. In addition, the two-scale expansion method is unsuited for high reaction rates (Bourbatache et al., 2021), i.e., when vapor deposits easily on the ice. Thus, the domain of applicability of the macroscopic model proposed by Calonne et al. (2014) is bounded to cases in which the crystal growth velocity due to the deposition (or sublimation) of water vapor molecules at the ice-pore interface is limited by the characteristic time of reaction rather than by the diffusion one (i.e., a low Damköhler number; Bourbatache et al., 2021). Under these two assumptions, the macroscopic system of

equations is as follows:

$$ \left\{ \begin{array}{l} (\rho C_p)^{\text{eff}} \partial_t T - \nabla \cdot (k^{\text{eff}} \nabla T) = L_m c \\ (1 - \Phi_i) \partial_t \rho_v - \nabla \cdot (D^{\text{eff}} \nabla \rho_v) = -c \end{array} \right. , \qquad (1) $$

where T is the temperature, $\rho_v$ is the water vapor density, $\Phi_i$ is the ice volume fraction, $L_m$ is the specific latent heat of sublimation, $(\rho C_p)^{\text{eff}}$ is the effective heat capacity, $k^{\text{eff}}$ is the effective heat conductivity, and $D^{\text{eff}}$ is the effective vapor diffusion coefficient. The deposition rate c is given by the following expression:

$$ c = s \alpha v_{\text{kin}} (\rho_v - \rho_v^{\text{eq}}), \qquad (2) $$

where s is the surface area density per unit volume, which is assumed to be constant, consistent with the implicit invariant-microstructure assumption; $v_{\text{kin}} = \sqrt{(k_B T)/(2\pi m_{H_2O})}$ is the kinetic velocity related to the velocity of water molecules in the pore space ($k_B$ is the Boltzmann constant and $m_{H_2O}$ is the mass of a water molecule); $\rho_v^{\text{eq}}$ is the saturation vapor density; and $\alpha$ is the sticking coefficient of water molecules on the ice surface. Referring back to the considerations raised above regarding the domain of applicability of this macroscopic model, it appears that the latter is valid provided that $\alpha \approx 10^{-3}$ or less (Fourteau et al., 2021).

As the saturation water vapor density $\rho_v^{\text{eq}}$ is a nonlinear function of T, the macroscopic model proposed by Calonne et al. (2014) amounts to a system of two, two-way-coupled nonlinear diffusion-reaction equations that must be solved for T and $\rho_v$. Because we use parameterizations that neglect the dependency of the effective parameters on T (Appendix A), the only nonlinearity of the system arises from the source terms through $\rho_v^{\text{eq}}$ and $v_{\text{kin}}$. In what follows, we refer to this model as the Calonne model.

### 2.2 The Hansen model

Using mixture theory, Hansen and Foslien (2015) derived another macroscopic heat and water vapor conservation model. Contrary to Calonne et al. (2014), they assumed the deposition/sublimation of water vapor to be fast enough so that small oversaturation/undersaturation events in the pore space are corrected almost instantaneously by the deposition/sublimation of water molecules. Mathematically, such a situation arises when the product $\alpha v_{\text{kin}}$ becomes very large. In this case, the deposition rate corresponds to the rate at which this deposition/sublimation of molecules must occur so that the water vapor density is permanently and instantaneously restored to its saturated value, i.e., so that $\rho_v(T) = \rho_v^{\text{eq}}(T)$ at any time. As for the Calonne model, the equations of Hansen and Foslien (2015) are derived at constant microstructure. Based on these two assumptions, the respective macroscopic energy conservation and water vapor mass balance can be
written as follows:

$$ \left\{ \begin{array}{l} (\rho C_p)^{\text{eff}} \partial_t T - \nabla \cdot (k^{\text{eff}} \nabla T) = L_m c \\ (1 - \Phi_i) \frac{d\rho_v^{\text{eq}}}{dT} \partial_t T - \nabla \cdot \left( D^{\text{eff}} \frac{d\rho_v^{\text{eq}}}{dT} \nabla T \right) = -c \end{array} \right. , \quad (3) $$

with the underlying assumption that

$$ \rho_v = \rho_v^{\text{eq}}. \quad (4) $$

In Hansen and Foslien (2015), the effective parameters $k^{\text{eff}}$ and $D^{\text{eff}}$ are intuited from synthetic microstructures and are not a direct by-product of the mixture theory. The system of Eq. (3) can be cast into a single expression:

$$ \left( (\rho C_p)^{\text{eff}} + (1 - \Phi_i) L_m \frac{d\rho_v^{\text{eq}}}{dT} \right) \partial_t T - \nabla \cdot \left[ \left( D^{\text{eff}} L_m \frac{d\rho_v^{\text{eq}}}{dT} + k^{\text{eff}} \right) \nabla T \right] = 0. \quad (5) $$

In the end, the model of Hansen and Foslien (2015) consists of a single prognostic equation for $T$ that has the form of a nonlinear diffusion equation. The nonlinearity comes from the dependence of the derivative of the saturation water vapor density, which appears in the apparent heat capacity and apparent thermal conductivity, on temperature. The deposition/sublimation rate $c$ can then be diagnosed from one of the expressions in Eq. (3). In what follows, we refer to this model as the Hansen model.

## 2.3 General form of the equations governing the heat and vapor budgets in evolving pore space

Despite being derived under different assumptions, the Calonne and Hansen models comprise similarities. In both of these works, the total energy budget accounts for the contribution of vapor transport and of the latent heat that is released or absorbed whenever water vapor deposits or sublimates, which affects the water vapor mass balance in return. As noted by Schürholt et al. (2022b), if the differences regarding the parameterizations of the effective parameters are put aside, the system of Eq. (3) can be derived from the system of Eq. (1) under the assumption that $\rho_v = \rho_v^{\text{eq}}$. Fundamentally, both systems of equations consist of an equation for the conservation of heat that includes a source term proportional to the deposition rate and an equation for the conservation of vapor that includes the same deposition rate as a sink term. However, the systems of equations differ in how they are closed: the Calonne model is closed by computing the deposition rate $c$ as a first-order reaction, whereas the Hansen model is closed by assuming water vapor saturation.

Both models have been derived assuming an invariant microstructure. This restriction can be lifted by realizing that the first terms of both the heat and vapor conservation equations correspond to the time derivatives of the total heat and

total vapor content, respectively. These terms can then simply be rewritten to include the effect of an evolving ice volume fraction on the temporal evolution of these two quantities. A more subtle issue lies in the fact that both works neglected the settlement of the snowpack and the resulting advection of material quantities. This problem is usually circumvented by solving the settlement and heat/vapor budget problems separately. This corresponds to the Eulerian-Lagrangian framework described by Simson et al. (2021): (i) the settlement of the snowpack is inferred adopting a Lagrangian point of view, (ii) the computed settlement is used to deform the numerical mesh, and (iii) the heat and vapor equations are solved on the obtained mesh using an Eulerian approach. With this procedure, the contribution of advection in the evolution of the temperature, vapor, and deposition rate fields is accounted for during the mesh deformation step and, thus, vanishes from the heat/vapor conservation equations. In view of these considerations, the general form of the equations governing the heat, vapor, and ice budgets in snow in an Eulerian framework can be written as follows:

$$ \left\{ \begin{array}{l} \partial_t h_{\text{ice}} - \nabla \cdot (k^{\text{eff}} \nabla T) = L_m c \\ \partial_t [(1 - \Phi_i) \rho_v] - \nabla \cdot (D^{\text{eff}} \nabla \rho_v) = -c \end{array} \right. , \quad (6) $$

where $h_{\text{ice}}$ is the heat content of snow and $\rho_i$ is the density of ice. We stress that, contrary to Eq. (1) or Eq. (3), the accumulation term of the heat equation is written in terms of heat content and not in terms of temperature. This heat content is related to temperature as follows:

$$ \frac{\partial h_{\text{ice}}}{\partial T} = (\rho C_p)^{\text{eff}}. \quad (7) $$

One might be tempted to express $\partial_t h_{\text{ice}}$ as $(\rho C_p)^{\text{eff}} \partial_t T$ by means of the chain rule. As thoroughly explained in publications such as Tubini et al. (2021), while this is valid in the continuous partial differential equation (PDE), this is not necessarily the case for the discrete domain. Specifically, application of the chain rule if $(\rho C_p)^{\text{eff}}$ is not constant during the discrete solution of the equation will result in the non-conservation of energy (e.g., Celia et al., 1990; Casulli and Zanolli, 2010; Tubini et al., 2021). Therefore, we choose to explicitly keep the heat content $h_{\text{ice}}$ in the heat budget and, thus, to express this equation in its so-called “mixed form” (similarly to Celia et al., 1990).

The system of Eqs. (6) and (7) contains five unknowns ($h_{\text{ice}}$, $T$, $\rho_v$, $\Phi_i$, and $c$) and can be closed either with Eq. (2) from the Calonne model or with Eq. (4) from the Hansen model. We stress that this system of equations alone is not sufficient to compute the evolution of the ice volume fraction. The contribution of mechanical settlement of snow on the latter must also be accounted for (Simson et al., 2021). Although this process could be described in an Eulerian framework through the use of an advection term (i.e., $\partial_t \Phi_i + \nabla \cdot$
($v\Phi_i$) = 0, where $v$ is the settling velocity), as mentioned above, adopting a Lagrangian point of view to describe the deformation of the snowpack provides a more natural framework to compute the inherent ice volume fraction increase. In the absence of phase change, the evolution of $\Phi_i$ is only due to settlement. If, additionally, the dependence of the settling rate on temperature through viscosity is neglected, the equations governing heat and vapor and the equations governing the ice volume fraction become partly independent and can be solved sequentially.

Finally, the compaction of snow leads to a reduction in the pore space. As the air is free to escape during the settlement of the snowpack, this decrease in porosity results in a net loss of dry air and vapor. If the effective heat capacity of snow ($\rho C_p$)⁵⁷ is computed considering the heat carried by dry air, this loss of dry air leads to a net loss of energy. While physically relevant, this loss of energy though air ejection can be overlooked by considering that only the ice phase carries energy in snow (Appendix A3). This assumption is justified by the small volumetric heat capacity of dry air compared with that of ice. In contrast, the net loss of vapor is fully part of the considered problem and should, therefore, be kept in mind when closing the vapor budget from one time step to the next.

## 3 Numerical implementation

In this section, we go through the most important features of our numerical implementation and underline the main differences compared with the published approaches of Simson et al. (2021) and Schürholt et al. (2022b). All variables, parameters, and constants used in our model are summarized in Table 1.

### 3.1 Numerical strategy

The spatial discretization of the PDEs of the model is based on the finite-element method (FEM). For the temporal discretization, we use a standard first-order implicit Euler method. Generalities about the FEM and the temporal discretization scheme are presented in Appendix B. This includes a description of the implementation of BCs as well as the calculation of residuals that is necessary to close the energy budget.

#### 3.1.1 Coupling scheme

One of the main goals of our work is to highlight the importance of a coupled numerical solution of the heat and vapor transport processes. Details of the various approaches that have been implemented and compared to reach this conclusion are given in Sect. 3.2. In contrast, the ice mass conservation is systematically solved in a separated step following the same first-order operator-splitting approach as Simson et al. (2021) and Schürholt et al. (2022b). However, as

demonstrated in Sect. 2, heat and vapor transfer and deposition are two-way coupled: perturbation of the $T$ and/or $\rho_v$ fields affects the deposition rate field, which in turns changes the distribution of $\Phi_i$; this feeds back on the distributions of the $\Phi_i$-dependent parameters ($\rho C_p$)⁵⁷, $k$⁵⁷ and $D$⁵⁷, impacting the fields of $T$ and $\rho_v$. Therefore, one has to be aware that solving these two processes sequentially will necessarily introduce some error. Specifically, as the energy and vapor mass budget are solved assuming a constant microstructure, the consecutive modification of the latter through deposition/sublimation breaks the previously computed energy and vapor mass budgets. This results in nonphysical energy and vapor mass sources/sinks that we referred to as energy/-mass leakage. However, in most of the natural configurations, the variation in the ice mass related to deposition within targeted time steps in the order of 15 min is expected to remain negligible compared with the total ice mass of the simulated domain. It follows that the energy leakage remains limited, and its associated error in the global energy budget is well identified. Note also that, because our settlement scheme is designed to conserve the ice mass perfectly in the absence of phase change (see Sect. 3.3), there is no such problem of energy leakage for the settlement-induced evolution of $\Phi_i$. For all these reasons, we consider that the operator-splitting approach is acceptable with respect to the ice mass conservation equation.

Figure 1a summarizes the general structure of the model. Within one time step, we first solve the heat and vapor transfer process, which can be modeled using either the system of Calonne or that of Hansen (Fig. 1b). Independently of the considered system, this is done with the distribution of $\Phi_i$ from the previous time step. Next, the nodal field of stress is updated from the weight contained in all overlying elements, with the latter being calculated from the distribution of $\Phi_i$ from the previous time step. Then, the settlement solver is executed in order to consistently update the mesh node positions and the element-wise field of $\Phi_i$. Finally, a solver is executed to diagnose the amount of energy contained in the domain from the nodal fields of $T$ and $\rho_v$ (for the system of Calonne and the system of Hansen in its $T$ form, i.e., when using $T$ as the prognostic variable) or from the nodal field of enthalpy (for the system of Hansen in its mixed form). At the very end of the time step, the global ice mass and energy budgets are evaluated to check for conservation. Note that any of the solvers can be activated or deactivated depending on the processes of interest. In particular, the settlement-related and deposition-related evolution of $\Phi_i$ can easily and independently be switched on or off.

#### 3.1.2 Computational domain and notation

As illustrated in Fig. 2, the snowpack is vertically discretized on a one-dimensional finite-element grid. The $z$ axis is oriented upward, with $z = 0$ corresponding to the soil-snow interface. The initial positions of nodes are based on the user-
<table><caption>**Table 1.** List of variables, parameters, and constants used in our model.</caption><thead><tr><th>Symbol</th><th>Name</th><th>Equation/Value</th><th>Unit</th></tr></thead><tbody><tr><td colspan="4"><strong>Nodal variables</strong></td></tr><tr><td>c</td><td>Deposition rate</td><td>Eqs. (2), (3)</td><td>kg m<sup>-3</sup> s<sup>-1</sup></td></tr><tr><td>H</td><td>Enthalpy content</td><td>Eq. (13)</td><td>J m<sup>-3</sup></td></tr><tr><td>T</td><td>Temperature</td><td>Eqs. (1), (5), (13)</td><td>K</td></tr><tr><td>ρ<sub>v</sub></td><td>Vapor density</td><td>Eqs. (1), (4)</td><td>kg m<sup>-3</sup></td></tr><tr><td>σ</td><td>Stress</td><td>Eq. (16)</td><td>Pa</td></tr><tr><td colspan="4"><strong>Element-wise variables</strong></td></tr><tr><td>E</td><td>Energy per element</td><td>Eq. (12)</td><td>J m<sup>-2</sup></td></tr><tr><td>Φ<sub>i</sub></td><td>Ice volume fraction</td><td>Eq. (19)</td><td>–</td></tr><tr><td colspan="4"><strong>Parameters</strong></td></tr><tr><td>D<sup>eff</sup></td><td>Effective vapor diffusion coefficient</td><td>Eqs. (A1), (A2)</td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>k<sup>eff</sup></td><td>Effective heat conductivity</td><td>Eqs. (A3), (A4)</td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>ν<sub>kin</sub></td><td>Kinetic velocity</td><td>Eq. (2)</td><td>m s<sup>-1</sup></td></tr><tr><td>s</td><td>Surface area density per unit volume</td><td>3770</td><td>m<sup>-1</sup></td></tr><tr><td>α</td><td>Sticking coefficient (value by default)</td><td>5 × 10<sup>-3</sup></td><td>–</td></tr><tr><td>η</td><td>Effective viscosity</td><td>Eq. (A7)</td><td>Pa s<sup>-1</sup></td></tr><tr><td>(ρC<sub>p</sub>)<sup>eff</sup></td><td>Effective heat capacity</td><td>Eq. (A5)</td><td>J m<sup>-3</sup> K<sup>-1</sup></td></tr><tr><td>ρ<sub>v</sub></td><td>Saturation water vapor density</td><td>Eq. (A6)</td><td>kg m<sup>-3</sup></td></tr><tr><td colspan="4"><strong>Constants</strong></td></tr><tr><td>C<sub>a</sub></td><td>Air heat capacity</td><td>1005</td><td>J m<sup>-3</sup> K<sup>-1</sup></td></tr><tr><td>C<sub>i</sub></td><td>Ice heat capacity</td><td>2000</td><td>J m<sup>-3</sup> K<sup>-1</sup></td></tr><tr><td>D<sub>0</sub></td><td>Vapor diffusion coefficient in air</td><td>2 × 10<sup>-5</sup></td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>g</td><td>Gravitational constant</td><td>9.81</td><td>m s<sup>-2</sup></td></tr><tr><td>k<sub>a</sub></td><td>Air heat conductivity</td><td>0.024</td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>k<sub>i</sub></td><td>Ice heat conductivity</td><td>2.3</td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>k<sub>B</sub></td><td>Boltzmann constant</td><td>1.38 × 10<sup>-23</sup></td><td>J K<sup>-1</sup></td></tr><tr><td>L<sub>m</sub></td><td>Specific latent heat of sublimation</td><td>2835333</td><td>J kg<sup>-1</sup></td></tr><tr><td>m<sub>H₂O</sub></td><td>Mass of a water molecule</td><td>2.991507 × 10<sup>-26</sup></td><td>kg</td></tr><tr><td>ρ<sub>i</sub></td><td>Ice density</td><td>917</td><td>kg m<sup>-3</sup></td></tr></tbody></table>

prescribed initial snowpack height so that the mesh is initially uniform. In the presence of mechanical settlement, the mesh will deform nonuniformly. Note that the application of the FEM to nonuniform mesh is straightforward. In this work, the problem of re-meshing has not been investigated, but it is discussed in Sect. 5. We distinguish between variables defined at nodes and variables defined in an element-wise manner (Table 1). In what follows, the former are identified with the subscripts *k*, corresponding to the node number, whereas the latter are identified with the subscripts *k + 1/2*, corresponding to the element number. The element *k + 1/2* is located between nodes *k* and *k + 1*. The total number of nodes is denoted as *N*<sub>*z*</sub>.

## 3.2 Numerical solution of heat and water vapor transfer

### 3.2.1 Solution of the system of Calonne

As illustrated in Fig. 1b, there are two possible approaches to solve the system of Eqs. (1) and (2): either the whole system can be solved in a coupled way (green boxes in Fig. 1b) or the three equations can be solved sequentially (red boxes in Fig. 1b). Two strategies are investigated for the coupled approach. The first strategy, denoted as CC_3DOF, consists of solving the whole system (Eqs. 1 and 2) at once for the solution vector **u** = (T, ρ<sub>v</sub>, c). This strategy corresponds to the one adopted by Schürholt et al. (2022b) to treat the Calonne system. The second strategy, denoted as CC_2DOF, consists of injecting Eq. (2) into the right-hand-side terms of Eq. (1) and solving the latter system for the solution vector **u** = (T, ρ<sub>v</sub>). The *c* field is then diagnosed from the obtained
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Panel (a) shows the general structure of the model. Panel (b) illustrates the three approaches implemented to compute the fields of temperature, water vapor, and deposition rate for the coupled solution of the Calonne system, the decoupled solution of the Calonne system, and the solution of the Hansen system. For each of these approaches, two strategies have been considered and are described in the text and Appendix C.</figcaption></figure>

$\rho_v$ field through the water vapor mass balance equation. Two strategies are also considered for the sequential treatment of the system. In the strategy denoted as CD_PC, the heat equation is solved first with a source term that is fixed from the $c$ field computed at the previous time step. Then, the water vapor mass balance is solved under its diffusion-reaction form: the right-hand side of the equation is replaced by its literal expression (Eq. 2), which introduces a reaction term on $\rho_v$, while the previously evaluated $T$ field is used to compute the source term stricto sensu, i.e., $s\alpha v_{\text{kin}}(T)\rho_v^{\text{eq}}(T)$. Finally, the $c$ field is updated as the closure of the water vapor mass balance from the computed $\rho_v$ field. The second strategy, denoted as CD FD, is similar to the CD_PC strategy, except that both the heat and water vapor mass balance equations have their right-hand side fixed from the $c$ field computed at the previous time step. The latter is updated in a third step from Eq. (2) using the obtained $T$ and $\rho_v$ fields. The resulting matrix systems are summarized in Appendix C.

In order to manipulate standard mass matrices in which all nonzero terms correspond to the product of test and shape functions, FEM software packages are prone to dealing with

a particular form of the heat equation in which the thermal diffusivity $k^{\text{eff}}/(\rho C_p)^{\text{eff}}$ is assigned to the flux divergence operator, rather than keeping the heat capacity $(\rho C_p)^{\text{eff}}$ as a factor of the accumulation term. Similarly, it can be tempting to divide the water vapor mass balance equation by the factor $(1 - \Phi_i)$ assigned to the accumulation term in order to deal with the generic form of a diffusion-reaction equation. These operations are performed, for example, in the FEniCS code of Schürholt et al. (2022a). Practically, this corresponds to assigning the inverse of these factors to the stiffness matrix and force vector rather than keeping them in the mass matrix. However, as soon as these factors are nonuniform, which is the case for snow in general as $\Phi_i = f(z)$ and $(\rho C_p)^{\text{eff}} = f(\Phi_i)$, the conservative forms of Eq. (1) are not equivalent to these reformulations. Therefore, in order to preserve the conservative properties (at given microstructure) of the system of Calonne in the discrete domain, we assign the factors $(\rho C_p)^{\text{eff}}$ and $(1 - \Phi_i)$ directly to the mass matrix (Appendix C). In the following, we refer to this as the proper form of the mass matrix.
<figure><img src="image_4.png" /><figcaption>Figure 2. Illustration of the computational domain. The nature of the various model variables is summarized in Table 1.</figcaption></figure>

By nature, the system of Eq. (1) respects the maximum principle in the continuous domain. From a physical perspective, the maximum principle states that, in the absence of phase change, the maximum value of T (or $\rho_v$) is reached at either the initial time or at the boundaries (Protter and Weinberger, 2012). In order to get a uniform convergence of the numerical solution to the exact one and avoid nonphysical extrema in the interior of the domain, it can be shown that a discrete counterpart of this principle, the so-called discrete maximum principle (DMP), must be fulfilled (Ciarlet, 1973). However, in some situations, the FEM is inclined to violate the DMP due to, among others, the treatment of time derivatives (e.g., Celia et al., 1990) and/or of reaction terms (e.g., John and Schmeyer, 2008). A sufficient condition for the DMP to be respected is that the matrix A of the system (Eq. B2) has the following properties: (i) all diagonal terms are positive, (ii) all off-diagonal terms are negative or zero, and (iii) the row sums are positive (John and Schmeyer, 2008). Thus, a commonly used method to enforce the matrix system to satisfy the DMP without adding further constraints on the time step is to lump the mass matrix and/or the reactive part of the stiffness matrix. This operation consists of replacing the diagonal term in each row of the considered matrix with the sum of all terms in the row while all off-diagonal terms are simultaneously forced to zero (e.g., Milly, 1985; Celia et al., 1990; John and Schmeyer, 2008; Thomée, 2015).

In the present work, this method is adopted whenever spurious oscillations show up. As we will show, although the obtained solution fields are very slightly modified, this method enables the removal of spurious oscillations without affecting energy conservation.

As mentioned in Sect. 2.1, the nonlinearity of the system of Calonne is related to the dependence of the source terms, through the parameters $\rho_v^{eq}$ and $v_{kin}$, on temperature (Appendix A4). As a consequence, this nonlinearity vanishes whenever these parameters are computed from a T field obtained in a separated step. This is the case for both strategies based on the decoupled solution for the system. In contrast, as soon as a coupled solution of the system is considered, a linearization procedure is required. For both strategies of the coupled approach, we implement a linearization algorithm which mixes a Picard method for $v_{kin}$ and a Newton method for $\rho_v^{eq}$. Practically, when solving the system (Eq. B2) at iteration $k+1$ for $u^{k+1}$, the value of $v_{kin}^{k+1}$ is fixed using the temperature field obtained at previous iteration, i.e., $v_{kin}^{k+1} = v_{kin}(T^k)$. Within the same iteration step, the value of $\rho_v^{eq,k+1}$ is evaluated as follows:

$$ \rho_v^{eq,k+1} = \left. \frac{d\rho_v^{eq}}{dT} \right|_{T=T^k} (T^{k+1} - T^k) + \rho_v^{eq}(T^k). \quad (8) $$

The first iteration corresponds to the solution vector obtained at the previous time step $T^n$. The convergence criterion of the linearization algorithm reads as follows:

$$ 2 \frac{\|u^{k+1} - u^k\|}{\|u^{k+1} + u^k\|} < \varepsilon, \quad (9) $$

where the tolerance criterion is set to $\varepsilon = 10^{-5}$ by default. This convergence criterion corresponds to the one adopted by default in the open-source Elmer multi-physics finite-element code (https://www.csc.fi/web/elmer, last access: 25 September 2023). Among many other applications, Elmer is one of the most popular codes for numerical simulations in the field of glaciology (Gagliardini et al., 2013). A number of other measures of the change in the solution between two consecutive iterations could be used instead of Eq. (9). We stress that the convergence criterion is used only to stop the nonlinear iterations, it has no impact on the intermediate iterations. In all simulations that are presented in Sect. 4, the systems at stake are not strongly nonlinear and only one to three nonlinear iterations are needed to satisfy the convergence criterion. Most importantly, the convergence is smooth, and the solution is already very close to its final (converged) value after the first nonlinear iteration (i.e., the change in the solution between, e.g., iterations one and two are very small compared with the change between the first guess and iteration one). It follows that the obtained solution fields are unlikely to show significant sensitivity to the formulation of the convergence criterion for the numerical experiments presented in this study. This point has been con-
firmed by dedicated sensitivity tests (not shown in this paper).

Once the linearization algorithm has converged, the residuals of the system are evaluated to assess the sensible heat $\Delta t F_T^{\partial\Gamma}$ and water vapor mass $\Delta t F_P^{\partial\Gamma}$ that have entered or left the domain over the time step. The energy budget is evaluated at the very end of the time step by comparing the evolution of the energy contained in the domain since the previous time step $\Delta E_{\Omega}$ to the aforementioned boundary fluxes. In practice, we compute the following:

$$E_{\text{leak}} = \Delta E_{\Omega} - \Delta t (F_T^{\partial\Gamma} + L_m F_P^{\partial\Gamma}) - \Delta t L_m F_{P,\text{set.}}^{\partial\Gamma}, \quad (10)$$

where $E_{\text{leak}}$ is the energy leakage, which is zero if energy is conserved. As mentioned in Sect. 2, because the contribution of dry air in the effective heat capacity of snow is neglected, the dry air expelled from the snowpack in response to settlement does not affect the energy budget. On the contrary, the latent energy loss due to water vapor leaving the domain as the snowpack settles must be accounted for through the term $\Delta t L_m F_{P,\text{set.}}^{\partial\Gamma}$. The water vapor mass leaving the snowpack over $\Delta t$ is diagnosed in the settlement solver as follows:

$$\Delta t F_{P,\text{set.}}^{\partial\Gamma} = \sum_{k=1}^{N_{\text{t}-1}} \left[ \int_{\zeta=z_k^{n+1}}^{z_{k+1}^{n+1}} \rho_v^{n+1} d\zeta - \int_{\zeta=z_k^{n+1}}^{z_{k+1}^{n+1}} \rho_v^{n+1} d\zeta \right]. \quad (11)$$

The total energy contained in the domain $E_{\Omega}$ is diagnosed as the sum of the energies contained in each element in the last executed solver (Fig. 1a). The energy contained in the element $k+1/2$ is calculated as follows:

$$E_{k+1/2} = \int_{\zeta=z_k^{n+1}}^{z_{k+1}^{n+1}} \left[ (\rho C)_{\text{eff}}^{n+1} (T^{n+1} - T_0) + L_m (1 - \Phi_i^{n+1}) \rho_v^{n+1} \right] d\zeta, \quad (12)$$

where $T_0$ is an arbitrary reference temperature, which we set to $T_0 = 273$ K.

### 3.2.2 Solution of the system of Hansen

For the system of Hansen, two implementation strategies are also investigated (yellow boxes in Fig. 1b). In the first strategy, denoted as H_MF, the prognostic equation that relates to T evolution is treated using its mixed form. Concretely, this consists of solving the following system for the solution vector $u = (H, T)$:

$$\left\{ \begin{array}{l} \partial_t H - \partial_z \left[ \left(D^{\text{eff}} L_m \frac{d\rho_v^{\text{eq}}}{dT} + k^{\text{eff}}\right) \partial_z T \right] = 0 \\ H = (\rho C_p)^{\text{eff}} (T - T_0) + (1 - \Phi_i) L_m \rho_v^{\text{eq}}, \end{array} \right. \quad (13)$$

where H is the enthalpy content defined as a nodal variable. Because H is a $\Phi_i$-dependent prognostic variable, it must

be updated within the settlement solver after $\Phi_i$ has been updated. In order to illustrate the problem of violation of energy conservation when a chain rule is performed in the discrete domain (Sect. 2.2), we consider a second strategy, denoted as H_TF, which consists of solving the prognostic equation of the Hansen system using its usual T form (Eq. 5). The resulting matrix systems are summarized in Appendix C. The H_TF strategy is the one adopted by Simson et al. (2021) and Schürholt et al. (2022b). Independently of the chosen strategy, the field of $\rho_v$ is deduced from the obtained T field assuming that water vapor density is saturated everywhere. Then, the second equation of (Eq. 3) is solved for the field of c.

For the T form of the prognostic equation, we are careful to assign the apparent heat capacity directly to the mass matrix (proper form of the mass matrix). Note that ambiguities do not arise for the mixed form of the equation, as it does not include any multiplying factor in the accumulation term. Furthermore, as for the Calonne system, spurious oscillations related to the violation of the DMP can occur. Again, this difficulty is overcome by lumping the mass matrices of both the prognostic equation that relates to T evolution (or H evolution for the mixed form) and the diagnostic equation of c when necessary. As mentioned in Sect. 2.2, both the mixed form and T form of the prognostic equation are nonlinear due to the dependence of $d\rho_v^{\text{eq}}/dT$ on T (Appendix A4). This nonlinearity is treated via the implementation of a Picard linearization loop. As for the coupled Calonne system, the tolerance criterion is set to $\epsilon = 10^{-5}$ by default. Note that the approach of Simson et al. (2021), which consists of linearizing the equation by fixing the value of $d\rho_v^{\text{eq}}/dT$ from the T field obtained at the previous time step, is equivalent to performing a single Picard iteration. In this case, the obtained solution does not correspond to the implicit solution of the problem and, thus, does not necessarily possess its stability features.

As for the Calonne system, the residuals are evaluated immediately after the convergence of the linearization algorithm. Note that, because the Hansen system contains only one prognostic equation that mixes up the contributions of latent and sensible heat to the energy budget, the obtained residuals correspond to the total enthalpy fluxes at boundaries, without any reference to how these fluxes are split between sensible and latent heat fluxes. Again, these boundary fluxes are used to close the energy budget. The total energy contained in the domain is directly obtained through integration of the nodal enthalpy over the elements for the mixed-form case. For the T-form case, it is evaluated from T based on Eq. (12) using the assumption of saturated water vapor density. As for the Calonne system, the energy budget accounts for energy loss related to the water vapor mass leaving the snowpack when settlement is activated through Eq. (11).

Although the c field is diagnosed a posteriori from the obtained T field based on the second equation of Eq. (3), this solution is also based on the FEM, and a boundary flux term
naturally appears in the force vector through integration by parts of the divergence operator. In the absence of any prescription from the user, the natural BC applies and this term is forced to be zero. However, the existence of a boundary flux strongly affects the deposition rate at the corresponding boundary node, as the local oversaturation/undersaturation governs the magnitude of the deposition/sublimation that is necessary to maintain vapor saturation. It follows that the prescription of proper vapor fluxes at boundaries is an integral part of the physical problem and should ideally arise from a vapor budget at interfaces, similarly to the standard surface energy budget. On the other hand, fixing c to an arbitrary value at boundaries (Dirichlet type BCs), as done by Simson et al. (2021), implicitly implies the prescription of a water vapor mass flux that adjusts itself so that saturation is maintained with a fixed contribution of deposition/sublimation. This boundary vapor flux then corresponds to the residuals of the vapor mass budget equation. This topic is further explored in Sect. 4.6.

## 3.3 A mass-conservative Lagrangian settlement scheme

### 3.3.1 General considerations

Our settlement scheme is based on the method of characteristics as presented in Simson et al. (2021) with some corrections to guarantee ice mass conservation as well as a more explicit formulation of the element-wise nature of this conservation. The general idea of this approach is that the mesh nodes should move at the velocity of the ice matrix as the latter settles so that all ice mass fluxes related to settlement wipe out in the working frame, thereby resorting to a Lagrangian coordinate system. Concretely, this allows for the elimination of the advection term $v \cdot \nabla \Phi_i$ from the continuity equation, which is then in one dimension:

$$ \partial_t \Phi_i + \Phi_i \partial_z v_z = \frac{c}{\rho_i}, \qquad (14) $$

where $\partial_z v_z$ is the divergence of the settling velocity which, in one dimension, directly corresponds to the vertical strain rate $\dot{\epsilon}_{zz}$. The vertical strain rate relates to the vertical stress $\sigma_{zz}$ through the constitutive law of snow, which we take as a simple linear viscous law:

$$ \dot{\epsilon}_{zz} = \frac{1}{\eta} \sigma_{zz}, \qquad (15) $$

where $\eta$ is the effective viscosity of snow and $\sigma_{zz}$ is the vertical stress. The effective viscosity is an important and poorly constrained snow property that depends, among others, on microstructure and temperature (e.g., Wiese and Schneebeli, 2017). Here, we use the viscosity parameterization of Vionnet et al. (2012) where the viscosity increases exponentially with density and decreases exponentially with temperature according to Eq. (A7). Note that Simson et al. (2021) also considered the case of a nonlinear Glen's flow law as well as

the case of a constant viscosity. As the goal of this study is not an assessment of the model sensitivity to parameterization choices but rather to the details of the numerical treatment of the equations, we do not consider these cases here. Neglecting the contribution of air in the effective density of snow, the momentum conservation equation relates the distribution of the vertical stress to that of the ice volume fraction as follows:

$$ \partial_z \sigma_{zz} = \rho_i g \Phi_i(z), \qquad (16) $$

where $g$ is gravity.

As underlined by Simson et al. (2021), Eq. (14) and its Lagrangian perspective corresponds to the strategy employed in all available detailed snowpack models to represent settlement, although it is usually not explicitly stated. Concretely, the constitutive law of snow is used to relate the stress supported by a layer, calculated as the cumulative weight of all overlying layers, to its total deformation. This total deformation is then used to either directly update the layer thickness defined as a state variable (e.g., Jordan, 1991; Vionnet et al., 2012) or to update the mesh coordinates using the fact that the node at the soil-snow interface does not move (e.g., Bader and Weilenmann, 1992; Bartelt and Lehning, 2002). The effective density, defined as a layer property, is then simply updated so that, in the absence of phase change, the mass contained within each layer remains the same before and after settlement. The layer-based nature of this scheme is obvious, in that conservation is not fulfilled locally but rather on average over finite space intervals referred to as control volumes (e.g., Jordan, 1991), elements (e.g., this study, Bartelt and Lehning, 2002), or layers (e.g., Bader and Weilenmann, 1992; Bartelt and Lehning, 2002; Vionnet et al., 2012), depending on the numerical method chosen to solve the other PDEs of the model. In contrast, as Eq. (14) expresses the local ice mass conservation, one may assume that its actual solution would enable a switch from this traditional layer-based vision to a more continuous description of the snowpack. We stress that this impression relies on a confusion between numerical and physical layers. Indeed, although all physical quantities involved in Eqs. (14)-(16) are continuous functions of z and can be calculated anywhere in the continuous domain, as soon as we go through the necessary numerical discretization step, the continuous vision breaks and we step back to a discrete description in which the ice phase is conserved on average over finite space intervals that can be seen as numerical layers. The settlement scheme then consists of two numerical operations that must be done in parallel and in a consistent way to ensure that conservative properties of Eq. (14) are preserved: (i) the update of the mesh node positions so that the variation in the mass contained in each element over one time step is entirely due to phase change and not at all to settlement and (ii) the update of the ice volume fraction defined as an element-wise state variable accounting for the possible source term due to phase change.
### 3.3.2 Update of mesh node positions

The displacement $\Delta z_{k+1}$ of node $k+1$ between $t$ and $t + \Delta t$ corresponds to the displacement $\Delta z_k$ of the underlying node $k$ to which the total deformation of the space interval between $k$ and $k+1$ over $\Delta t$ must be added. This condition can be written as follows:

$$ \Delta z_{k+1} = \Delta z_k + \int_{z_k^n}^{z_{k+1}^n} \dot{\epsilon}_{zz}^{n+1}(z) dz \Delta t, \quad (17) $$

where $\dot{\epsilon}_{zz}^{n+1}$ is the vertical strain rate that, in one dimension, directly corresponds to the divergence of the vertical velocity field $\partial_z v_z$. By exploiting the fact that the node located at the soil-snow interface is immobile, we can trace the displacement of all nodes of the domain through Eq. (17). A difficulty arises in the numerical treatment of the integral term of Eq. (17) after spatial discretization. Indeed, the numerical integration requires an assumption regarding the spatial variation in the integrand in between nodes. While this assumption is made a priori via the choice of the shape functions in the FEM, only the nodal values of the fields are considered to constitute the solution in the finite-difference method employed by Simson et al. (2021), without any explicit reference to how these fields should vary between nodes (Patankar, 1980). In their mesh deformation procedure, Simson et al. (2021) do not directly use the strain rate $\dot{\epsilon}_{zz}$ to evaluate the total deformation of the space intervals in between nodes; rather, they perform a numerical integration to recover the settling velocity $v_z$ and use the latter to move the nodes. A careful inspection of Eq. (17) of Simson et al. (2021) reveals that this numerical integration implicitly assumes that the strain rate calculated at node $k$ from the stress and viscosity evaluated at node $k$ actually applies to the whole space interval between the node $k$ and the node $k-1$ (assuming that the $\Delta z_j^n = z_{j+1}^n - z_j^n$ stated in Eq. 17 of Simson et al., 2021, actually means $\Delta z_j^n = z_j^n - z_{j-1}^n$, which makes more sense and is in accordance with their published code). Similarly, the numerical integration of the momentum conservation equation yielding the nodal stresses from the distribution of $\Phi_i$ as expressed in Eq. (18) of Simson et al. (2021) implicitly assumes that the value of $\Phi_i$ stored at node $k$ actually applies to the whole space interval between nodes $k$ and $k+1$. In other words, in the approach of Simson et al. (2021), the distribution of $\Phi_i$ is piecewise constant over numerical layers, but the information is only assigned to the bottom node of the considered numerical layer. As detailed in Appendix D1, the fact that $\Phi_i$ at node $k$ relates to the snow mass contained between nodes $k$ and $k+1$, whereas the strain rate calculated at node $k$ is used to deform the space interval between nodes $k$ and $k-1$, leads to an inconsistency that hampers mass conservation.

In contrast, in our approach, $\Phi_i$ is defined in an element-wise manner, whereas the stress is a nodal variable. The latter

is calculated at each node $k$ as the cumulative weight of all overlying elements:

$$ \sigma_{zz,k}^{n+1} = \sum_{j=k}^{N_z-1} g \Phi_{i,j+1/2}^{n} \rho_i (z_{j+1}^n - z_j^n). \quad (18) $$

This equation is formally equivalent to Eq. (18) of Simson et al. (2021) except that $\Phi_{i,j+1/2}$ explicitly corresponds to an averaged quantity applying to the whole element $j+1/2$ between nodes $j$ and $j+1$. The motion of all nodes can then be determined directly by solving Eq. (17) where the integral term is treated through the Gaussian quadratures.

### 3.3.3 Update of ice volume fraction

It can be shown that conservation of the ice mass is guaranteed only if the temporal discretization of Eq. (14) is based on an implicit numerical scheme (Appendix E). Therefore, in our approach, we replace the first-order Euler explicit temporal discretization given in Eq. (15) of Simson et al. (2021) with the following:

$$ \Phi_{i,k+1/2}^{n+1} = \frac{\Phi_{i,k+1/2}^{n} + \Delta t \frac{c_{k+1/2}^{n+1}}{\rho_i}}{1 + \Delta t \epsilon_{zz,k+1/2}^{n+1}}, \quad (19) $$

where the mean strain rate $\dot{\epsilon}_{zz,k+1/2}^{n+1}$ and mean deposition rate $c_{k+1/2}^{n+1}$ of element $k+1/2$ are calculated as the integral of the corresponding field over the element using Gaussian quadratures divided by the length of the considered element before deformation.

## 4 Numerical simulations

In this section, we compare the capabilities of the various numerical treatments introduced above to provide solutions to the coupled problem of heat transport, vapor transport, and settlement in snow that are stable, accurate, and respect energy and mass conservation. To this end, we use some of the experiments proposed by Schürholt et al. (2022b) and Simson et al. (2021) as numerical benchmarks. The main features of the numerical setup and model configurations of all experiments described in the following are summarized in Table 2. All simulations are run with the parameterization (Eq. A6) for the saturation vapor density and with the Calonne parameterizations of effective parameters (given by Eqs. A1–A3), unless stated otherwise. Note that the Hansen parameterizations of effective parameters given by Eqs. (A2)–(A4) have also been implemented. In what follows, we use the absolute root-mean-square deviations (RMSDs) as a metric when comparing two solution fields produced with two different implementations. For interpretation, it is important to relate these RMSDs to the typical orders of magnitude of the considered fields shown in respective figures.
<table><caption>**Table 2.** Summary of numerical setups and model configurations. Scenarios 2 and 3 are taken from Schürholt et al. (2022b). Cases 6 and 7 are taken from Simson et al. (2021). Dirichlet BCs for *T* are *T*₀ = 273 K and *T*<sub>h</sub> = 253 K. Dirichlet BCs for *ρ*<sub>v</sub> are *ρ*<sub>v,0</sub> = *ρ*<sub>v</sub><sup>eq</sup>(*T*₀) and *ρ*<sub>v,h</sub> = *ρ*<sub>v</sub><sup>eq</sup>(*T*<sub>h</sub>).</caption><thead><tr><th>Section</th><th>Setup</th><th>N<sub>z</sub></th><th>Δt (min)</th><th>BC for T</th><th>BC for ρ<sub>v</sub></th><th>Configuration</th><th>Mass matrix<sup>a</sup></th><th>α</th><th>Dep.<sup>a</sup></th><th>Settl.<sup>a</sup></th></tr></thead><tbody><tr><td rowspan="4">4.1</td><td rowspan="4">Scenario 2</td><td rowspan="4">201</td><td rowspan="4">15</td><td rowspan="4">Dirichlet</td><td rowspan="4">Dirichlet</td><td>FEniCS<sup>b</sup></td><td>I/N.L.</td><td rowspan="4">5 × 10<sup>-3</sup></td><td rowspan="4">On</td><td rowspan="4">Off</td></tr><tr><td>CC_3DOF</td><td>I/N. L.</td></tr><tr><td>CC_3DOF</td><td>P/N. L.</td></tr><tr><td>CC_3DOF</td><td>P/L.</td></tr><tr><td rowspan="6">4.2</td><td rowspan="6">Scenario 2</td><td rowspan="6">201</td><td rowspan="6">15, 5</td><td>Dirichlet</td><td>Dirichlet</td><td>CC_2DOF</td><td>P/L.</td><td>5 × 10<sup>-3</sup></td><td rowspan="6">On</td><td rowspan="6">Off</td></tr><tr><td>Dirichlet</td><td>Dirichlet</td><td>CC_3DOF</td><td>P/L.</td><td>5 × 10<sup>-3</sup></td></tr><tr><td>Dirichlet</td><td>Dirichlet</td><td>CD_PC</td><td>P/L.</td><td>5 × 10<sup>-3</sup></td></tr><tr><td>Dirichlet</td><td>Dirichlet</td><td>CD_FD</td><td>P/L.</td><td>5 × 10<sup>-3</sup></td></tr><tr><td>Dirichlet</td><td>No flux</td><td>H_MF</td><td>-/L.</td><td>-</td></tr><tr><td>Dirichlet</td><td>No flux</td><td>H_TF</td><td>P/L.</td><td>-</td></tr><tr><td rowspan="2">4.3</td><td rowspan="2">Scenario 2</td><td rowspan="2">201</td><td rowspan="2">15</td><td>Dirichlet</td><td>No flux</td><td>CC_3DOF</td><td>P/L.</td><td>0 → 1</td><td rowspan="2">On</td><td rowspan="2">Off</td></tr><tr><td>Dirichlet</td><td>No flux</td><td>H_MF</td><td>-/L.</td><td>-</td></tr><tr><td rowspan="4">4.4</td><td rowspan="4">Scenario 2</td><td rowspan="4">201</td><td rowspan="4">15, 5</td><td rowspan="4">No flux</td><td rowspan="4">No flux</td><td>CC_3DOF</td><td>P/L.</td><td>5 × 10<sup>-3</sup></td><td rowspan="4">On, Off</td><td rowspan="4">Off</td></tr><tr><td>CD_PC</td><td>P/L.</td><td>5 × 10<sup>-3</sup></td></tr><tr><td>H_MF</td><td>-/L.</td><td>-</td></tr><tr><td>H_TF</td><td>P/L.</td><td>-</td></tr><tr><td>4.5</td><td>Case 6</td><td>11, 51, 101</td><td>15</td><td>-</td><td>-</td><td>Simson code<sup>c</sup><br/>Our code</td><td>-</td><td>-</td><td>Off</td><td>On</td></tr><tr><td rowspan="3">4.6</td><td rowspan="3">Case 7</td><td rowspan="3">101</td><td>Adaptive</td><td>Dirichlet</td><td>c = 0</td><td>Simson code<sup>c</sup></td><td>P/-</td><td>-</td><td rowspan="3">On</td><td rowspan="3">On</td></tr><tr><td>15</td><td>Dirichlet</td><td>c = 0</td><td>H_MF</td><td>-/L.</td><td>-</td></tr><tr><td>15</td><td>Dirichlet</td><td>No flux</td><td>H_MF</td><td>-/L.</td><td>-</td></tr><tr><td rowspan="2">4.7</td><td rowspan="2">Scenario 3</td><td rowspan="2">1001</td><td rowspan="2">1</td><td rowspan="2">Dirichlet</td><td rowspan="2">Dirichlet</td><td>FEniCS<sup>b</sup></td><td>I/N.L.</td><td rowspan="2">5 × 10<sup>-3</sup></td><td rowspan="2">On</td><td rowspan="2">Off</td></tr><tr><td>CC_3DOF</td><td>I/N. L.</td></tr><tr><td rowspan="2">4.8</td><td rowspan="2">Scenario 3</td><td rowspan="2">1001</td><td rowspan="2">1</td><td rowspan="2">Dirichlet</td><td rowspan="2">Dirichlet</td><td>CC_3DOF</td><td>P/L.</td><td rowspan="2">5 × 10<sup>-3</sup></td><td rowspan="2">On</td><td>On</td></tr><tr><td>CC_3DOF</td><td>P/L.</td><td>Off</td></tr></tbody></table>

## 4.1 On the form of the mass matrices

In this part, we reproduce Scenario 2 presented in Schürholt et al. (2022b) to illustrate the importance of dealing with proper mass matrices. Concretely, we consider a 1 m thick snowpack with a piecewise linear initial density profile mimicking a stratified snowpack containing a crust as well as an ice layer at the bottom (Eq. 16 of Schürholt et al., 2022b). All BCs are Dirichlet-type conditions and are constant over time: the bottom and top temperatures are fixed to *T*₀ = 273 K and *T*<sub>h</sub> = 253 K, respectively, while the bottom and top water vapor densities are set assuming saturation at boundaries, i.e., *ρ*<sub>v,0</sub> = *ρ*<sub>v</sub><sup>eq</sup>(*T*₀) and *ρ*<sub>v,h</sub> = *ρ*<sub>v</sub><sup>eq</sup>(*T*<sub>h</sub>). The initial temperature profile is linear between boundary values, and the initial water vapor density profile is deduced from the latter assuming that the water vapor density is initially saturated (solid black lines in Fig. 3). To facilitate a comparison with results presented by Schürholt et al. (2022b), all simulations described in this part are run with settlement deactivated and deposition feedback on Φᵢ activated. The time step is set to Δ*t* = 15 min. The mesh is uniform and contains 200

elements. For the simulations with the Calonne model, the sticking coefficient is set to a default value of α = 5 × 10⁻³ so that our *s*α*ν*<sub>kin</sub> is comparable to the *ρ*ᵢs/β*ρ*<sub>v</sub><sup>eq</sup> adopted by Schürholt et al. (2022b).

Figure 3 shows the vertical profiles of temperature, vapor density, and deposition rate produced after 38 h of simulation by FEniCS and by different versions of CC_3DOF. The FEniCS run is based on the code of Schürholt et al. (2022a), with slight adaptations. In particular, the original Crank–Nicolson time scheme adopted by Schürholt et al. (2022b) is replaced by a fully implicit scheme. Despite slight differences in the numerical implementation (e.g., Φᵢ is a nodal variable in FEniCS and an element-wise variable in our model, variable-dependent parameters are evaluated at nodes in FEniCS and directly at Gaussian points in our model, and the expressions of the source terms are not strictly equivalent in FEniCS and in our model), the solution fields produced by FEniCS are very close to those obtained with CC_3DOF when the mass matrix has the improper form and no lumping is performed (superimposed solid green and dotted red lines). This can be quantified by calculating the RMSDs between the solutions

<footnote><sup>a</sup> The following abbreviations are used in the table: I. – improper (Sect. 3.2); P. – proper (Sect. 3.2); N.L. – not lumped (Sect. 3.2); L. – lumped (Sect. 3.2); Dep. – deposition; and Settl. – settlement. <sup>b</sup> Code of Schürholt et al. (2022a). <sup>c</sup> Code of Simson and Kowalski (2021).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of temperature (a), water vapor density (b), and deposition rate (c) fields produced by FEniCS, CC_3DOF without lumping and with the improper mass matrix, CC_3DOF without lumping but with the proper mass matrix, and CC_3DOF with lumping and the proper mass matrix after 38 h of simulation for Scenario 2 of Schürholt et al. (2022b). The dashed blue and solid orange lines are superimposed in panels (a) and (b). The solid black lines correspond to the initial conditions.</figcaption></figure>

obtained with the two approaches: it is $3.7 \times 10^{-4}$ K for $T$ and $4.4 \times 10^{-8}$ kg m$^{-3}$ for $\rho_v$. A more noticeable difference is found with respect to the amplitude of the oscillations occurring on the field of $c$ at the top boundary of the domain. Note that oscillations are also observed at the bottom boundary of the domain as well as at the boundaries of the deposition/sublimation peaks, but they are of similar amplitude in both approaches, except at the very bottom node. Thus, the RMSD between the fields of $c$ produced by the two approaches is $3.7 \times 10^{-6}$ kg m$^{-3}$ s$^{-1}$ when the whole domain is considered, whereas it is reduced to $2.4 \times 10^{-8}$ kg m$^{-3}$ s$^{-1}$ when the first bottom node and four top nodes are omitted.

In contrast, adopting the proper form of the mass matrix significantly changes the $T$ and $\rho_v$ fields in places where peaks in sublimation/deposition are observed (dashed blue lines hidden below the solid orange lines in Fig. 3a and b): the RMSDs between the solution fields of $T$ and $\rho_v$ obtained

with CC_3DOF using the improper mass matrix and those obtained with CC_3DOF using the proper mass matrix are increased to $0.42$ K and $6.5 \times 10^{-5}$ kg m$^{-3}$, respectively. The sublimation/deposition peaks also become sharper (dashed blue line in Fig. 3c). The oscillations are mostly reproduced, with larger amplitudes at deposition/sublimation peaks and smaller amplitudes at domain boundaries, especially at the bottom boundary. These spurious oscillations are related to the violation of the DMP. Note that oscillations of the same nature also occur when FEniCS is run with a Crank–Nicolson time scheme (not shown). As expected, these oscillations vanish when the mass matrix is lumped (solid orange line in Fig. 3).
## 4.2 On the different solver options

Here, we use the same numerical setup as above in order to compare the relative performance of the various numerical treatments of the Calonne system that are summarized in Fig. 1b. All simulations are run with the proper form of the mass matrices and with lumping activated. In addition to the time step of $\Delta t = 15$ min, all simulations are also run with a time step reduced to $\Delta t = 5$ min. A first important result is that the boundary values taken by the field of $c$ depart strongly from their distribution in the interior of the domain and that they are highly sensitive to the chosen numerical approach (Fig. 4c). This topic is treated in detail in Sect. 4.6. Secondly, it turns out that the two strategies investigated for the coupled solution of the Calonne system, i.e., CC_2DOF and CC_3DOF, produce solution fields that are very close. After 24 h of simulation with $\Delta t = 15$ min, the RMSDs for the $T$, $\rho_v$, and $c$ fields produced with the two methods amount to $1.8 \times 10^{-4}$ K, $1.6 \times 10^{-8}$ kg m$^{-3}$, and $1.2 \times 10^{-6}$ kg m$^{-3}$ s$^{-1}$, respectively. The higher (relative) deviation in $c$ is mostly due to the values obtained at the two boundary nodes, which are slightly sensitive to the chosen strategy: if the very top and bottom nodes are omitted, the RMSD for the $c$ field is reduced to $3.7 \times 10^{-9}$ kg m$^{-3}$ s$^{-1}$. In addition, the solution fields produced by CC_2DOF and CC_3DOF show very little sensitivity to the time step size (superimposed green and blue lines in Fig. 4) with similar RMSDs as above for all three fields when $\Delta t = 5$ min. In contrast, the decoupled solution of the Calonne system leads to significantly different behaviors. First, the CD FD approach is highly unstable, with the $\rho_v$ and $T$ fields rapidly diverging towards unrealistic values. This continues to be the case even when the time step size is decreased down to $\Delta t = 1$ s. In contrast, the CD PC approach gives steady solutions that are mostly independent of the time step. These steady solutions are not very different from those obtained with the coupled approaches (solid red and orange lines in Fig. 4). Concretely, after 24 h, the RMSDs between the solution fields produced with CC_3DOF and CD PC are $4.2 \times 10^{-2}$ K, $4.8 \times 10^{-6}$ kg m$^{-3}$, and $4.3 \times 10^{-5}$ kg m$^{-3}$ s$^{-1}$ ($2.4 \times 10^{-6}$ kg m$^{-3}$ s$^{-1}$ when the top and bottom nodes are omitted) for the $T$, $\rho_v$, and $c$ fields, respectively. However, the transient solution fields produced by CD PC show oscillations (even when the mass matrices are lumped) and high sensitivity to the time step size (dotted and dashed red and orange lines in Fig. 4). For example, after 2 h of simulation, RMSDs between CD PC solutions obtained with $\Delta t = 15$ min and with $\Delta t = 5$ min amount to $1.3$ K, $1.7 \times 10^{-4}$ kg m$^{-3}$, and $7.2 \times 10^{-4}$ kg m$^{-3}$ s$^{-1}$ ($2.4 \times 10^{-5}$ kg m$^{-3}$ s$^{-1}$ when the top and bottom nodes are omitted) for the $T$, $\rho_v$, and $c$ fields, respectively. This behavior has to be compared to the transient solutions produced by CC_3DOF/CC_2DOF (blue lines in Fig. 4) that are smooth and smoothly converge to the steady solutions.

The much higher stability of CD PC compared with CD FD is due to the proper treatment of the reaction term in the water vapor mass balance equation, which acts as an attractor of $\rho_v$ toward $\rho_v^{eq}$. Nevertheless, the results presented in this paragraph underline the major importance of having a coupled solution of the coupled heat and water vapor diffusion equations. Although the CD PC approach gives steady solutions that are close to those obtained with the coupled approaches, the accuracy of the transient responses is essential, as the external forcings are constantly evolving in time for the vast majority of real-world configurations.

The same experiment has been run with the H_MF and H_TF configurations (not shown). The two strategies produce solution fields of $T$ (and, consequently, of $\rho_v$ and $c$, as both fields are diagnosed from the obtained $T$ field) that are very close to each other and that show very little sensitivity to the time step size: the RMSD between both strategies after 24 h of simulation with $\Delta t = 5$ min ($\Delta t = 15$ min) is $6.2 \times 10^{-3}$ K ($6.2 \times 10^{-3}$ K), $5.6 \times 10^{-7}$ kg m$^{-3}$ ($5.5 \times 10^{-7}$ kg m$^{-3}$), and $7.3 \times 10^{-9}$ kg m$^{-3}$ s$^{-1}$ ($7.3 \times 10^{-9}$ kg m$^{-3}$ s$^{-1}$) for the $T$, $\rho_v$, and $c$ fields, respectively. As we shall see in Sect. 4.4, only the H_MF approach shows perfect energy conservation.

## 4.3 Comparing the Hansen and Calonne systems for high $\alpha$

Figure 5 shows the solutions produced after 38 h with our model using the following configurations: (1) H_MF and (2) CC_3DOF for various values of the sticking coefficient $\alpha$. All simulations are run with the proper form of the mass matrices from the numerical setup introduced above, except for the BCs for $\rho_v$ which are forced to no flux in all cases (Table 2) so that both configurations are easily comparable (see Sect. 4.6). For both configurations, effective parameters are computed from the Calonne parameterizations (Appendix A). We recall that the Calonne system has been derived via the two-scale expansion method that is valid for low reaction rates only, i.e., $\alpha \sim 10^{-3}$ or less (Sect. 2.1). Therefore, neither the Calonne system nor the Calonne parameterization of effective parameters is expected to be valid for higher values of $\alpha$. However, in a similar vein to Schürholt et al. (2022b), our goal here is simply to compare the mathematical behaviors of the Hansen and Calonne systems when they are run with the same parameterizations, independently of the physical soundness of obtained results.

The solution fields produced by CC_3DOF progressively converge to those produced by H_MF as $\alpha$ increases and a higher reaction rate forces the vapor to saturation. This behavior was expected. Indeed, if effective parameters are calculated the same way and BCs are identical, both systems are formally equivalent when water vapor density is assumed to always be at saturation (Sect. 2.3). To quantify this behavior, Table 3 gathers all RMSDs between solutions produced
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of temperature (a), water vapor density (b), and deposition rate (c) fields produced by CC_3DOF and CD_PC with $\Delta t = 5$ min and $\Delta t = 15$ min. Solutions are shown after 2, 4, and 24 h of simulation. A zoomed-in view of the deposition peaks is included in panel (c). Note that the green and blue lines are almost superimposed in all panels.</figcaption></figure>

at the end of the simulation by (1) H_MF and (2) CC_3DOF for all tested values of $\alpha$.

The difference in $\rho_v$ and $c$ between H_MF and CC_3DOF becomes less than 1% when $\alpha$ becomes higher than $10^{-6}$ and $10^{-3}$, respectively. Below these values of $\alpha$, the difference varies from a few percent to ~100% as $\alpha$ tends to 0. In contrast, $T$ is less affected by the value of $\alpha$. Even for the lowest values tested, the differences between the fields of $T$ produced by H_MF and CC_3DOF are below 0.1% (all lines are superimposed in Fig. 5a).

All together, these results seem to contradict those presented in Fig. 2 of Schürholt et al. (2022b): in their case, the solution fields (especially c) produced by the Hansen and Calonne models differ significantly, even when both models rely on the Calonne parameterizations for effective parameters. We think that this is due to the improper treatment of the mass matrices in their FEniCS implementation. This conclusion is supported by the fact that “existing continuum-mechanical models derived through homogenization (i.e., Calonne) or mixture theory (i.e., Hansen) yield sim-
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of temperature (a), water vapor density (b), and deposition rate (c) fields produced after 38 h of simulation with H_MF and with CC_3DOF for various values of $\alpha$. Note that the case $\alpha = 1$ is not represented, as it would not be distinguishable from the case $\alpha = 10^{-1}$.</figcaption></figure>

<table><caption>**Table 3.** RMSDs between solutions obtained with H_MF and with CC_3DOF for the various values of $\alpha$.</caption><thead><tr><th>RMSD for</th><th><i>T</i> (K)</th><th><i>ρ</i><sub>v</sub> (kg m<sup>-3</sup>)</th><th><i>c</i> (kg m<sup>-3</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>α = 0</td><td>4.3 × 10<sup>-2</sup></td><td>1.1 × 10<sup>-3</sup></td><td>1.5 × 10<sup>-6</sup></td></tr><tr><td>α = 10<sup>-8</sup></td><td>2.1 × 10<sup>-2</sup></td><td>6.8 × 10<sup>-5</sup></td><td>1.4 × 10<sup>-6</sup></td></tr><tr><td>α = 10<sup>-6</sup></td><td>8.1 × 10<sup>-3</sup></td><td>1.6 × 10<sup>-6</sup></td><td>9.3 × 10<sup>-7</sup></td></tr><tr><td>α = 10<sup>-4</sup></td><td>1.1 × 10<sup>-2</sup></td><td>1.0 × 10<sup>-6</sup></td><td>6.7 × 10<sup>-8</sup></td></tr><tr><td>α = 10<sup>-1</sup></td><td>1.1 × 10<sup>-2</sup></td><td>1.0 × 10<sup>-6</sup></td><td>9.4 × 10<sup>-9</sup></td></tr><tr><td>α = 1</td><td>1.1 × 10<sup>-2</sup></td><td>1.0 × 10<sup>-6</sup></td><td>9.4 × 10<sup>-9</sup></td></tr></tbody></table>
<table><caption>Table 4. Total energy leakage (J m⁻²) after 5 d of simulation.</caption><thead><tr><th rowspan="2">Deposition</th><th colspan="2">Off</th><th colspan="2">On</th></tr><tr><th>5 min</th><th>15 min</th><th>5 min</th><th>15 min</th></tr></thead><tbody><tr><td>Δt =</td><td></td><td></td><td></td><td></td></tr><tr><td>CC_3DOF</td><td>0</td><td>0</td><td>-296.3</td><td>-295.0</td></tr><tr><td>CC_PC</td><td>-2.8</td><td>-6.4</td><td>-299.1</td><td>-301.4</td></tr><tr><td>H_MF</td><td>0</td><td>0</td><td>-296.3</td><td>-295.0</td></tr><tr><td>H_TF</td><td>0.14</td><td>0.02</td><td>-296.2</td><td>-295.0</td></tr></tbody></table>

ilar results for homogeneous snowpacks of constant density" (Schürholt et al., 2022b), as demonstrated from their Scenario 1. Indeed, when density is uniform, the multiplying factors of the accumulation terms are independent of space and can be arbitrarily affected directly to the mass matrix or their inverses to the stiffness matrix and force vector.

## 4.4 On energy conservation

In order to evaluate how the different numerical strategies behave with respect to energy conservation, we run a slightly modified version of the numerical setup introduced in Sect. 4.1: the original Dirichlet BCs are replaced by no-flux BCs for both heat and vapor. Because settlement is deactivated (no vapor expelled from the snowpack to the atmosphere), the system is closed and all recorded energy leakage should be considered to be a numerical artifact. We run the experiment for 5 d, testing all strategies of the three options presented in Sect. 3.2 and summarized in Fig. 1b, except CD FD which yields unrealistic results even for a time step as low as $\Delta t = 1$ s. All experiments are run twice: once with $\Delta t = 15$ min and another time with $\Delta t = 5$ min. Then, the feedback of deposition on $\Phi_i$ is deactivated and the procedure is repeated. The obtained energy leakage is represented as a function of time in Fig. 6. Again, results produced with CC_3DOF are very close to those obtained with CC_2DOF, and only the former are reported. The total energy leakage at the end of the simulation is summarized in Table 4.

As expected, decoupling the deposition-related evolution of $\Phi_i$ (source term of the ice mass conservation equation) from the heat/vapor transfer equation induces an artificial energy loss of ~300 J m⁻² for all considered cases at the end of the simulation. This energy loss remains very limited compared with the total energy contained in the system but could become noticeable for simulations on seasonal timescales and/or for configurations associated with stronger deposition rates. On the contrary, when deposition is omitted, energy leakage become negligible. Nevertheless, we stress that, contrary to the CC_2DOF, CC_3DOF, and H_MF strategies that are rigorously conservative, the CD_PC and H_TF strategies induce tiny but nonzero energy leakage, in line with the demonstration of Tubini et al. (2021). This illustrates that a rigorously energy-conservative solution of the problem of heat conduction and water vapor diffusion in

snow implies solving the heat diffusion, water vapor diffusion, and ice mass conservation equations in a coupled way for the Calonne system and opting for the mixed form of the heat equation for the Hansen system.

## 4.5 On mass conservation

In this part, we reproduce the numerical setup corresponding to the Case 6 designed by Simson et al. (2021): we consider a snowpack with an initial thickness of 50 cm split into two equally thick snow layers of uniform initial densities, i.e., 150 kg m⁻³ for the lower one and 75 kg m⁻³ for the top one. All simulations are run considering only the settlement process. Our goal is to illustrate how the modifications made to the settlement scheme of Simson et al. (2021) affect the mass conservation and settlement rates. For this, we run the code of Simson and Kowalski (2021) with three bug fixes that are described in Appendix D. In particular, the strain rate governing the deformation of the whole space interval comprised between nodes $k$ and $k+1$ is calculated from the stress and viscosity evaluated at node $k$ (and not at node $k+1$ as in the published version of the code; Appendix D1). Simulations are run for 20 d, with $\Delta t = 15$ min, and with three initially uniform meshes containing 11, 51, and 101 nodes, respectively. For all runs, viscosity is calculated from Eq. (A7) with a fixed temperature set to $T = 263$ K.

We first run the simulations using the original explicit time discretization scheme implemented by Simson and Kowalski (2021) to update $\Phi_i$. As illustrated in Fig. 7a, mass is not conserved in this case: there is an artificial mass loss that gets higher for coarser meshes. In the present case, this mass loss is in the order of 20 g after 20 d of simulation for a total initial mass of 56.25 kg. Although this mass loss tends to stabilize after a few days as viscosity increases exponentially with increasing $\Phi_i$, it will become more significant when integrated over a full season characterized by regular snowfalls of low initial densities. As soon as the explicit time integration scheme is replaced by an implicit one, i.e., following Eq. (19) and Appendix D4, mass is perfectly conserved regardless of the mesh size. Note that this is the case when running both a corrected version of the code published by Simson and Kowalski (2021), in which the explicit settlement scheme has been replaced by an implicit one, and our model (green line in Fig. 7a). In Fig. 7b, we compare the settlement rates obtained with the code of Simson and Kowalski (2021) using the implicit time integration scheme to the one obtained with our code for an 11-node mesh. Using the approach of Simson and Kowalski (2021) corrected following Appendix D1 to ensure conservation of mass, the deformation of the whole numerical layer between $k$ and $k+1$ is implicitly calculated from the stress evaluated at the bottom node of the layer, where it is maximal. In contrast, we assess the deformation occurring between nodes $k$ and $k+1$ by integrating the ratio between the stress and viscosity fields along the element $k+1/2$. As a consequence, the former treatment
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Energy leakage as a function of simulation time for the Calonne (a) and Hansen (b) systems, with and without deposition. For the Calonne system, we test the CC_3DOF and CD_PC strategies. For the Hansen system, we test the H_MF and H_TF strategies. Each simulation is run with $\Delta t = 5$ min and $\Delta t = 15$ min. Note that orange and red lines are almost superimposed in panel (b). Both panels contain a zoomed-in view of the region around $\Delta E_{\Omega} = 0$ J m$^{-2}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Panel (a) shows the evolution of the total mass over time when $\Phi_i$ is updated using the original explicit time discretization scheme implemented by Simson and Kowalski (2021) for meshes containing 11, 51, and 101 nodes. As soon as the time discretization scheme becomes implicit, mass is conserved independently of mesh refinement (green line). Panel (b) shows the absolute settling velocity profile obtained on the 11-node mesh with our model and with the corrected code of Simson et al. (2021) using an implicit time discretization scheme after 1, 2, 5, and 20 d of simulation.</figcaption></figure>
tends to slightly overestimate the settlement rate compared with ours. This sensitivity is stronger at the beginning of the simulation when settlement rates are higher, and it tends to decrease over time. As expected, when the mesh resolution is increased, the two approaches converge to the same settlement rates (not shown).

## 4.6 On the boundary conditions for water vapor density

In this part, we want to highlight the sensitivity of the solution fields to the treatment of water vapor at boundaries. To this end, we reproduce Case 7 proposed by Simson et al. (2021). Case 7 corresponds to the same numerical setup as that of Case 6 introduced above, but the heat and water vapor diffusion processes are included. Furthermore, the feedback of temperature on viscosity is taken into account via Eq. (A7).

Figure 8 shows the $c$ and $\Phi_i$ fields obtained after 1, 2, 5, and 10 d of simulation on a 101-node mesh when running the code published by Simson and Kowalski (2021) with an implicit settlement scheme and with the bug fixes described in Appendix D (black lines in Fig. 8a). The values of $c$ are always zero at the boundaries because they are forced as such in the model by the authors. Despite slight differences in the numerical implementation (e.g., element-wise vs. nodal nature of $\Phi_i$, treatment of the strain rate, and mixed form vs. $T$ form for the Hansen system), running our model in its H_MF configuration gives solution fields that are consistent with those produced by the code of Simson and Kowalski (2021) if we also force the boundary values of $c$ to zero (orange lines in Fig. 8a). Concretely, the RMSDs for $c$ vary between $2.8 \times 10^{-8}$ and $3.4 \times 10^{-8}$ kg m$^{-3}$ s$^{-1}$, corresponding to differences in the order of ~1%.

As stated in Sect. 3.2, because the deposition rate is derived as the closure of the water vapor mass balance, imposing a vanishing deposition rate at a boundary is equivalent to imposing a boundary vapor flux that adjusts itself so that local disequilibrium between $\rho_v$ and $\rho_v^{eq}$ are entirely compensated for without any requirements for phase change. This choice is arbitrary and cannot be used as the generic BC for snowpack models. Ideally, the BCs for the vapor mass budget should be derived through a vapor budget at interfaces (akin to the surface energy budget routinely used as BCs for the energy equation). The importance of the BCs for vapor can be illustrated by running the same H_MF simulation as above but removing all constraints on $c$ and assuming no vapor fluxes at boundaries. The resulting field of $c$ is then characterized by peak values at the two boundary nodes, which depart significantly from the distribution of $c$ in the interior of the domain (blue lines in Fig. 8a). These high deposition rates at domain boundaries have direct implications for the field of $\Phi_i$: a significant part of the ice mass of the bottom element sublimates and is transported upwards. This leads to a situation in which the bottom element is less dense than

the one immediately above, a tendency that is exacerbated over time (blue lines in Fig. 8b). This behavior creates strong local density gradients that seem to trigger self-amplifying oscillations at the bottom boundary: a sublimation peak at the bottom node is followed by a deposition peak at the node right above. This oscillatory pattern propagates inwards over a number of nodes that grows in time. On the top, strong mass gain occurs associated with a high deposition peak over the top element but does not trigger oscillations at this point in time. We stress that the peaks in $c$ at boundaries and oscillations at the bottom boundary are not numerical artifacts but are truly part of the solution of Eq. (3) without vapor fluxes at the boundaries. They must be seen as the deposition rates that are required so that the Hansen system assumption regarding the permanent and instantaneous restoration of the water vapor density to its saturated value is fulfilled when no contribution can be expected from water vapor mass fluxes at the domain boundaries to bridge the gap.

## 4.7 On the formation of density wave instabilities

In their work, Schürholt et al. (2022b) used a dedicated numerical setup to illustrate that both the Hansen and Calonne systems are prone to producing self-amplifying spatial oscillations in the $c$ and $\Phi_i$ fields when regions of very high density gradients are present, while linear stability analysis was utilized to confirm their results. These oscillations are true mathematical features related to the dependence of the effective heat and vapor diffusion coefficients on the ice volume fraction (Schürholt et al., 2022b). Here, we reproduce their experiment to investigate how this instability materializes when $\Phi_i$ is treated as an element-wise variable rather than a nodal variable. The numerical setup consists of a 2 cm thick snowpack with an initial ice volume fraction mimicking a Gaussian crust following Eq. (20) of Schürholt et al. (2022b) (solid black line in Fig. 9c). The initial conditions and BCs for $T$ and $\rho_v$ are the same as those used in Sect. 4.1. For these simulations, the time step is decreased to $\Delta t = 1$ min.

Figure 9 shows the $c$ and $\Phi_i$ fields obtained after 48 h of simulation on a 1000-node mesh. Let us compare first the solutions produced by FEniCS with an implicit time scheme (solid green line) and CC_3DOF with the same (improper) form of the mass matrix as for the FEniCS run and without lumping (dotted dark red line). Consistently with results already reported in Sect. 4.1, solutions obtained in both simulations are very close to each other. Concretely, RMSDs for $c$ and $\Phi_i$ are $1.7 \times 10^{-5}$ kg m$^{-3}$ s$^{-1}$ and $4.5 \times 10^{-4}$, respectively. In addition, both approaches produce smooth wave patterns – for both $c$ and $\Phi_i$ – encompassing several nodes at the bottom boundary. A closer look in this area (Fig. 9b, d) reveals slight differences. In particular, the amplitude of the oscillations in $\Phi_i$ are larger for the FEniCS solution than for the CC_3DOF one. This is not surprising given the nodal nature of $\Phi_i$ in the FEniCS approach: a deposition/sublimation peak at a given node translates directly into an ice volume
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison of deposition rate (a) and ice volume fraction (b) fields produced after 1, 2, 5, and 10 d of simulation on a 101-node mesh running (i) the corrected version of the code published by Simson et al. (2021), (ii) H_MF with boundary values of *c* forced to zero, and (iii) H_MF with no constraint on boundary values of *c* but with a no-flux BC on vapor. Note that, as it is an element-wise variable in our approach, *Φ*ᵢ is represented as steps for the two H_MF cases.</figcaption></figure>

fraction increase/decrease at the same node. In contrast, in our approach, the averaging of the nodal *c* over the elements when updating *Φ*ᵢ (Sect. 3.3) tends to limit the amplitude of these oscillations, but it does not erase them. We also note that, contrary to the results reported in the previous section, no strong sublimation/deposition peaks are observed at the very bottom/top nodes. In fact, the deposition rate is very slightly positive at both nodes. This is because the Dirichlet BCs $ρ_v = ρ_v^{eq}$ applied to the water vapor mass balance equation of the Calonne system imply water vapor fluxes at boundaries that contribute to maintaining the saturation at these nodes, limiting the necessity for deposition/sublimation to bridge the oversaturation/undersaturation.

Running the same CC_3DOF simulation with the proper form of the mass matrix (solid orange line in Fig. 9) considerably changes the profiles of *Φ*ᵢ and *c* obtained after 2 d of simulation. In particular, large oscillations in *Φ*ᵢ related

to large peaks in sublimation and deposition are observed on the sublimating (cold) side of the Gaussian crust. In fact, the general pattern of formation and propagation of the instability wave is not changed compared to the two previous simulations (as illustrated in Fig. 5 of Schürholt et al., 2022b), but oscillations set up at a much faster pace in this simulation. Again, this shows that treating mass matrices in their improper form can dramatically affect the obtained solutions at a given point in time.

We run an additional simulation that again consists of the CC_3DOF strategy with the proper form of the mass matrix and lumping activated but instead uses constant values for the effective parameters *k*<sup>eff</sup> and *D*<sup>eff</sup>. More specifically, we set *k*<sup>eff</sup> = 0.2 W m<sup>-1</sup> K<sup>-1</sup> and *D*<sup>eff</sup> = 1.066 × 10<sup>-5</sup> m<sup>2</sup> s<sup>-1</sup>. This corresponds to the Calonne parameterization of these effective parameters for a *Φ*ᵢ fixed to its initial averaged value, i.e., *Φ*ᵢ = 0.318. The obtained solution shows a nonuniform
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Comparison of deposition rate (a, b) and ice volume fraction (c, d) fields produced after 2 d of simulation on a 1000-node mesh by FEniCS, CC_3DOF without lumping and with the improper mass matrix, and CC_3DOF with lumping and with the proper mass matrix. The latter is run a second time with constant $D^{\text{eff}}$ and $k^{\text{eff}}$. The blue boxes in panels (a) and (c) highlight the zoomed-in areas depicted in panels (b) and (d). The solid black lines correspond to the initial conditions. Note that, despite the element-wise nature of $\Phi_i$ in our model, the field of $\Phi_i$ is represented as lines drawn from elemental values affected to the middle of elements in order to ease comparison with the FEniCS solution.</figcaption></figure>

deposition over the whole domain (solid red line in Fig. 9a and c). This is because the water vapor fluxes at boundaries implied by the previously mentioned Dirichlet BCs for water vapor bring enough water vapor mass from the exterior so that the whole layer is oversaturated for the implemented values of $k^{\text{eff}}$ and $D^{\text{eff}}$, which causes deposition. In contrast, when the same simulation is run with a no-flux condition on vapor at boundaries, we observe a very strong peak in sublimation at the bottom node and a less pronounced peak in deposition at the top node (not shown). The nonuniform deposition leads to an increase in $\Phi_i$ over the whole domain (more pronounced on the lower half than on the upper half). Logically, the quasi-advection of the Gaussian crust towards the warm boundary that is observed in all other simulations and thoroughly analyzed by Schürholt et al. (2022b) does not occur here. Indeed, the aforementioned quasi-advection of the Gaussian crust towards the warm boundary is due to continuous deposition on the lower side of the crust associated with continuous sublimation on the upper side, which does

not happen in this case. Another remarkable feature is the absence of the wave instability pattern in this case, in line with the prediction of Schürholt et al. (2022b).

All together, these observations confirm the assertion of Schürholt et al. (2022b): the wave patterns on c and $\Phi_i$ are intrinsic features of the mathematical models rather than due to, for example, the form of the implemented mass matrix or a violation of the DMP; they are triggered by strong density gradients and are due to the dependence of $k^{\text{eff}}$ and $D^{\text{eff}}$ on the ice volume fraction. There are slightly smoothened, but not removed, when the $\Phi_i$ is treated as an element-wise quantity rather than as a nodal quantity.

Note that we have also run the same experiment with the H_MF model using the Hansen parameterization of effective parameters (not shown). Again, our results are in line with those of Schürholt et al. (2022b): the wave instability appears much faster and with a much stronger amplitude than for the Calonne model with the Calonne parameterization of effective parameters.
## 4.8 On the effect of settlement on the density wave instabilities

We take advantage of our developments to run the same experiment as above but with settlement activated. As the considered domain is only 2 cm thick and the initial mean density is $291 \text{ kg m}^{-3}$, the stresses are too low and the initial viscosity is too high (viscosity increases exponentially with density) to induce significant settlement in affordable computational time if the parameterization (Eq. A7) of viscosity is taken as such. Therefore, we decide to keep this parameterization but to divide the obtained viscosity by a factor $10^4$. Given the linear nature of the implemented constitutive law, the general pattern of deformation is expected to be kept but with a strongly enhanced settlement rate. We limit the simulations to two runs, both performed with the CC_3DOF strategy with the proper form of the mass matrix and with lumping activated: one is run with settlement on, whereas the other is run with settlement off.

As illustrated in Fig. 10, including settlement has a stabilizing effect on the wave instability: oscillations in $c$ and $\Phi_i$ are still observed but with strongly reduced amplitudes. For the considered numerical setup, the ratio between the maximum amplitudes of the oscillations observed on $c$ ($\Phi_i$) when settlement is included and when it is omitted is ~40% (~30%) after 1 d and ~20% (~25%) after 2 d. These oscillations also take longer to set up when considering settlement: after 1 d of simulation, they are only starting to affect the $c$ field at the bottom nodes, and their feedback on the $\Phi_i$ field in this area is hardly perceptible. In contrast, when settlement is not included, the oscillatory pattern at the bottom boundary is already well established on both fields after 1 d of simulation. This stabilizing effect could be expected from the $\Phi_i$ dependence of viscosity: the denser the element, the less it deforms, and vice versa. This tends to homogenize the density of the snowpack and competes with phase change which, as shown previously, tends to further densify the denser layers and to further deplete the hollower layers. Another interesting feature is that, while only deposition is observed far enough from the Gaussian crust when settlement is off, these areas show sublimation when settlement is included.

## 5 Implications for future developments

The numerical issues that have been highlighted in the present work are relevant beyond the particular processes investigated here. First of all, it is a very general rule that the discretization of a continuous PDE must preserve its properties. First, this supposes the inclusion of any factor affecting the accumulation term in the corresponding mass matrix, but not its inverse in the stiffness matrix, as soon as such factor is nonuniform. Second, this contraindicates the use of the chain rule on the accumulation term during time discretization, as

it might break the conservative properties of the continuous PDE. In fact, in many situations, it is preferable to treat such a PDE in its mixed form, i.e., with the time derivative of the accumulation term directly applied to the conserved quantity. This is a generic result for conservation equations and has already been illustrated for several other processes (e.g., Celia et al., 1990; Casulli and Zanolli, 2010; Tubini et al., 2021). We have also shown that two-way-coupled PDEs need to be solved as single monolithic process to preserve their conservative properties. In addition, depending on the relative dynamics of the subprocesses at stake, solving the latter sequentially might require one to considerably decrease the time step in order to guarantee acceptable transient solutions. A monolithic treatment of coupled equations does not fundamentally increase the complexity and numerical cost of the solution. Indeed, solving two decoupled tridiagonal systems of equations or solving one coupled pentadiagonal system of equations both come with the same complexity (i.e., $O(n)$; El-Mikkawy and Atlan, 2014; Jia and Jiang, 2015). The possibility of using larger time steps with a monolithic treatment is, thus, beneficial in terms of computational cost. Again, this result holds as soon as two-way-coupled PDEs are at stake, in snow modeling or elsewhere.

More generally, coming up with a numerical implementation that guarantees energy and mass conservation is not straightforward. Therefore, a rigorous assessment of the evolution of energy and mass within the domain is a golden rule that should guide any numerical developments. This requires knowing the energy and mass fluxes in and out of the domain and, thus, explicitly computing the residuals of the algebraic system. In the eventuality that constraint on numerical cost requires numerical implementations violating the energy/mass conservation, the associated leakage should remain negligible, be well identified, and be tracked carefully. This is what we have done for the energy leakage associated with the sequential treatment of the feedback of deposition/sublimation on the ice mass conservation. Not only is this required to ensure that the model produces physically sound results but it is also required from the perspective of latter couplings to other components of Earth system models (ESMs). Indeed, the introduction of an artificial sink/source of energy/mass within ESMs obviously contributes to undesired model drift (e.g., Gupta et al., 2013; Hobbs et al., 2016). A good identification and a correct quantification of the spurious energy/mass leakage makes it easier to correct this drift.

More specific to the problem of water vapor diffusion is the treatment of water vapor at boundaries. We have shown (in Sect. 4.6) how sensitive the solution fields can be to this treatment. Previous works modeling macroscopic water vapor diffusion in snow have only considered two types of BCs for vapor: a no-flux condition (e.g., Touzeau et al., 2018) and a Dirichlet-type BC enforcing saturation of vapor at boundaries (e.g., Schürholt et al., 2022b). Jafari et al. (2020) used a mix of both, i.e., no flux at the bottom BC and saturation at the top BC. We also recall that forcing the deposition rate to
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Comparison of deposition rate (a) and ice volume fraction (b) fields produced after 1 and 2 d of simulation on a 1000-node mesh by CC_3DOF with lumping and the proper mass matrix in two configurations: without and with settlement. The solid black lines correspond to the initial conditions. We opt for the same representation of $\Phi_i$ as in Fig. 9.</figcaption></figure>

zero at boundary nodes when solving the Hansen system, as done by Simson et al. (2021), is a roundabout way of arbitrarily prescribing vapor fluxes at boundaries that perfectly compensate for the oversaturation/undersaturation at these nodes. In natural configurations, vapor fluxes are expected at both boundaries. It follows that the proper way of proceeding for real-world applications is to solve a surface water vapor mass balance in order to derive well-defined boundary water vapor fluxes, similarly to what is normally done for energy. In contrast, cold-laboratory experiments can impose a no-flux condition on vapor at boundaries by using impermeable plates. Such experiments reveal strong sublimation at the bottom BC and strong deposition at the top BC (e.g., supplement of Hagenmuller et al., 2019; Bouvet et al., 2023). These observations are consistent with the peaks in sublimation/deposition obtained at the bottom/top boundaries when vapor is forced to no flux, as presented in Sect. 4.6.

These peaks in sublimation/deposition and the related local increase in density gradients are a trigger for wave instability. This behavior is carefully analyzed in the study of Schürholt et al. (2022b), who also discuss the implications for the numerical solution. The existence of such wave instability was already proposed by Adams and Brown (1990) from theoretical arguments. It is due to the dependence of the effective thermal conductivity and vapor diffusion coefficient on the ice volume fraction. Physically, we understand the mechanism as follows: water vapor diffuses less readily across regions with a higher ice volume fraction; because the thermal conductivity of air is about 100 times less than that of ice, temperature gradients are stronger in regions with a low ice volume fraction, and vapor fluxes are enhanced;

these two phenomena combine to induce local accumulation of water vapor in places where the water vapor encounters an abrupt increase in the ice volume fraction in its flow towards decreasing temperatures; and this local accumulation leads to local oversaturation that is resorbed through deposition, further increasing the ice volume fraction at the transition. Similarly, when water vapor flows across an abrupt decrease in the ice volume fraction, enhanced water fluxes in the downstream low-density region combined with the difficulty in the upstream high-density region with respect to supplying water vapor lead to local undersaturation, which causes sublimation. Our developments have highlighted two stabilizing effects: one that is physical, i.e., the inclusion of settlement, and one that is numerical, i.e., the element-wise treatment of $\Phi_i$. Given these two effects, and acknowledging the fact that the numerical setup considered here is rather extreme (macroscopic temperature gradient of 1000 K m$^{-1}$ and 1000 elements for a 2 cm thick snowpack), we think that this instability wave is unlikely to become a major modeling difficulty for many of the natural settings. This point is more hazardous for situations such as seasonal-scale simulations of thin arctic snowpacks, and more simulations covering a large range of realistic configurations are needed.

As discussed by Simson et al. (2021), the surface ice mass balance related to new precipitation or to sublimation (for dry snow) also has to be addressed. We think that our choice regarding the element-wise treatment of conserved quantities will greatly facilitate this implementation. For simulations on a seasonal timescale, occasional re-meshing will be necessary to limit the computational cost. Again, the element-wise nature of conserved quantities will enable unequivocal redis-
tribution of the conserved quantities over merged or split elements. More generally, we draw the reader's attention to the widespread confusion between physical and numerical layers. Even when one claims to be adopting a continuous representation of the snowpack, the numerical solution of PDEs always requires discretization of the domain, which amounts to the definition of numerical layers. As soon as some integration is at stake (e.g., assessment of snowpack deformation from stress and viscosity), some assumption must necessarily be made about the evolution of fields within the numerical layers. Therefore, it is obviously preferable to have finer mesh in areas where gradients of relevant physical quantities are strong and a coarser mesh elsewhere. This is equivalent to saying that there must be some kind of correspondence between numerical and physical layers, if we define the latter as sections of the domain characterized by low gradients of physical properties.

The present work is a further step towards the development of the next generation of detailed snowpack models, associating sound and universal physics with a robust numerical treatment. So far, the implemented processes are limited to the coupled heat conduction, water vapor diffusion, and settlement in dry snow, but the numerical subtleties that have been highlighted are relevant for many other processes, in snow modeling and beyond. Our model is now ready to reproduce well-controlled laboratory experiments, provided that BCs are well constrained. In contrast, modeling natural settings will require more work, in particular regarding the energy, water vapor mass, and ice mass balances at boundaries. For longer-term developments, it will be necessary to implement other processes, notably those related to the liquid phase of water (melting/refreezing and percolation) as well as those regarding snow metamorphism.

# Appendix A: Parameterizations used in our model

## A1 Effective vapor diffusion coefficient

All simulations presented in the paper are run using the parameterization of the effective vapor diffusion coefficient proposed by Simson et al. (2021). It is based on the parameterization derived by Calonne et al. (2014), but it uses the Heaviside function $Θ$ to hinder vapor diffusion for ice volumes above two-thirds. It is expressed as follows:

$$D^{\text{eff}} = D_0 \left(1 - \frac{3}{2} \Phi_i\right) \Theta \left(\frac{2}{3} - \Phi_i\right), \quad (\text{A1})$$

with $D_0$ given in Table 1. Additional simulations mentioned in the paper were run using the parameterization of Hansen and Foslien (2015), which is expressed as follows:

$$D^{\text{eff}} = \Phi_i(1 - \Phi_i)D_0 + (1 - \Phi_i) \left( \frac{k_i D_0}{\Phi_i \left( k_a + L_m D_0 \frac{d\rho_{v,sat}}{dT} + (1 - \Phi_i)k_i \right)} \right), \quad (\text{A2})$$

with all constant values listed in Table 1.

## A2 Effective thermal conductivity

All simulations presented in the paper are run using the parameterization of Calonne et al. (2011) for the effective thermal conductivity, which is expressed as follows:

$$k^{\text{eff}} = k_0 + k_1 \rho_i \Phi_i + k_2 (\rho_i \Phi_i)^2, \quad (\text{A3})$$

where $k_0 = 0.024 \text{ W m}^{-1} \text{ K}^{-1}$, $k_1 = -1.23 \times 10^{-4} \text{ W m}^2 \text{ K}^{-1} \text{ kg}^{-1}$, and $k_2 = 2.5 \times 10^6 \text{ W m}^5 \text{ K}^{-1} \text{ kg}^{-2}$. Additional simulations mentioned in the paper were run using the parameterization of Hansen and Foslien (2015), which is expressed as follows:

$$k^{\text{eff}} = \Phi_i ((1 - \Phi_i)k_a + \Phi_i k_i) + (1 - \Phi_i) \left( \frac{k_i k_a}{\Phi_i \left( k_a + L_m D_0 \frac{d\rho_{v,sat}}{dT} + (1 - \Phi_i)k_i \right)} \right), \quad (\text{A4})$$

with all constant values listed in Table 1.

## A3 Effective heat capacity

Contrary to Simson et al. (2021) and Schürholt et al. (2022b), our parameterization of the effective heat capacity neglects the heat carried by dry air, which is then expressed as follows:

$$(\rho C_p)^{\text{eff}} = \rho_i C_i \Phi_i, \quad (\text{A5})$$

with all constant values given in Table 1. This assumption is justified by the small volumetric heat capacity of dry air compared with that of ice. This choice enables one to close the energy budget without having to track the dry air leaving the snowpack as it settles.

## A4 Saturation water vapor density

All simulations presented in this paper are run using the parameterization of Libbrecht (1999) for the saturation water vapor density:

$$\rho_v^{\text{eq}} = \frac{\exp(-T_r/T)}{fT} \left( a_0 + a_1(T - T_m) + a_2(T - T_m)^2 \right), \quad (\text{A6})$$

where $T_r = 6150 \text{ K}$, $f = 461.31 \text{ J K}^{-1} \text{ kg}^{-1}$, $a_0 = 3.6636 \times 10^{12} \text{ Pa}$, $a_1 = -1.3086 \times 10^8 \text{ Pa K}^{-1}$, $a_2 = -3.3793 \times 10^6 \text{ Pa K}^{-2}$, and $T_m = 273 \text{ K}$.

## A5 Effective viscosity

All simulations presented in the paper are run using the viscosity parameterization of Vionnet et al. (2012):

$$\eta(\Phi_i, T) = f \eta_0 \frac{\rho_i \Phi_i}{c_\eta} \exp(a_\eta(T_0 - T) + b_\eta \rho_i \Phi_i), \quad (\text{A7})$$
with a fusion temperature of $T_0 = 273$ K; the constant parameters $\eta_0 = 7.62237 \times 10^6$ kg s$^{-1}$, $a_\eta = 0.1$ K$^{-1}$, $b_\eta = 0.0.023$ m$^3$ kg$^{-1}$, and $c_\eta = 250$ kg m$^{-2}$; and an additional parameter $f$, which normally accounts for snow microstructure properties, that was set to $f = 1$ in the present study.

## Appendix B: Spatial and temporal discretization

The spatial discretization is based on the FEM. Indeed, the FEM enables a clear distinction between element-wise and nodal variables. The former are constant over elements and typically correspond to conserved quantities (e.g., ice volume fraction), whereas the latter are defined at mesh nodes, are continuous in space, and typically correspond to physical quantities driving the fluxes in between elements (e.g., temperature). As explained in further detail in Sect. 3.3, such a vision is consistent with the fact that the conservation of energy and mass are fulfilled on average per spatial interval rather than locally. In contrast to the work of Schürholt et al. (2022b), who used a Python-based FEM platform, we code every step of the method internally except the inversion of the linear system. This allows us to have complete control over the numerical implementation. Furthermore, this enables us to design a code structure that is well suited to snowpack modeling: each physical process corresponds to a module that can be easily activated or deactivated at the user's convenience. Below, we briefly recall the basics of the FEM. For more detail, we refer the reader to one of the many books that have been written on the subject (e.g., Pepper and Heinrich, 2005).

As the exact (analytical) solution $u$ of a PDE is usually out of reach, the idea of the FEM is to find an approximate (numerical) solution $\tilde{u}$ in the following form:

$$ \tilde{u}(\mathbf{x}, t) = \sum_{j=1}^{N_{\text{DOF}}} \varphi_j(\mathbf{x}) u_j(t) \quad \forall \mathbf{x} \in \Gamma, \qquad (\text{B1}) $$

where $\Gamma$ is the considered domain and the $u_j$ represents the discrete unknown scalar values of the problem to be solved. The functions $\varphi_j$ are linearly independent functions of space, referred to as the shape functions. Here, we follow the classical approach and adopt first-order Lagrangian polynomials as shape functions: the shape function $\varphi_j$ is a continuous piecewise linear function, whose value equals 1 at the $j$th node and 0 at all other nodes. In this case, the unknown scalar value $u_j$ simply corresponds to the value of $\tilde{u}$ at the $j$th node, and the shape functions $\varphi_j$ can be viewed as linear interpolators in between nodes. The restriction of shape functions to first-order polynomials is motivated by the results of Schürholt et al. (2022b), who reported that, in their experiments, increasing the polynomial order was equivalent to increasing the mesh resolution with the corresponding number of nodes.

The PDE is then rewritten in its weak form, which consists of multiplying the PDE by any arbitrary test function and

then integrating this product over $\Gamma$. Integration by parts is carried out on the divergence term to weaken the differentiation requirement of the solution field, which naturally makes boundary normal flux integrals arise. In order to obtain a closed set of discrete equations, the arbitrary test function is replaced by a finite set of test functions $\varphi_i$ that are taken equal to the shape functions (standard Galerkin procedure). The problem is then reduced to solving $N_{\text{DOF}}$ algebraic equations that can be cast into matrix forms, with matrix terms defined as integrals over space (see Appendix C). These integrals are evaluated using Gaussian quadratures. The Gaussian quadratures consist of replacing a continuous integral by a weighted sum of function values at specific points, the so-called Gaussian points, located within the elements. In particular, model parameters that depend on model variables are evaluated directly at Gaussian points at which nodal variables are interpolated through shape functions. When the considered model parameter is a nonlinear function of some model variable, this is different from evaluating the model parameter at nodes and then using shape functions to interpolate the obtained values at Gaussian point as done in studies such as Schürholt et al. (2022b). We use a default value of two Gaussian points per element, but this setting could easily be changed.

While the FEM provides a spatial discretization scheme, the obtained matrix equation also needs to be discretized in time. The $\theta$ method is the most commonly adopted. It consists of expressing all quantities as a weighted average between their values at the previous time step and their values at the current time step. In general, the matrix form of the resulting algebraic system is as follows:

$$ \left[ \mathbf{M} + \theta \Delta t \mathbf{K}^{n+1} \right] U^{n+1} = \underbrace{\left[ \mathbf{M} + (\theta - 1) \Delta t \mathbf{K}^n \right] U^n + \Delta t \left[ \theta \mathbf{F}^{n+1} + (1 - \theta) \mathbf{F}^n \right]}_{\text{right-hand side vector B}} + \Delta t \mathbf{F}^{\partial \Gamma}, \qquad (\text{B2}) $$

where the superscript $n+1$ ($n$) is applied to quantities evaluated at the current (previous) time step, $\Delta t$ is the time step size, and $\mathbf{U}$ the solution vector. The matrices $\mathbf{M}$ and $\mathbf{K}$ are referred to as the mass and stiffness matrices, respectively. The vector $\mathbf{F}$ is called the force vector, and the vector $\mathbf{F}^{\partial \Gamma}$ corresponds to boundary fluxes entering/leaving the domain during the time step. The particular cases of $\theta = 0$, $\theta = 0.5$, and $\theta = 1$ correspond to the first-order explicit Euler, the Crank-Nicolson, and the first-order implicit Euler methods, respectively. It can be shown that, whenever $0.5 \le \theta \le 1$, the timestepping algorithm is L2 stable: the error between the continuous time derivative and its discrete counterpart remains bounded without any requirement with respect to the time step size. However, if the considered PDE is nonlinear, the system (Eq. B2) also becomes nonlinear for any $\theta > 0$. Such a system must be linearized. Linearization algorithms, such
as the Picard or Newton techniques, are iterative methods that require the prescription of an initial guess, which must be sufficiently close to the solution to ensure convergence of the algorithm. The fulfillment of this condition hampers the prescription of an arbitrary large time step, even when the time discretization method is said to be unconditionally stable. On the other hand, while every case with $\theta \neq 0.5$ is first order in time (i.e., the discretization error decreases linearly with $\Delta t$, as $\Delta t$ converge to zero), the Crank–Nicolson method is second order in time (i.e., the discretization error decreases quadratically with $\Delta t$, as $\Delta t$ converges to zero). This higher accuracy of the Crank–Nicolson method and its stability are often invoked to justify its use (e.g., Bader and Weilenmann, 1992; Decharme et al., 2011; Schürholt et al., 2022b). However, the obtained solution can be affected by spurious (decaying) oscillations if the time step is too large or the element size too small. For this reason, we prefer to use the standard first-order implicit Euler method ($\theta = 1$). Although potentially less accurate than the Crank–Nicolson method, the latter is more stable and less prone to oscillations (Formaggia and Scotti, 2011).

Specific BCs can be implemented via the appropriate modification of the matrix system (Eq. B2). More specifically, Neumann BCs can be directly implemented through the boundary flux vector $F^{\partial\Gamma}$. A Dirichlet BC at a boundary node can be implemented by replacing its corresponding line in the left-hand-side matrix A of Eq. (B2) with the line of the identity matrix and setting the Dirichlet value in the right-hand-side vector B. Robin BCs must be seen as a weighted combination of Neumann and Dirichlet BCs in which the flux at the considered boundary is specified as a linear function of the value of the solution field at the corresponding boundary node. Implementing a Robin BC then consists of (i) adding the prescribed flux to the corresponding line of the boundary flux vector $F^{\partial\Gamma}$ and (ii) modifying the left-hand-side matrix A as for a Dirichlet BC, except that the diagonal term of the modified line, which is set to 1 when implementing a Dirichlet BC, must be set to the multiplying factor applying to the solution field in this case. The right-hand-side vector B is left unchanged.

Once the solution vector has been obtained, it is possible to diagnose the fluxes entering/leaving the system by calculating the matrix system residual, defined as $\Delta t F^{\partial\Gamma} = AU - B$. Note that this operation applies for all types of BCs, notably Dirichlet conditions for which the boundary fluxes are not explicitly prescribed but, nonetheless, exist to maintain the boundary solution at its prescribed value. Assessing the boundary fluxes is necessary to verify the closure of the energy budget and can be required for latter coupling to external models.

## Appendix C: Discrete forms of the systems of equations

This appendix presents the different FEM discretizations of the Calonne and Hansen systems for all considered strategies that are summarized in Fig. 1. In what follows, the functions $\varphi_i$ and $\varphi_j$ correspond to the test and shape functions, respectively.

### C1 CC_3DOF

The Calonne system is solved for the solution vector $u = (T, \rho_v, c)$. The discrete system evaluated at nonlinear iteration $k+1$ can be expressed in matrix form as follows:

$$ 
\left\{ 
\begin{aligned}
\mathbf{M}_T \dot{\mathbf{T}} + \mathbf{K}_{T,T} \mathbf{T} + \mathbf{K}_{T,C} \mathbf{C} &= \mathbf{F}_T^{\partial\Gamma} \\
\mathbf{M}_P \dot{\mathbf{P}} + \mathbf{K}_{P,P} \mathbf{P} + \mathbf{K}_{P,C} \mathbf{C} &= \mathbf{F}_P^{\partial\Gamma} \\
\tilde{\mathbf{M}}_C \mathbf{C} + \mathbf{K}_{C,T} \mathbf{T} + \mathbf{K}_{C,P} \mathbf{P} &= \mathbf{F}_C.
\end{aligned}
\right. 
\quad \text{(C1)} 
$$

Here,

$$ M_T^{i,j} = \int_{\Gamma} (\rho C_p)^{\text{eff}} \varphi_i \varphi_j \, dV, \quad \text{(C2)} $$

$$ M_P^{i,j} = \int_{\Gamma} (1 - \Phi_i) \varphi_i \varphi_j \, dV, \quad \text{(C3)} $$

$$ \tilde{M}_C^{i,j} = \int_{\Gamma} \varphi_i \varphi_j \, dV, \quad \text{(C4)} $$

$$ K_{T,T}^{i,j} = \int_{\Gamma} k^{\text{eff}} \nabla \varphi_i \cdot \nabla \varphi_j \, dV, \quad \text{(C5)} $$

$$ K_{P,P}^{i,j} = \int_{\Gamma} D^{\text{eff}} \nabla \varphi_i \cdot \nabla \varphi_j \, dV, \quad \text{(C6)} $$

$$ K_{T,C}^{i,j} = \int_{\Gamma} -L_m \varphi_i \varphi_j \, dV, \quad \text{(C7)} $$

$$ K_{P,C}^{i,j} = \int_{\Gamma} \varphi_i \varphi_j \, dV, \quad \text{(C8)} $$

$$ K_{C,T}^{i,j} = \int_{\Gamma} s \alpha v_{\text{kin}} (T^k) \left. \frac{d\rho_v^{\text{eq}}}{dT} \right|_{T=T^k} \varphi_i \varphi_j \, dV, \quad \text{(C9)} $$

$$ K_{C,P}^{i,j} = \int_{\Gamma} -s \alpha v_{\text{kin}} (T^k) \varphi_i \varphi_j \, dV, \quad \text{(C10)} $$

$$ F_C^i = \int_{\Gamma} s \alpha v_{\text{kin}} (T^k) \left( \left. \frac{d\rho_v^{\text{eq}}}{dT} \right|_{T=T^k} T^k - \rho_v^{\text{eq}} (T^k) \right) \varphi_i \, dV. \quad \text{(C11)} $$

In the above expressions, the saturation vapor density $\rho_v^{\text{eq}}$ has been linearized through Eq. (8), while $v_{\text{kin}}$ is fixed from the temperature field $T^k$ obtained at the previous nonlinear iteration. Vectors $F_T^{\partial\Gamma}$ and $F_P^{\partial\Gamma}$ correspond to normal boundary fluxes of sensible heat and vapor, respectively. Additional internal sources of energy besides the phase-change-related latent heat effect could easily be added as a force vector in the
energy budget equation. Finally, note that the matrix $\tilde{M}_C$ is not a true mass matrix in the sense that it does not apply to a time derivative; therefore, it must be treated carefully during the time-stepping assembly (Eq. B2).

## C2 CC_2DOF

The Calonne system is solved for the solution vector $u = (T, \rho_v)$. The discrete system evaluated at nonlinear iteration $k+1$ can be expressed in matrix form as follows:

$$ \left\{ \begin{aligned} \mathbf{M}_T \dot{\mathbf{T}} + \mathbf{K}_{T,T} \mathbf{T} + \mathbf{K}_{T,P} \mathbf{P} &= \mathbf{F}_T + \mathbf{F}_T^{\partial \Gamma} \\ \mathbf{M}_P \dot{\mathbf{P}} + \mathbf{K}_{P,P} \mathbf{P} + \mathbf{K}_{P,T} \mathbf{T} &= \mathbf{F}_P + \mathbf{F}_P^{\partial \Gamma}. \end{aligned} \right. \quad (C12) $$

Here,

$$ K_{T,T}^{i,j} = \int_{\Gamma} \left[ k^{\text{eff}} \nabla \varphi_i \cdot \nabla \varphi_j + L_m s \alpha v_{\text{kin}}(T^k) \left. \frac{d\rho_v^{\text{eq}}}{dT} \right|_{T=T^k} \varphi_i \varphi_j \right] dV, \quad (C13) $$

$$ K_{P,P}^{i,j} = \int_{\Gamma} \left[ D^{\text{eff}} \nabla \varphi_i \cdot \nabla \varphi_j + s \alpha v_{\text{kin}}(T^k) \varphi_i \varphi_j \right] dV, \quad (C14) $$

$$ K_{T,P}^{i,j} = \int_{\Gamma} -L_m s \alpha v_{\text{kin}}(T^k) \varphi_i \varphi_j dV, \quad (C15) $$

$$ K_{P,T}^{i,j} = \int_{\Gamma} -s \alpha v_{\text{kin}}(T^k) \left. \frac{d\rho_v^{\text{eq}}}{dT} \right|_{T=T^k} \varphi_i \varphi_j dV, \quad (C16) $$

$$ F_T^i = \int_{\Gamma} L_m s \alpha v_{\text{kin}}(T^k) \left( \left. \frac{d\rho_v^{\text{eq}}}{dT} \right|_{T=T^k} T^k - \rho_v^{\text{eq}}(T^k) \right) \varphi_i dV, \quad (C17) $$

$$ F_P^i = \int_{\Gamma} -s \alpha v_{\text{kin}}(T^k) \left( \left. \frac{d\rho_v^{\text{eq}}}{dT} \right|_{T=T^k} T^k - \rho_v^{\text{eq}}(T^k) \right) \varphi_i dV. \quad (C18) $$

In the above expressions, $\mathbf{M}_T$ and $\mathbf{M}_P$ are expressed as in the CC_3DOF case. Again, additional internal sources of energy could easily be added in the force vector $\mathbf{F}_T$ if needed, and the nonlinearity due to $\rho_v^{\text{eq}}$ and $v_{\text{kin}}$ is treated as for the CC_3DOF case. The field of $c$ is then diagnosed in a next step as the closure of the water vapor mass balance equation, where $\rho_v$ is set to the obtained solution field.

## C3 CD_PC

The equations of the Calonne system are solved sequentially. In a first step, the heat equation is solved for $T$ with a source term that is fixed from the $c$ field computed at the previous time step. The discrete counterpart of this equation can be expressed in matrix form as follows:

$$ \mathbf{M}_T \dot{\mathbf{T}} + \mathbf{K}_{T,T} \mathbf{T} = \mathbf{F}_T + \mathbf{F}_T^{\partial \Gamma}, \quad (C19) $$

where

$$ F_T^i = \int_{\Gamma} L_m c \varphi_i dV. \quad (C20) $$

In the above expressions, $\mathbf{M}_T$ and $\mathbf{K}_{T,T}$ are expressed as in the CC_3DOF case, and there is the possibility of adding additional internal sources of energy in $\mathbf{F}_T$. The obtained solution field is used to fix the distributions of $v_{\text{kin}}$ and $\rho_v^{\text{eq}}$. The water vapor mass balance equation can then be solved for $\rho_v$ under its diffusion-reaction form in a second step. The discrete counterpart of this equation can be expressed in matrix form as follows:

$$ \mathbf{M}_P \dot{\mathbf{P}} + \mathbf{K}_{P,P} \mathbf{P} = \mathbf{F}_P + \mathbf{F}_P^{\partial \Gamma}, \quad (C21) $$

where

$$ F_P^i = \int_{\Gamma} s \alpha v_{\text{kin}}(T^{k+1}) \rho_v^{\text{eq}}(T^{k+1}) \varphi_i dV. \quad (C22) $$

In the above expressions, $\mathbf{M}_P$ and $\mathbf{K}_{P,P}$ are expressed as in the CC_2DOF case. We stress that, because $v_{\text{kin}}$ and $\rho_v^{\text{eq}}$ are fixed from the field $T^{k+1}$ obtained in a previous (separated) step, no linearization is required in this case. The field of $c$ is computed in a third step as the closure of the water vapor mass balance equation, where $\rho_v$ is set to the obtained solution field.

## C4 CDFred

The CD Fred strategy is similar to the CD PC strategy except that the $c$ field computed at the previous time step is used to fix the source terms of both the heat and water vapor mass balance equations. It follows that the discrete counterparts of both equations are written in exactly the same manner as for the CD PC case, except for $\mathbf{K}_{P,P}$ and $\mathbf{F}_P$ which are expressed as follows:

$$ K_{P,P}^{i,j} = \int_{\Gamma} D^{\text{eff}} \nabla \varphi_i \cdot \nabla \varphi_j dV \quad (C23) $$

and

$$ F_P^i = \int_{\Gamma} -c \varphi_i dV. \quad (C24) $$

As for the CD PC case, no linearization is required in this approach. The field of $c$ is computed in a third step by solving Eq. (2) at all nodes, where $\rho_v$ and $\rho_v^{\text{eq}}(T)$ are set to the obtained solution fields.

## C5 H_MF

In the Hansen system, noting that $\mathbf{H}$ is the vector of enthalpy content, the total (i.e., sensible plus latent) energy budget can be expressed in matrix form:

$$ \left\{ \begin{aligned} \mathbf{M}_H \dot{\mathbf{H}} + \mathbf{K}_{T,T} \mathbf{T} &= \mathbf{F}_T^{\partial \Gamma} \\ \tilde{\mathbf{M}}_H \mathbf{H} + \mathbf{K}_{H,T} \mathbf{T} &= \mathbf{F}_H. \end{aligned} \right. \quad (C25) $$
Here,

$$M_H^{i,j} = \int_{\Gamma} \varphi_i \varphi_j dV, \qquad (C26)$$

$$\tilde{M}_H^{i,j} = \int_{\Gamma} \varphi_i \varphi_j dV, \qquad (C27)$$

$$K_{T,T}^{i,j} = \int_{\Gamma} \left( k^{\text{eff}} + D^{\text{eff}} L_m \left. \frac{d\rho_v^{\text{eq}}}{dT} \right|_{T=T^k} \right) \nabla \varphi_i \cdot \nabla \varphi_j dV, \qquad (C28)$$

$$K_{H,T}^{i,j} = \int_{\Gamma} \left( L_m(1 - \Phi_i) \left. \frac{d\rho_v^{\text{eq}}}{dT} \right|_{T=T^k} + (\rho C_p)^{\text{eff}} \right) \varphi_i \varphi_j dV, \qquad (C29)$$

$$F_H^i = \int_{\Gamma} \left( L_m(1 - \Phi_i)(\rho_v^{\text{eq}}(T^k) - \left. \frac{d\rho_v^{\text{eq}}}{dT} \right|_{T=T^k} T^k) - (\rho C_p)^{\text{eff}} T_0 \right) \varphi_i dV. \qquad (C30)$$

In the above expressions, the saturation vapor density $\rho_v^{\text{eq}}$ has been linearized through Eq. (8). The vector $F_T^{\partial\Gamma}$ corresponds to the normal heat flux (sensible and latent) at the boundary. This matrix representation is referred to as the mixed form, as the heat budget equation includes two unknowns, the total energy $H$ and the temperature $T$, which are related through a nonlinear constitutive equation. As for the CC_3DOF case, note that the matrix $\tilde{M}_H$ is not a true mass matrix and must be treated carefully during the time-stepping assembly (Eq. B2). The field of $c$ is then diagnosed in a next step as the closure of the water vapor mass balance equation with $\rho_v = \rho_v^{\text{eq}}$, where $\rho_v^{\text{eq}}$ is computed from the obtained field of $T$.

## C6 H_TF

By applying the chain rule in the Hansen system of equations, one can eliminate the total energy $H$ to express the system in terms of $T$ only. The equivalent matrix form is given by the following:

$$M_T \dot{T} + K_{T,T} T = F_H^{\partial\Gamma}, \qquad (C31)$$

where

$$M_T^{i,j} = \int_{\Gamma} \left( (\rho C_p)^{\text{eff}} + L_m(1 - \Phi_i) \left. \frac{d\rho_v^{\text{eq}}}{dT} \right|_{T=T^k} \right) \varphi_i \varphi_j dV \qquad (C32)$$

and the expression of $K_{T,T}$ is the same as for the H_MF case. Again, the field of $c$ is then diagnosed in a next step as the closure of the water vapor mass balance equation with $\rho_v = \rho_v^{\text{eq}}$, where $\rho_v^{\text{eq}}$ is computed from the obtained field of $T$.

## Appendix D: Modifications made to the published code of Simson et al. (2021)

After running the code of Simson and Kowalski (2021), three bugs (i.e., inconsistencies in the code that hamper mass conservation)

have been identified. Here, we detail these three bugs and explain the changes that were made to the lines of code that caused the problems, in order to fix these bugs. In addition to these bugs, we recall that the explicit temporal discretization of the continuity equation performed by Simson et al. (2021) also causes a violation of mass conservation (Sect. 3.3 and Appendix E). Although we consider this last point to be of a conceptual order rather than a bug (the problem appears both in the paper and in the code), we also present the modifications that are needed in the published code to switch from the original explicit time discretization to the implicit one.

## D1 Spatial inconsistency in the strain rate computation

We recall that the settlement scheme consists of two operations that must be done consistently to guarantee mass conservation (Sect. 3.3): (i) the update of the ice volume fraction field and (ii) the update of the mesh. The strain rate is required in both operations. In the approach of Simson et al. (2021), the strain rate is evaluated at nodes: the strain rate at node $k$ is calculated from the viscosity and stress at node $k$. As explained in Sect. 3.3, the stress at node $k$ is calculated from the overburden snow mass under the implicit assumption that the ice volume fraction $\Phi_{i,k}$ stored at node $k$ actually applies to the whole space interval located above, i.e., between nodes $k$ and $k+1$. On the other hand, in the published version of the code, the authors were updating the mesh using the strain rate evaluated at node $k$ to compute the deformation of the space interval below, i.e., between the nodes $k$ and $k-1$. This operation was done in the Model/velocity.py file via the following code lines (l. 214–220):

```
D_rate = 0

v[0] = D_rate[0] * dz[0] # local velocity at the lowest node v=0 [ms-1]
v[1:] = np.cumsum(
         D_rate[1:] * dz[:]
) # Integrate deformation rates [s-1] in space to derive velocity [ms-1]
```

The fact that the ice volume fraction stored at node $k$ relates to the snow mass contained in the element above, whereas the strain rate calculated at node $k$ is used for the deformation of the element below, leads to inconsistencies when ice volume fraction is updated. Therefore, we modified the code so that the strain rate used to deform the space interval between nodes $k$ and $k-1$ during the mesh update step is the one evaluated at node $k-1$ instead of the one evaluated at node $k$. Concretely, the code lines reported above are replaced by the following:

```
v[0] = 0 # local velocity at the lowest node v=0 [ms-1]
v[1:] = np.cumsum(
         D_rate[:-1] * dz :
) # Integrate deformation rates [s-1] in space to derive velocity [ms-1].
```
## D2 Inconsistency in the sequence of tasks

In the published version of the code, the strain rate was updated in between the ice volume fraction update and mesh update steps, leading to inconsistency between the two operations. This inconsistency was corrected by replacing the code lines

```
phi_new = update Phi(c, dt, nz, phi, v_dz)
```
```
 coord_new, dz, v_dz_new, v_new, sigma) = updateCoord(
    T, c, dt, nz, phi, coord, SetVel, v_opt, viscosity
)
```

of the Model/coupled_updatephi_coors.py file (l. 13–16) with

```
 (coord_new, dz, v_dz_new, v_new, sigma) = updateCoord(
    T, c, dt, nz, phi, coord, SetVel, v_opt, viscosity
)
```

```
phi_new = update phi(c, dt, nz, phi, v_dz_new)
```

## D3 Inconsistency in the relationship between $\Phi_i$ and the effective snow density

In the code of Simson and Kowalski (2021), which is consistent with Eq. (18) of their article (Simson et al., 2021), the stress is calculated from the overburden snow mass neglecting the contribution of the air mass, i.e., assuming that the effective snow density is $\rho_{\text{snow}} = \rho_i \Phi_i$. However, in the code, the initial field of the ice volume fraction $\Phi_i^0$ is retrieved from the user-prescribed initial field of snow density $\rho_{\text{snow}}^0$ as follows:

$$ \Phi_i^0 = \frac{\rho_{\text{snow}}^0 - \rho_a}{\rho_i - \rho_a}, \qquad (\text{D1}) $$

where $\rho_a = 1.335 \text{ kg m}^{-3}$ is the air density. To correct this inconsistency, line 19 of the Model/phi_from_rho_eff.py file, which was

is replaced with

## D4 Explicit vs. implicit time integration of the continuity equation

In the published version of the code, the continuity equation is discretized using an explicit time integration scheme. This is done in the Model/coupled_update phi coord.py

file via the following code line:

```
phi_new = phi + dt * (c / rho_i - v_dz * phi)
```

## Appendix E: Mass-conservative temporal discretization of the continuity equation

Let us consider the numerical layer $k+1/2$ undergoing settlement without any phase change between $t$ and $t+\Delta t$. If we denote $\Phi_{i,k+1/2}^n$ and $L_{k+1/2}^n$ ($\Phi_{i,k+1/2}^{n+1}$ and $L_{k+1/2}^{n+1}$), the ice volume fraction and length of the numerical layer before (after) settlement, and if we neglect the contribution of air in the snow mass, the mass conservation can be simply expressed as follows:

$$ \Phi_{i,k+1/2}^{n+1} L_{k+1/2}^{n+1} = \Phi_{i,k+1/2}^{n} L_{k+1/2}^{n}, \qquad (\text{E1}) $$

which can be rewritten as

$$ \Phi_{i,k+1/2}^{n+1} = \Phi_{i,k+1/2}^{n} \frac{L_{k+1/2}^{n}}{L_{k+1/2}^{n+1}}. \qquad (\text{E2}) $$

The total deformation of the layer $k+1/2$ over $\Delta t$ is as follows:

$$ L_{k+1/2}^{n+1} - L_{k+1/2}^{n} = \int_{z=z_k^n}^{z_k^n} \dot{\epsilon}_{zz}^{n+1}(z) dz \Delta t. \qquad (\text{E3}) $$

Thus, the spatially averaged strain can be expressed as follows:

$$ \frac{L_{k+1/2}^{n+1} - L_{k+1/2}^{n}}{L_{k+1/2}^{n}} = \frac{\int_{z=z_k^n}^{z_k^n} \dot{\epsilon}_{zz}^{n}(z) dz}{L_{k+1/2}^{n}} \Delta t = \dot{\epsilon}_{zz,k+1/2}^{n+1} \Delta t. \qquad (\text{E4}) $$

It follows that

$$ \frac{L_{k+1/2}^{n+1}}{L_{k+1/2}^{n}} = 1 + \dot{\epsilon}_{zz,k+1/2}^{n+1} \Delta t. \qquad (\text{E5}) $$

Combining Eqs. (E2) and (E5) gives

$$ \Phi_{i,k+1/2}^{n+1} = \frac{\Phi_{i,k+1/2}^{n}}{1 + \dot{\epsilon}_{zz,k+1/2}^{n+1} \Delta t}. \qquad (\text{E6}) $$

Equation (E6) guarantees mass conservation by construction (provided that the air contribution in the snow mass is neglected). Equation (19) is an extension of Eq. (E6) that also includes the contribution of the mean deposition within layer $k+1/2$. When phase change is present, the deposition/sublimation effectively changes the ice volume fraction affected to the layer.

*Code and data availability.* The code source files are provided at https://doi.org/10.5281/zenodo.7941767 to guarantee the permanent reproducibility of results (Brondex et al., 2023); however, we recommend that potential future users and developers access the code from its Git repository (https://github.com/jbrondex/ivori_model_homemadefem, last access: 25 September 2023) to benefit from the last versions of the code. The version used in this work is tagged as v0.1.0. For setting up the environment and running the simulations, please follow the instructions described in the README file present in the GitHub repository.
*Author contributions.* MD obtained funding and supervised the work. JB and FT reviewed the existing literature. JB implemented the model with constant support from KF and insights from MD, NC, PH, and HL. JB ran the numerical simulations. JB interpreted the results with help from all co-authors. JB wrote the manuscript with contributions from all co-authors.

*Competing interests.* The contact author has declared that none of the authors has any competing interests.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

*Acknowledgements.* We thank Clément Cancès for fruitful discussions on numerical methods. We thank Louis Le Toumelin for his help with managing the Git repository. CNRM/CEN is part of Labex OSUG (ANR10 LABX56). Julien Brondex, Kevin Fourteau, Marie Dumont, and Francois Tuzet have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation program (IVORI; grant no. 949516).

*Financial support.* This research has been supported by the European Research Council, H2020 European Research Council (grant no. 949516).

*Review statement.* This paper was edited by Ludovic Räss and reviewed by Mohammad Afzal Shadab and one anonymous referee.

## References

Adams, E. E. and Brown, R. L.: A mixture theory for evaluating heat and mass transport processes in nonhomogeneous snow, Continuum Mech. Therm., 2, 31–63, https://doi.org/10.1007/BF01170954, 1990.

Albert, M. R. and McGilvary, W. R.: Thermal effects due to air flow and vapor transport in dry snow, J. Glaciol., 38, 273–281, https://doi.org/10.1017/S0022143000003683, 1992.

Bader, H.-P. and Weilenmann, P.: Modeling temperature distribution, energy and mass flow in a (phase-changing) snowpack. I. Model and case studies, Cold Reg. Sci. Technol., 20, 157–181, https://doi.org/10.1016/0165-232X(92)90015-M, 1992.

Barrere, M., Domine, F., Decharme, B., Morin, S., Vionnet, V., and Lafaysse, M.: Evaluating the performance of coupled snow-soil models in SURFEXv8 to simulate the permafrost thermal regime at a high Arctic site, Geosci. Model Dev., 10, 3461–3479, https://doi.org/10.5194/gmd-10-3461-2017, 2017.

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

Bourbatache, M. K., Le, T. D., Millet, O., and Moyne, C.: Limits of Classical Homogenization Procedure for Coupled Diffusion-Heterogeneous Reaction Processes in Porous Media, Transport Porous Med., 140, 437–457, https://doi.org/10.1007/s11242-021-01683-2, 2021.

Bouvet, L., Calonne, N., Flin, F., and Geindreau, C.: Heterogeneous grain growth and vertical mass transfer within a snow layer under a temperature gradient, The Cryosphere, 17, 3553–3573, https://doi.org/10.5194/tc-17-3553-2023, 2023.

Brondex, J., Fourteau, K., Dumont, M., Hagenmuller, P., Calonne, N., Tuzet, F., and Löwe, H.: Supplementary to “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)”, Zenodo [code and data set], https://doi.org/10.5281/zenodo.7941767, 2023.

Brun, E., Martin, E., Simon, V., Gendre, C., and Coléou, C.: An energy and mass model of snow cover suitable for operational avalanche forecasting, J. Glaciol., 35, 333–342, https://doi.org/10.1017/S0022143000009254, 1989.

Brun, E., David, P., Sudul, M., and Brunot, G.: A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting, J. Glaciol., 38, 13–22, https://doi.org/10.1017/S0022143000009552, 1992.

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, Geophys. Res. Lett., 38, L23501, https://doi.org/10.1029/2011gl049234, 2011.

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

Casulli, V. and Zanolli, P.: A nested Newton-type algorithm for finite volume methods solving Richards’ equation in mixed form, SIAM J. Sci. Comput., 32, 2255–2273, https://doi.org/10.1137/100786320, 2010.

Celia, M. A., Bouloutas, E. T., and Zarba, R. L.: A General Mass-Conservative Numerical Solution for the Unsaturated Flow Equation, Water Resour. Res., 26, 1483–1496, https://doi.org/10.1029/WR026i007p01483, 1990.

Ciarlet, P.: Maximum principle and uniform convergence for the finite element method, Comput. Method. Appl. M., 2, 17–31, https://doi.org/10.1016/0045-7825(73)90019-4, 1973.

Decharme, B., Boone, A., Delire, C., and Noilhan, J.: Local evaluation of the Interaction between Soil Biosphere Atmosphere soil multilayer diffusion scheme using four pedotransfer functions, J. Geophys. Res.-Atmos., 116, D20126, https://doi.org/10.1029/2011JD016002, 2011.

Domine, F., Barrere, M., and Sarrazin, D.: Seasonal evolution of the effective thermal conductivity of the snow and the soil in high Arctic herb tundra at Bylot Island, Canada, The Cryosphere, 10, 2573–2588, https://doi.org/10.5194/tc-10-2573-2016, 2016.

Domine, F., Picard, G., Morin, S., Barrere, M., Madore, J.-B., and Langlois, A.: Major Issues in Simulating Some Arctic Snowpack Properties Using Current Detailed Snow Physics Models: Consequences for the Thermal Regime and Water
Budget of Permafrost, J. Adv. Model. Earth Sy., 11, 34–44,
https://doi.org/10.1029/2018MS001445, 2019.

El-Mikkawy, M. and Atlan, F.: Algorithms for solving linear systems of equations of tridiagonal type via transformations, Applied Mathematics, 5, 413–422, https://doi.org/10.4236/am.2014.53042, 2014.

Essery, R., Rutter, N., Pomeroy, J., Baxter, R., Stähli, M., Gustafsson, D., Barr, A., Bartlett, P., and Elder, K.: SNOWMIP2: An Evaluation of Forest Snow Process Simulations, B. Am. Meteorol. Soc., 90, 1120–1135, https://doi.org/10.1175/2009BAMS2629.1, 2009.

Etchevers, P., Martin, E., Brown, R., Fierz, C., Lejeune, Y., Bazile, E., Boone, A., Dai, Y.-J., Essery, R., Fernandez, A., Gusev, Y., Jordan, R., Koren, V., Kowalczyk, E., Nasonova, N. O., Pyles, R. D., Schlosser, A., Shmakin, A. B., Smirnova, T. G., Strasser, U., Verseghy, D., Yamazaki, T., and Yang, Z.-L.: Validation of the energy budget of an alpine snowpack simulated by several snow models (Snow MIP project), Ann. Glaciol., 38, 150–158, https://doi.org/10.3189/172756404781814825, 2004.

Formaggia, L. and Scotti, A.: Positivity and conservation properties of some integration schemes for mass action kinetics, SIAM J. Numer. Anal., 49, 1267–1288, https://doi.org/10.1137/100789592, 2011.

Fourteau, K., Domine, F., and Hagenmuller, P.: Macroscopic water vapor diffusion is not enhanced in snow, The Cryosphere, 15, 389–406, https://doi.org/10.5194/tc-15-389-2021, 2021.

Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a new-generation ice sheet model, Geosci. Model Dev., 6, 1299–1318, https://doi.org/10.5194/gmd-6-1299-2013, 2013.

Gupta, A. S., Jourdain, N. C., Brown, J. N., and Monselesan, D.: Climate Drift in the CMIP5 Models, J. Climate, 26, 8597–8615, https://doi.org/10.1175/JCLI-D-12-00521.1, 2013.

Hagenmuller, P., Flin, F., Dumont, M., Tuzet, F., Peinke, I., Lalalus, P., Dufour, A., Roulle, J., Pézard, L., Voisin, D., Ando, E., Rolland du Roscoat, S., and Charrier, P.: Motion of dust particles in dry snow under temperature gradient metamorphism, The Cryosphere, 13, 2345–2359, https://doi.org/10.5194/tc-13-2345-2019, 2019.

Hansen, A. C. and Foslien, W. E.: A macroscale mixture theory analysis of deposition and sublimation rates during heat and mass transfer in dry snow, The Cryosphere, 9, 1857–1878, https://doi.org/10.5194/tc-9-1857-2015, 2015.

Hobbs, W., Palmer, M. D., and Monselesan, D.: An Energy Conservation Analysis of Ocean Drift in the CMIP5 Global Coupled Models, J. Climate, 29, 1639–1653, https://doi.org/10.1175/JCLI-D-15-0477.1, 2016.

Irving, D., Hobbs, W., Church, J., and Zika, J.: A Mass and Energy Conservation Analysis of Drift in the CMIP6 Ensemble, J. Climate, 34, 3157–3170, https://doi.org/10.1175/JCLI-D-20-0281.1, 2021.

Jaafar, H. and Picot, J. J. C.: Thermal Conductivity of Snow by a Transient State Probe Method, Water Resour. Res., 6, 333–335, https://doi.org/10.1029/WR006i001p00333, 1970.

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 trans-

port on snow profiles in deep and shallow snow covers and on sea ice, Frontiers in Earth Science, 8, 249, https://doi.org/10.3389/feart.2020.00249, 2020.

Jia, J. and Jiang, Y.: Two symbolic algorithms for solving general periodic pentadiagonal linear systems, Comput. Math. Appl., 69, 1020–1029, https://doi.org/10.1016/j.camwa.2015.03.009, 2015.

John, V. and Schmeyer, E.: Finite element methods for time-dependent convection-diffusion-reaction equations with small diffusion, Comput. Method. Appl. M., 198, 475–494, https://doi.org/10.1016/j.cma.2008.08.016, 2008.

Jordan, R.: A One-Dimensional Temperature Model for Snow Cover, Technical Documentation for SNTHERM, Vol. 89, Special Report 91-16, US Army Corps of Engineers, 1991.

Krinner, G., Derksen, C., Essery, R., Flanner, M., Hagemann, S., Clark, M., Hall, A., Rott, H., Brutel-Vuilmet, C., Kim, H., Ménard, C. B., Mudryk, L., Thackeray, C., Wang, L., Arduini, G., Balsamo, G., Bartlett, P., Boike, J., Boone, A., Chéruy, F., Colin, J., Cuntz, M., Dai, Y., Decharme, B., Derry, J., Ducharne, A., Dutra, E., Fang, X., Fierz, C., Ghattas, J., Gusev, Y., Haverd, V., Kontu, A., Lafaysse, M., Law, R., Lawrence, D., Li, W., Marke, T., Marks, D., Ménégoz, M., Nasonova, O., Nitta, T., Niwano, M., Pomeroy, J., Raleigh, M. S., Schaedler, G., Semenov, V., Smirnova, T. G., Stacke, T., Strasser, U., Svenson, S., Turkov, D., Wang, T., Wever, N., Yuan, H., Zhou, W., and Zhu, D.: ESM-SnowMIP: assessing snow models and quantifying snow-related climate feedbacks, Geosci. Model Dev., 11, 5027–5049, https://doi.org/10.5194/gmd-11-5027-2018, 2018.

Libbrecht, K. G.: Physical properties of ice, http://www.cco.caltech.edu/~atomic/snowcrystals/ice/ice.htm, (last access: 2023-03-29), 1999.

Magnusson, J., Wever, N., Essery, R., Helbig, N., Winstral, A., and Jonas, T.: Evaluating snow models with varying process representations for hydrological applications, Water Resour. Res., 51, 2707–2723, https://doi.org/10.1002/2014WR016498, 2015.

Menard, C. B., Essery, R., Krinner, G., Arduini, G., Bartlett, P., Boone, A., Brutel-Vuilmet, C., Burke, E., Cuntz, M., Dai, Y., Decharme, B., Dutra, E., Fang, X., Fierz, C., Gusev, Y., Hagemann, S., Haverd, V., Kim, H., Lafaysse, M., Marke, T., Nasonova, O., Nitta, T., Niwano, M., Pomeroy, J., Schädler, G., Semenov, V. A., Smirnova, T., Strasser, U., Swenson, S., Turkov, D., Wever, N., and Yuan, H.: Scientific and Human Errors in a Snow Model Intercomparison, B. Am. Meteorol. Soc., 102, E61–E79, https://doi.org/10.1175/BAMS-D-19-0329.1, 2021.

Milly, P. C. D.: A mass-conservative procedure for time-stepping in models of unsaturated flow, Adv. Water Resour., 8, 32–36, https://doi.org/10.1016/0309-1708(85)90078-8, 1985.

Morin, S., Horton, S., Techel, F., Bavay, M., Coléou, C., Fierz, C., Gobiet, A., Hagenmuller, P., Lafaysse, M., Ližar, M., Mitterer, C., Monti, F., Müller, K., Olefs, M., Snook, J. S., van Herwijnen, A., and Vionnet, V.: Application of physical snowpack models in support of operational avalanche hazard forecasting: A status report on current implementations and prospects for the future, Cold Reg. Sci. Technol., 170, 102910, https://doi.org/10.1016/j.coldregions.2019.102910, 2020.

Patankar, S. V.: Numerical heat transfer and fluid flow, 1st edn, CRC Press, https://doi.org/10.1201/9781482234213, 1980.

Pepper, D. W. and Heinrich, J. C.: The finite element method: basic concepts and applications, Taylor & Francis, https://doi.org/10.1201/9781315395104, 2005.
Pfeffer, W. T. and Mrugala, R.: Temperature gradient and initial snow density as controlling factors in the formation and structure of hard depth hoar, J. Glaciol., 48, 485–494, 2002.

Pinzer, B. R., Schneebeli, M., and Kaempfer, T. U.: Vapor flux and recrystallization during dry snow metamorphism under a steady temperature gradient as observed by time-lapse micro-tomography, The Cryosphere, 6, 1141–1155, https://doi.org/10.5194/tc-6-1141-2012, 2012.

Protter, M. H. and Weinberger, H. F.: Maximum principles in differential equations, Springer Science & Business Media, ISBN 9781461252825, 2012.

Riche, F. and Schneebeli, M.: Thermal conductivity of snow measured by three independent methods and anisotropy considerations, The Cryosphere, 7, 217–227, https://doi.org/10.5194/tc-7-217-2013, 2013.

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.

Schürholt, K., Kowalski, J., and Löwe, H.: A numerical solver for heat and mass transport in snow based on FEniCS, EnviDat [code and data set], https://doi.org/10.16904/enividat.298, 2022a.

Schürholt, K., Kowalski, J., and Löwe, H.: Elements of future snowpack modeling – Part 1: A physical instability arising from the nonlinear coupling of transport and phase changes, The Cryosphere, 16, 903–923, https://doi.org/10.5194/tc-16-903-2022, 2022b.

Simson, A. and Kowalski, J.: Geo-fluid-dynamics/Eulerian_Lagrangian_snowSolution: final paper submission, TC, Zenodo [code], https://doi.org/10.5281/zenodo.5588308, 2021.

Simson, A., Löwe, H., and Kowalski, J.: Elements of future snow-pack modeling – Part 2: A modular and extendable Eulerian-Lagrangian numerical scheme for coupled transport, phase changes and settling processes, The Cryosphere, 15, 5423–5445, https://doi.org/10.5194/tc-15-5423-2021, 2021.

Slater, A. G., Schlosser, C. A., Desborough, C. E., Pitman, A. J., Henderson-Sellers, A., Robock, A., Vinnikov, K. Y., Mitchell, K., Boone, A., Braden, H., Chen, F., Cox, P. M., de Rosnay, P., Dickinson, R. E., Dai, Y. J., Duan, Q., Entin, J., Etchevers, P., Gedney, N., Gusev, Y. M., Habets, F., Kim, J., Koren, V., Kowalczyk, E. A., Nasonova, O. N., Noilhan, J., Schaake, S., Shmakin, A. B., Smirnova, T. G., Verseghy, D., Wetzel, P., Xue, Y., Yang, Z. L., and Zeng, Q.: The Representation of Snow in Land Surface Schemes: Results from PILPS 2(d), J. Hydrometeorol., 2, 7–25, https://doi.org/10.1175/1525-7541(2001)002<0007:TROSIL>2.0.CO;2, 2001.

Sturm, M. and Benson, C. S.: Vapor transport, grain growth and depth-hoar development in the subarctic snow, J. Glaciol., 43, 42–59, https://doi.org/10.1017/S0022143000002793, 1997.

Sturm, M. and Johnson, J. B.: Thermal conductivity measurements of depth hoar, J. Geophys. Res.-Sol. Ea., 97, 2129–2139, https://doi.org/10.1029/91JB02685, 1992.

Thomée, V.: On positivity preservation in some finite element methods for the heat equation, in: Numerical Methods and Applications: 8th International Conference, NMA 2014, Borovets, Bulgaria, 20–24 August 2014, Revised Selected Papers 8, Springer, 13–24, https://doi.org/10.1007/978-3-319-15585-2_2, 2015.

Touzeau, A., Landais, A., Morin, S., Arnaud, L., and Picard, G.: Numerical experiments on vapor diffusion in polar snow and firn and its impact on isotopes using the multi-layer energy balance model Crocus in SURFEX v8.0, Geosci. Model Dev., 11, 2393–2418, https://doi.org/10.5194/gmd-11-2393-2018, 2018.

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.

van Pelt, W. J. J., Oerlemans, J., Reijmer, C. H., Pohjola, V. A., Pettersson, R., and van Angelen, J. H.: Simulating melt, runoff and refreezing on Nordenskiöldbreen, Svalbard, using a coupled snow and energy balance model, The Cryosphere, 6, 641–659, https://doi.org/10.5194/tc-6-641-2012, 2012.

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

Wiese, M. and Schneebeli, M.: Snowbreeder 5: a Micro-CT device for measuring the snow-microstructure evolution under the simultaneous influence of a temperature gradient and compaction, J. Glaciol., 63, 355–360, https://doi.org/10.1017/jog.2016.143, 2017.

Yosida, Z., Oura, H., Kuroiwa, D., Huzioka, T., K., K., Aoki, S.-I., and Kinosita, S.: Physical Studies on Deposited Snow. I.: Thermal Properties, Contributions from the Institute of Low Temperature Science, 7, 19–74, http://hdl.handle.net/2115/20216 (last access: 29 November 2023), 1955.