# Some insights from the second principle of thermodynamics for snowpack modeling

Kévin Fourteau, Kaoane Jondeau, Clément Cancès, Marie Dumont

► To cite this version:

Kévin Fourteau, Kaoane Jondeau, Clément Cancès, Marie Dumont. Some insights from the second principle of thermodynamics for snowpack modeling. 2025. hal-05384623

HAL Id: hal-05384623

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

Preprint submitted on 27 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Some insights from the second principle of thermodynamics for snowpack modeling

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

**Correspondence:** Kévin Fourteau (kevin.fourteau@meteo.fr)

**Abstract.** Entropy and the second principle of thermodynamics are regularly used as an analysis tool in applied mathematics for physics-based numerical models. In essence, this approach states that the second principle (i.e. the non-destruction of entropy) is closely related to stability. Consequently, numerical models complying with the second principle are expected to be more robust than models that do not. A notable advantage of this method is its straight-forward generalization to nonlinear physics and to systems of coupled equations. The goal of this work is to thus investigate the added-value of such an entropy-based analysis to the case of snowpack modelling. For that, we study the conditions under which the physics describing snowpacks respects the second principle and the numerical schemes that preserve this compliance after temporal and spatial discretization. Specifically, we consider three cases of increasing complexity: (i) a dry snowpack governed by heat conduction only (meant to be an example of the method for unfamiliar readers), (ii) a system composed of a canopy and a snowpack exchanging heat, and (iii) a dry snowpack with heat conduction, vapor diffusion, and ice-vapor phase changes. Overall, the analysis shows that to comply with the second principle, numerical snowpack models should follow three main rules. First, physical variables should be co-localized. In other words, the temperature at a given point (and other intensive variables) should depend on the energy (and other extensive variables) only at that point. This property is naturally met with the finite volume method, but requires adaptation for the finite element method. Second, advected quantities, such as the enthalpy advected with vapor diffusion, should be numerically upstreamed. Finally, thermodynamical fluxes (for instance heat conduction or phase change rate) should be consistent with the end-of-timestep value of the thermodynamical gradients/differences that drive them (for instance temperature gradients or chemical potential differences). This can be achieved by employing a Backward Euler temporal integration and resolving the physical processes in a tightly-coupled manner. While proper compliance to the second principle is a good practice to build robust numerical models, it can however be cumbersome to achieve in practice. Therefore, we suggest to rather see this kind of entropy-based analysis as a tool, helping to diagnostic potential instability issues, rather than a rule to strictly follow.

## 1 Introduction

Snowpacks are a major component of the Earth system, impacting the energy and mass fluxes on continental surfaces (Estilow et al., 2015; Royer et al., 2021; Vavrus, 2007), on sea ice (Sturm and Massom, 2017; Webster et al., 2018), or on glaciarized

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

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

areas (Lenaerts et al., 2019). They are notably at the core of important climatic feedbacks (Thackeray et al., 2019). Numerical snowpack models are essential tools to quantitatively study the interaction between snow and the rest of its environment, from regional (Marty et al., 2017) to the global spatial scales (Decharme et al., 2016) and from the seasonal (Mott et al., 2023) to the millennial (Hoang et al., 2025) time scales. As such, on top of their ability to fairly represent the physics at play, numerical snowpack models are expected to run with relatively large timestep, from 15 min (Vionnet et al., 2012) to the day (Born et al., 2019). In particular, numerical instabilities and/or oscillations, which usually occur with large timesteps and small grid refinements (Butcher, 2008), but not only (Ciarlet and Raviart, 1973; Rank et al., 1983), can be a major drawback of models. It can indeed lead to physical inconsistencies or even abrupt crashes, jeopardizing the integrity of the whole simulation chain in which the snowpack models are integrated. This problem is exacerbated in some snowpack models, which aim to explicitly track layers of snow and thus have little control on their spatial grid in the case of very thin layers (such as an ice layer following a freezing rain event) present in the stratigraphy (e.g. the layering scheme described in Vionnet et al., 2012). For instance, the Crocus numerical model (Vionnet et al., 2012) can sometimes suffer from such instabilities (M. Lafaysse; personal communication, June 2025). An ad-hoc fix for this situation, is usually to decrease the timestep and to restart the simulation. While this approach is applicable in the case of seasonal and punctual simulations, it can rapidly become cumbersome for simulations over large areas and long periods of time. Thus, it is necessary to develop methodologies and techniques to favor the numerical robustness of numerical snowpack models.

30

35

40

Linear stability analysis can provide a first -yet still very informative- picture on the requirements for the development of robust numerical models, able to handle large timesteps (Butcher, 2008). Notably, this can be used to define L-stable methods, which are known to perform well in the case of stiff problems with large timesteps (Fazio, 2001). However, linear stability analysis naturally presents limitations when applied to non-linear system of equations. This is usually the case of snowpack models, that must represent various, potentially non-linear, mechanisms interacting together (Bartelt and Lehning, 2002; Vionnet et al., 2012). Another approach to study the stability (and further properties) of dissipative systems governed by partial differential equations (PDEs) is the so-called "entropy methods" (Jüngel, 2016). Roughly speaking, these methods build on the notion that the second law of thermodynamics (i.e. the non-destruction of entropy) is closely related to the notion of physical stability of material systems (Chapter V of Dafermos, 2010). The compliance to the second principle by a system of PDEs can be used to show that it is Lyapunov or asymptotically stable (Lyapunov, 1992). This methodology can then be extended to the discrete setting to study numerical approximations of such non-linear systems of PDEs (Chainais-Hillairet, 2014; Cancès and Guichard, 2016; Cancès and Guichard, 2017; Cancès and Gaudeul, 2020). In the past, several studies have discussed the consequences of the second principle in snow with a focus on metamorphism (Bartelt and Buser, 2004), avalanche flow (Bartelt et al., 2005, 2006), or energy and mass transfers at the snowpack scale (Bartelt et al., 2004). However, the application of the second principle to design numerical snowpack models has not been explored hitherto and is thus an interesting prospect.

45

50

55

The goal of this work is thus to explore the potential added value of applying entropy methods to the case of numerical snowpack models. A particular focus will be put on stability, even though entropy methods have further applications, for
instance to demonstrate the existence and uniqueness of solutions (Jüngel, 2015) or to investigate the long-time behavior of dissipative systems (Carrillo et al., 2001; Bessemoulin-Chatard and Chainais-Hillairet, 2017; Chainais-Hillairet and Herda, 2020). A slightly different point of view is to see conditions under which numerical snowpack models are consistent with thermodynamics and the discrete version of the second principle. This is by itself a desired property, as the second principle imposes constraints on what is physically possible or not.

To this purpose, the article is organized as follows. First, we briefly explain the rationale behind entropy methods (Sect. 2) and present a concrete application based on the heat equation (Sect. 3). While these two sections do not present original research, they are meant as a familiarization with entropy methods. Then, we study two cases relevant to snowpack modeling, namely the use of a snowpack model in a coupled snow-canopy situation (Sect. 4) and the modeling of a snowpack with internal vapor transport (Sect. 5). Finally, we discuss the added value of respecting the second principle and provide some general considerations for the development of robust numerical snowpack models.

# 2 Entropy methods: thermodynamics as a criterion for numerical robustness

## 2.1 Classical Irreversible Thermodynamics and the second principle

The framework in which entropy methods are deployed is usually that of Classical Irreversible Thermodynamics (CIT, sometimes also referred to as local equilibrium thermodynamics; Lebon et al., 2008). In this vision, macroscopically out-of-equilibrium systems are assumed to still follow standard equilibrium thermodynamics at the local scale (i.e. at a scale large enough to be a thermodynamical system, but small enough to locally equilibrate much faster than the whole system under consideration). For instance, a system could be characterized (among others) by an internal energy density field $u$ (expressed in J m⁻³) that is related to its temperature field $T$ by the same function $u = u(T)$ that hold in the case of overall equilibrium. Note that this is the same framework under which snowpack models are naturally built. Thus, current versions of snowpack models are directly suited to be analyzed with entropy methods.

In the CIT framework, the entropy density $s$ (expressed in J K⁻¹ m⁻³) is governed by a balance equation at the snowpack scale that reads

$$ \partial_t s + \nabla \cdot F_s = \psi_s + \psi_{ext} \quad (1) $$

where $\partial_t$ and $\nabla \cdot$ respectively stand for the partial derivative with respect to time and the spatial divergence, $F_s$ is the entropy flux, $\psi_s$ stands for the entropy production rate (also referred to as the internal entropy source in what follows), and $\psi_{ext}$ denotes an external entropy source. The second principle of thermodynamics simply states the internal entropy source must be non-negative, i.e. $\psi_s \ge 0$, everywhere and every time (Lebon et al., 2008). Note that the external entropy source $\psi_{ext}$ has been introduced to account for the distributed absorption of external energy directly within the snowpack (i.e. shortwave absorption), which is an entropy flux from the exterior but cannot be treated as a boundary value for $F_s$.
## 2.2 Lyapunov functionals and stability

As stated in the introduction, one of the applications of entropy methods is to use the second principle to demonstrate the stability of time-evolving systems. This is done in the context of Lyapunov stability and using Lyapunov functionals (Lyapunov, 1992).

Let's consider a dynamical (i.e. time-evolving) system characterized by the set of prognostic variables {$X$} (in the context of a snowpack this would for instance include the internal energy field $u$). A Lyapunov functional is a real-valued function $V(\{X\})$ such that

$$ 
\left\{ 
\begin{array}{l}
V(\{X_{eq}\}) = 0, \\
V(\{X\}) > 0 \quad \text{for } \{X\} \neq \{X_{eq}\}, \\
\frac{dV}{dt} \le 0.
\end{array}
\right. 
\qquad (2) $$

where {$X_{eq}$} is a point in the phase space of {$X$} that corresponds to an equilibrium state. The existence of such a Lyapunov functional directly confers stability to the dynamical system. Indeed, $V$ can be seen as a pseudo-distance in the phase space to the equilibrium point {$X_{eq}$}. Since $V$ cannot increase over time, the dynamical system cannot diverge away from the equilibrium. Moreover, in the case where $\frac{dV}{dt} < 0$, the system is drawn towards the equilibrium point over time.

Thermodynamical potentials, like e.g. negentropy or Helmholtz free energy (Kardar, 2007) naturally play the role of such Lyapunov functionals (Lebon et al., 2008). For instance, let us consider an isolated system, that is with no energy nor matter exchange with the exterior. In this case, the total negentropy $-S$ of the system is the building block of the Lyapunov functional. Indeed, it's time derivative is given by

$$ \frac{d(-S)}{dt} = \int_{\Omega} -\partial_t s \, dV \qquad (3) $$

with $\Omega$ the spatial domain of the system. Using Eq. (1) together with the divergence theorem and the no-flux boundary condition, directly yields that

$$ \frac{d(-S)}{dt} = \int_{\Omega} -\psi_s \, dV \le 0, \qquad (4) $$

the last inequality following from the second principle. Thus, $S(\{X_{eq}\}) - S$ is Lyapunov functional whose decay follows from the second principle and the non-destruction of entropy. Similarly, in the case of a system with constant temperature $T$ enforced on the boundary for instance, one can build a Lyapunov functional based on the Helmholtz free energy $A = U - TS$ (with $U$ and $S$ the total internal energy and entropy of the system), as $A - A(\{X_{eq}\})$. Again, the decreasing nature of the Lyapunov functional is related to the non-destruction of entropy.
## 2.3 Approach followed in the article

As briefly hinted above, the relevant thermodynamic potential -negentropy, Helmholtz free energy, etc- depends on the boundary conditions at play. However, in the case of more complex boundary conditions (for instance an inhomogeneous temperature field), the expression of a general thermodynamics potential remains an open question (Glansdorff and Prigogine, 1970; Keizer and Fox, 1974; Di Vita, 2013; Bodineau et al., 2014; Maes and Netočný, 2015). We stress that this does not necessarily mean that systems with inhomogeneous boundary conditions are not Lyapunov stable, as the ignorance of such thermodynamical potential in inhomogeneous boundary conditions does not mean that no Lyapunov functional exists. Therefore, for the rest of the article we do not attempt to build the proper thermodynamical potentials, consistent with a particular set of external forcings. Rather, we assume that the respect of the second principle -$\psi_s \ge 0$, everywhere, every time- is a sufficient condition to avoid thermodynamical instabilities (or at least that it is a condition that strongly favors stability; this last point will notably be discussed in Sect. 6.2).

Moreover, we believe that the use of $\psi_s > 0$ as our guiding criterion, also offers a slightly richer picture than the use of a Lyapunov functional. Indeed, $\psi_s \ge 0$ is a local criterion, and is thus sensitive to local violations of the second principle, while a Lyapunov functional is a global criterion which can in principle let local violations of the second principle occur (entropy can be locally destroyed as long as there are sufficient production elsewhere to overall decrease the Lyapunov functional). Notably, strict adherence to the local form of the second principle tends to hinder the occurrence of spontaneous local divergence, and oscillations, and overshoots as illustrated in Fig. 1. This is a desirable property for numerical models, in which artificial divergence, overshoots, and oscillations are to be avoided.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Example of spontaneous temporal evolutions over a time interval $\Delta t$ of a simple temperature field that are inconsistent with the second principle of thermodynamics, and where entropy has been locally destroyed. Systems respecting the second principle should not present such features.</figcaption></figure>
135 **3 Revisiting the heat equation**

We start by analyzing the heat equation (i.e. internal energy conservation without matter fluxes), as it is at the core of most snowpack models (e.g. Bartelt and Lehning, 2002; Vionnet et al., 2012; Essery, 2015; Sauter et al., 2020). While the entropy-based analysis of the heat equation has already been carried out in the literature (e.g. Vidal et al., 1994; Lebon et al., 2008), we believe that its brief presentation in the article offers an easy but concrete example of how the entropy analysis can be carried out, on which the rest of the article is built.

140

**3.1 The continuous model**

Based on classical results from the field of CIT, we start by analyzing the conditions under which the continuous (i.e. not discretized) heat equation complies with the second principle of thermodynamics. As stated above, the goal here is to briefly present the major steps in an entropy-based analysis. Note that a similar demonstration is for instance available in Lebon et al. (2008).

145

**3.1.1 Gibbs' equation**

The first step of the entropy method is to define the link between the entropy and the rest of the thermodynamical variables describing the snow material. This is done through Gibbs' equation, also sometimes referred to as the fundamental equation of thermodynamics (Kardar, 2007), which relates the differential of entropy to the differential of other extensive quantities. If one neglects mass transfers, the Gibbs' equation in snow viewed as a homogeneous medium simply writes

150

$$ds = \frac{du}{T}, \qquad (5)$$

where $s$ is the volumetric density of entropy (J K⁻¹ per m³ of snow), $u$ the volumetric density of internal energy (J m⁻³), and $T$ the temperature (K). Two observations can be made from Gibbs' equation. First, in this framework, $s$ can be expressed as a function of the internal energy (i.e. $s = s(u)$). Stability arguments require this function to be concave (Kardar, 2007). Second, the temperature $T$ is directly given by this entropy function as $\frac{1}{T} = \frac{ds}{du}$, and thus we also have $T = T(u)$. As will be done in Sect. 3.1.2 below, specifying how $s$ (and therefore $T$) depends on $u$, amounts to defining the material properties of snow. From Eq. (5), it follows that

155

$$\partial_t s = \frac{\partial_t u}{T}. \qquad (6)$$

**3.1.2 Defining the material law**

To close the thermodynamical description of the snow material, one as to define the relation between $s$ and $u$, and by the same operation define the relation between $T$ and $u$. Here, we assume that

160

$$dT = \frac{du}{c_p} \qquad (7)$$
with $c_p$ the thermal capacity (expressed in J K$^{-1}$ m$^{-3}$). We note that $c_p$ does not need to be a constant. The condition that $s(u)$ is a concave function simply requires the positivity of $c_p$, which makes physical sense: for thermodynamics stability, the temperature must increase with internal energy.

### 3.1.3 Prognostic equation: energy conservation

In the framework developed so far, snow is described by a single prognostic variable, namely the internal energy $u$. The description of a macroscopically out-of-equilibrium snowpack thus involves a prognostic equation for the internal energy, this equation simply being the balance of internal energy

$$ \partial_t u + \nabla \cdot (F_u) = Q \qquad (8) $$

where $F_u$ is the heat flux (expressed in W m$^{-2}$) and $Q$ is a volumetric energy source/sink provided by the exterior of the snowpack (e.g. penetrating shortwave radiations). It is this prognostic equation (complemented with material laws to compute $F_u$) that is concretely solved in snowpack models to simulate the evolution of a snowpack. Also, boundary conditions need to be prescribed for Eq. (8), but they do not need to be further discussed here as they do not impact the internal entropy source of the snowpack.

### 3.1.4 Expressing the internal entropy source

Multiplying Eq. (8) by $1/T$, using Eq. (6), and realizing that $\frac{1}{T}\nabla \cdot (F_u) = \nabla \cdot (\frac{F_u}{T}) - F_u \nabla \cdot (\frac{1}{T})$ provides the entropy balance equation

$$ \partial_t s + \nabla \cdot \left(\frac{F_u}{T}\right) = F_u \cdot \nabla \left(\frac{1}{T}\right) + \frac{Q}{T}, \qquad (9) $$

which allows us to identify the entropy flux $F_s = \nabla \cdot (\frac{F_u}{T})$, the external entropy source $\psi_{ext} = \frac{Q}{T}$, and the internal entropy production $\psi_s = F_u \cdot \nabla (\frac{1}{T})$. We note that this entropy balance equation directly follows from the heat equation, and therefore does not constitute a second prognostic governing the evolution of the snowpack. Rather, it is a diagnostic equation, that provide information about the heat equation, but does not modify its behavior.

The second principle of thermodynamics requires the internal entropy source $\psi_s$ to be positive under all potential circumstances. This simply suggests that $F_u \propto \nabla (\frac{1}{T})$, which is readily the case if the heat flux obeys Fourier's law, i.e. $F_u = -\lambda \nabla T$, with $\lambda \ge 0$ the thermal conductivity of snow. Note that we do not need here to assume that $\lambda$ is constant in space or that is independent of temperature. The second principle, which encodes the dissipative nature of snow, is fulfilled as long as heat flows down the temperature gradients, might it be in a linear or non-linear fashion.

## 3.2 Time and space discretization

While in the continuous domain, the heat equation (complemented by Fourier's law) naturally fulfills the second principle, this does not necessarily translate when it is spatially and temporally discretized as it is the case in numerical models. The
goal of this section is thus to determine the discretization conditions under which the entropy structure of the heat equation is preserved, or in other words under which numerical schemes does the discretized entropy source remain positive.

We consider the finite volume method (FVM, Eymard et al., 2000) for the spatial discretization, as it is the numerical discretization used in the most snowpack models for solving the heat equation (e.g. Vionnet et al., 2012; Essery, 2015; Sauter et al., 2020). For the sake of simplicity, we develop our purpose in the simple one-dimensional in space setting, which is natural and of practical interest in the numerical modeling of snowpacks. The extension to several space dimension is straightforward as long as one restricts to so-called two-point flux approximations (Droniou, 2014; Cancès, 2018).

In the one-dimensional FVM framework, the snowpack is divided into $N$ cells (sometimes referred to as layers in snowpack models, even though they should not be confused with the physical strata building the snowpack) and an energy balance integrated over each cell. We thus obtain an equation for each cell reading

$$ \Delta z_k \frac{u_k^{n+1} - u_k^n}{\Delta t} + F_{k+\frac{1}{2}}^{n \to n+1} - F_{k-\frac{1}{2}}^{n \to n+1} = \Delta z_k Q_k^{n \to n+1}, \quad (10) $$

where $\Delta t = t^{n+1} - t^n$ the duration of the $n^{th}$ time step, $\Delta z_k$ is the thickness of the $k^{th}$ cell, $u_k$ its internal energy (the superscripts $n$ and $n+1$ indicates that the variables is taken at the start or end of the timestep, respectively), $Q_k^{n \to n+1}$ the external volumetric energy source applied to cell $k$ during the timestep, and $F_{k-\frac{1}{2}}^{n \to n+1}$ and $F_{k+\frac{1}{2}}^{n \to n+1}$ are the heat fluxes between the cells $k-1$ and $k$ and cells $k$ and $k+1$ (or the heat fluxes from or towards the exterior of the snowpack if the cell $k \in \{1, N\}$ is at a boundary) during the timestep. Together with an expression for the fluxes $F_{k+\frac{1}{2}}^{n \to n+1}$ to be specified, this is the prognostic system of equations that is actually solved in the snowpack model and governs its evolution over time. We used a general time superscript for the fluxes in Eq. (10), letting free the choice of the time integration scheme so far.

As in the continuous case, an FVM entropy source can be derived from Eq. (10). For that, we first note that in the FVM framework the entropy is defined per cell and time step, and is simply expressed as $s_k^n = s(u_k^n)$, with $s(u)$ the same function as in the continuous domain. To obtain the discrete entropy source, we consider the variation of the total entropy of the system $S^n = \sum_{k=1}^N \Delta z_k s_k^n$. To this end, we multiply Eq. (10) by $1/T_k^{n+1}$ and sum over the all cells. Rearranging and using the fact that the heat flux between cells $k$ and $k+1$ appears both in the energy balance of cell $k$ and $k+1$, we obtain that

$$ \frac{1}{\Delta t} (S^{n+1} - S^n) = \Psi^{n+1} + \Psi_{\text{ext}}^{n+1} + \mathcal{F}^{n+1}, \quad (11) $$

with

$$ \Psi_{\text{ext}}^{n+1} = \sum_{\text{cells}} \Delta z_k \frac{Q_k^{n \to n+1}}{T_k^{n+1}}, $$

the entropy source due to external shortwave radiation,

$$ \mathcal{F}^{n+1} = \frac{F_{\text{bot}}^{n \to n+1}}{T_{\text{ext,bot}}} + \frac{F_{\text{top}}^{n \to n+1}}{T_{\text{ext,top}}}, $$
the total entropy flux at the boundary of the domain - with $F_{\text{bot}}^{n \to n+1}$ and $F_{\text{top}}^{n \to n+1}$ the bottom and top heat fluxes orientated towards the snowpack and $T_{\text{bot}}$ and $T_{\text{top}}$ the temperature of the external physical bodies responsible for the boundary heat fluxes -, and

$$ \Psi^{n+1} = \underbrace{\frac{1}{\Delta t} \sum_{\text{cells}} \Delta z_k \left( s_k^{n+1} - s_k^n - \frac{u_k^{n+1} - u_k^n}{T_k^{n+1}} \right)}_{\mathcal{A}^{n+1}} + \underbrace{\sum_{\text{internal edges}} F_{k+\frac{1}{2}}^{n \to n+1} \left( \frac{1}{T_{k+1}^{n+1}} - \frac{1}{T_k^{n+1}} \right)}_{\mathcal{B}^{n+1}} + \underbrace{F_{\text{bot}}^{n \to n+1} \left( \frac{1}{T_1^{n+1}} - \frac{1}{T_{\text{ext,bot}}} \right) + F_{\text{top}}^{n \to n+1} \left( \frac{1}{T_N^{n+1}} - \frac{1}{T_{\text{ext,top}}} \right)}_{c^{n+1}} \quad (12) $$

the discrete entropy production term. It is this production term that should be non-negative everywhere to respect the second principle. As seen in the equation above, it can be naturally split into three different contributions.

The $\mathcal{A}^{n+1}$ term corresponds to entropy created directly in each cell by the time integration scheme. It has a pure numerical nature and no continuous counterpart. It is positive due to the concavity of the entropy function $s(u)$, as for all $u$, there holds $s(u) \le s(u_k^{n+1}) + \frac{ds}{du}(u_k^{n+1})(u - u_k^{n+1})$.

Noting that $\frac{ds}{du}(u_k^{n+1}) = 1/T_{k+1}^{n+1}$ and taking $u = u_k^n$, it directly follows that $\mathcal{A}^{n+1} \ge 0$. We note that for the demonstration to hold, it is important for the variables $u$ and $s$ to be co-localized, i.e. that the entropy attached to a cell is a function of the energy of that cell only. While this property is natural in the FVM framework, it is not necessarily the case with the finite elements method (FEM), where the entropy at a given point depends on the energy at neighboring points. This potential entropy destruction through the term $\mathcal{A}^{n+1}$ is related to the potential violation of the discrete principle maximum of the FEM, and can be resolved through mass lumping techniques (Chatzipantelidis et al., 2015).

The $\mathcal{B}^{n+1}$ term is defined on every internal edge (i.e. the interfaces between layers) and directly corresponds to a discretized version of the entropy source $F_u \cdot \nabla (\frac{1}{T})$ in the continuous domain. It is composed of the product of the energy flux crossing each internal interface and of the difference in inverse temperature across the interface. It is positive as long as the numerical heat flux between two cells (i.e. $F_{k+\frac{1}{2}}^{n \to n+1}$) is orientated consistently with the temperature gradient at the end of the timestep. This consistency is naturally obtained with Fourier's law and a backward Euler temporal scheme, in which $F_{k+\frac{1}{2}}^{n \to n+1}$ is based on the temperatures $T_k^{n+1}$ and $T_{k+1}^{n+1}$. In particular, the Forward Euler or Crank-Nicholson temporal schemes cannot ensure the full consistency between the numerical heat fluxes and the temperature gradients at the end of the timestep, potentially resulting in local entropy destruction. We note that one might be tempted to prove the positivity of $\mathcal{B}^{n+1}$ by dividing Eq. (10) by $1/T_k^n$ (rather than $1/T_k^{n+1}$) in the derivation of Eq. (11). However, in this case, the $\mathcal{A}^{n+1}$ term would become nonpositive, leading to some entropy destruction. More advanced strategies to go to higher order in time while maintaining the compatibility with the second principle of thermodynamics exist (see for instance Cancès et al., 2024), but will not be further discussed in this paper, as they are seldom used in snow modelling.
Finally, the $C^{n+1}$ term corresponds to the entropy production due to the exchange of energy between the discrete snowpack and its surrounding. It is of the similar nature as $B^{n+1}$, and is positive as long as the energy flux at a boundary is physically consistent (i.e. energy flows from the warmer to the colder region) with the temperature difference between the snowpack and the exterior at the end of the timestep. As with $B^{n+1}$, this condition is naturally fulfilled in the case of a backward Euler temporal scheme, provided that the expression of the boundary heat flux is consistent with the physical notion that energy must flow from the warm to the cold regions.

In summary, there are two main criteria for the second principle to unconditionally hold when numerically solving the heat equation: (i) variables must be co-localized, which is directly the case with a standard FVM spatial discretization, and (ii) heat fluxes and temperature gradients/differences at the end of the timestep must be consistent, which is directly the case with a backward Euler temporal discretization. Note that to derive these results, we did not have to assume the system of equations to be linear.

# 4 Coupling of snowpack models with other components

In the analysis carried out so far, the system we studied consisted of a single snowpack, forced by independent external conditions (under the form of boundary conditions and shortwave absorption). However, snowpack models are regularly coupled with other numerical models, in order to represent more complex systems, of which the snowpack only represents a subcomponent. This is for instance the case when modelling the energy budget of the soil-snowpack-canopy continuum (Boone et al., 2017; Essery et al., 2025), which is essentially achieved by coupling a soil, a snowpack, and a canopy model.

The goal of the present section is to show that the entropy-based analysis can easily be extended to such coupled systems, and can similarly provide conditions to favor their numerical stability. For that, we consider the concrete example of a snowpack and a big leaf canopy, exchanging energy through longwave radiations and through the canopy air between them. This system is directly inspired from the ISBA-MEB (Boone et al., 2017) and FSM (Essery et al., 2025) models, from which the soil component has been omitted for simplicity.

## 4.1 Continuous domain

A schematic of the considered system is displayed in Fig. 2. Following Boone et al. (2017), the canopy is considered as having a single temperature, while the snowpack has a uni-dimensional spatially varying temperature. Finally, the canopy air is considered as having a zero thermal capacity. Thus, the system is notably described by the two energy conservation equations in the snowpack and in the canopy. For the snowpack, we have

$$ \partial_t u_{sp} + \nabla \cdot (F_{sp}) = Q_{sp} \quad \text{in } \Omega_{sp} \qquad (13) $$
with $u_{sp}$ the internal energy density of the snowpack, $F_{sp}$ the internal heat flux in the snowpack, $Q_{sp}$ the external energy source (i.e. shortwave absorption), and $\Omega_{sp}$ the spatial domain of the snowpack. This energy budget is complemented with boundary fluxes: $F_{sc}$ the net energy flux between the snow surface and the canopy, $F_{sa}$ the net energy flux between the snow surface and the canopy air, and $F_{bot}$ the energy flux at the bottom of the snowpack (all in W m$^{-2}$). Similarly, for the thermally homogeneous canopy we have

285

$$ \partial_t U_c = F_{sc} - F_{ca} + Q_c \qquad (14) $$

where $U_c$ is the internal energy of the canopy (directly expressed in J m$^{-2}$), $F_{ca}$ the net energy flux between the canopy and the canopy air, and $Q_c$ an external energy source (i.e. shortwave absorption). Note that because we consider a thermally homogeneous canopy, its energy conservation equation is already integrated over the vertical canopy space, and thus $U_c$ is not an internal energy density. Finally, the zero thermal capacity of the canopy air implies that the sum of the energy fluxes towards of the canopy air vanishes, i.e.

290

$$ F_{sa} + F_{ca} + F_{atm} = 0 \qquad (15) $$

295

where $F_{atm}$ is the energy flux coming from the external atmosphere orientated towards the canopy air.

By multiplying Eqs 13 and 14 by $\frac{1}{T_a}$ and $\frac{1}{T_c}$ respectively and using Gibbs' equations for snow and canopy (which are directly similar to Eq. (5)), we can recover equations for the snowpack entropy density $s_{sp}$ and canopy entropy $S_c$ (following the same procedure as in Sect. 3). The budget equation of the total entropy $S$ of the system is thus

300

$$ \begin{aligned} \partial_t S = \partial_t S_c + \int_{\Omega_{sp}} \partial_t s_{sp} d\Omega_{sp} \\ = \int_{\Omega_{sp}} F_{sp} \cdot \nabla \left( \frac{1}{T_{sp}} \right) d\Omega_{sp} \\ + F_{sc} \left( \frac{1}{T_c} - \frac{1}{T_{sp}^{surf}} \right) + F_{sa} \left( \frac{1}{T_a} - \frac{1}{T_{sp}^{surf}} \right) + F_{ca} \left( \frac{1}{T_a} - \frac{1}{T_c} \right) \\ + F_{atm} \left( \frac{1}{T_a} - \frac{1}{T_{atm}} \right) + F_{bot} \left( \frac{1}{T_{sp}^{bot}} - \frac{1}{T_{ground}} \right) \\ + \int_{\Omega_{sp}} \frac{Q_{sp}}{T_{sp}} d\Omega_{sp} + \frac{Q_c}{T_c} \\ + \frac{F_{atm}}{T_{atm}} + \frac{F_{bot}}{T_{ground}}, \end{aligned} \qquad (16) $$

with $T_{sp}^{surf}$ and $T_{sp}^{bot}$ the surface and bottom temperatures of the snowpack, $T_{atm}$ and $T_{ground}$ the temperatures of the ground below the snowpack and of the atmosphere, and where we have implicitly used Eq. (15). The temporal variation of the total entropy has thus five different origins: (i) an entropy production due to the heat fluxes within the snowpack, (ii) an entropy production due to the heat flux between the sub-parts of the system (e.g. $F_{sc}(\frac{1}{T_c} - \frac{1}{T_{sp}^{surf}})$ between the snowpack surface

305
310

and the canopy), (iii) an entropy production at the boundaries of the system due to energy exchange with the exterior (e.g. $F_{atm}$ ($\frac{1}{T_a}$ - $\frac{1}{T_{atm}}$) between the atmosphere and the canopy air), (iv) an external entropy source due to shortwave absorption, and (v) boundary entropy fluxes due to energy exchange at the boundaries of the system (e.g. $\frac{F_{atm}}{T_{atm}}$ the entropy flux at the atmosphere boundary).

315

The second principle of thermodynamics requires the three entropy production terms to be non-negative. This is readily the case when (i) the snowpack follows Fourier's law (in a linear or non-linear fashion; as in Sect. 3.1.4), and (ii) the energy fluxes between the different subcomponents of the system and the heat fluxes at the boundaries are consistent with the temperatures at play, i.e. that the energy fluxes are orientated from the warmer part towards the colder. These two points are of similar nature, and simply state that energy fluxes must be consistent with the temperature gradients/differences in the system that drive them.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Schematic of the configuration of the considered Snow-Canopy coupled model. The simulated system of interest is encapsulated in the dashed rectangle, and is composed a spatially in-homogenous snowpack (blue), a homogeneous canopy (green), exchanging energy through longwave radiation (red arrows) and through a canopy air space (orange dot) connected to the outside atmosphere. The outside external conditions (i.e. the temperatures of the atmosphere and below the snowpack) are assumed to be given.</figcaption></figure>
## 4.2 Numerical discretization

A similar analysis can be carried out in the discrete framework. Starting from the discrete energy conservation equations of the snowpack and the canopy (not shown for brevity since similar to Eq.(10)), one can recover a discrete equation for the total entropy variation by multiplying each cell equation by the inverse temperature at the end of the timestep $\frac{1}{T^{n+1}}$, summing over all cells, and re-arranging

320

$$ \frac{S^{n+1} - S^n}{\Delta t} = \mathcal{A}^{n+1} + \mathcal{B}^{n+1} + \mathcal{C}^{n+1} + \mathcal{D}^{n+1} + \Psi_{\text{ext}}^{n+1} + \mathcal{F}^{n+1}, \quad (17) $$

where we have set

$$ \mathcal{A}^{n+1} = \frac{1}{\Delta t} \sum_{\text{snow cells k}} \Delta z_k \left( s_k^{n+1} - s_k^n - \frac{u_k^{n+1} - u_k^n}{T_k^{n+1}} \right) + \frac{1}{\Delta t} \left( S_c^{n+1} - S_c^n - \frac{U_c^{n+1} - U_c^n}{T_c^{n+1}} \right), \quad (18) $$

325

$$ \mathcal{B}^{n+1} = \sum_{\text{internal snow edges}} F_{k+\frac{1}{2}}^{n \to n+1} \left( \frac{1}{T_{k+1}^{n+1}} - \frac{1}{T_k^{n+1}} \right), \quad (19) $$

$$ C^{n+1} = F_{\text{sc}}^{n \to n+1} \left( \frac{1}{T_c^{n+1}} - \frac{1}{T_{\text{surf}}^{n+1}} \right) + F_{\text{sa}}^{n \to n+1} \left( \frac{1}{T_a^{n+1}} - \frac{1}{T_{\text{surf}}^{n+1}} \right) + F_{\text{ca}}^{n \to n+1} \left( \frac{1}{T_a^{n+1}} - \frac{1}{T_c^{n+1}} \right), \quad (20) $$

330

$$ \mathcal{D}^{n+1} = F_{\text{atm}}^{n \to n+1} \left( \frac{1}{T_a^{n+1}} - \frac{1}{T_{\text{atm}}} \right) + F_{\text{bot}}^{n \to n+1} \left( \frac{1}{T_{\text{bot}}^{n+1}} - \frac{1}{T_{\text{ground}}} \right), \quad (21) $$

$$ \Psi_{\text{ext}}^{n+1} = \sum_{\text{snow cells k}} \frac{\Delta z_k Q_k^{n \to n+1}}{T_k^{n+1}} + \frac{Q_c^{n \to n+1}}{T_c^{n+1}}, \quad (22) $$

and

$$ \mathcal{F}^{n+1} = \frac{F_{\text{bot}}^{n \to n+1}}{T_{\text{ground}}} + \frac{F_{\text{atm}}^{n \to n+1}}{T_{\text{atm}}}. \quad (23) $$

335

As in the previous section, $z_k$, $u_k$, $s_k$, $T_k$, and $Q_k^{n \to n+1}$ stand for the thickness, internal energy density, entropy density, temperature, and absorbed energy of the k-th cell of the snowpack, and $F_{k+\frac{1}{2}}^{n \to n+1}$ for the internal heat flux between cells k and k+1 in the snowpack.

The terms $\Psi_{\text{ext}}^{n+1}$ and $\mathcal{F}^{n+1}$ correspond to the external entropy sources and boundary entropy fluxes, respectively. For the discrete entropy production, we see that it can be decomposed into four different contributions, denoted $\mathcal{A}^{n+1}$, $\mathcal{B}^{n+1}$, $\mathcal{C}^{n+1}$, $\mathcal{D}^{n+1}$ in Eq (17). As in Sect. 3.2, the term $\mathcal{A}^{n+1}$ is a purely discrete term, without any equivalent in the continuous domain.

340
It can be shown to be positive thanks to the concavity of both the snow and canopy entropy functions. Then, as in the single snowpack case, the term $\mathcal{B}^{n+1}$ corresponds to the discretized version of the internal snowpack entropy source, and is positive if the snowpack internal heat flux directions are consistent with the snowpack temperatures at the end of the time-step. The terms $C^{n+1}$ and $\mathcal{D}^{n+1}$ are of similar nature, and correspond to the entropy production due to energy fluxes between the different subcomponents of the system and between the system and the exterior. Again, they are positive if the heat fluxes between the subcomponents and the exterior are consistent with the temperature driving them.

## 4.3 The problem with a sequential treatment

In order to take advantages from modularity, the numerical simulation of a multi-components system is usually achieved by coupling pre-existing models developed for each subcomponent. In the simplest case, this coupling can be done with a fully sequential treatment. Such treatment for the snowpack-canopy system could for instance be: (i) the canopy temperature is first solved assuming fixed snowpack and canopy air temperatures from the end of the previous timestep. Note that for this stage, one needs to compute the energy flux between the canopy and the rest of the system (i.e. $F_{sc}^{n \to n+1}$ and $F_{ca}^{n \to n+1}$). Then, (ii) the snowpack temperature is solved using a snowpack model and a fixed canopy air temperature. To ensure energy conservation, the previously computed snowpack-canopy energy flux $F_{sc}^{n \to n+1}$ is used as such and directly injected as a Neumann boundary condition at the top of the snowpack. For this step, the flux $F_{sa}^{n \to n+1}$ needs to be computed. Finally, (iii) the canopy air temperature is solved, using the previously computed energy fluxes ($F_{sa}^{n \to n+1}$ and $F_{ca}^{n \to n+1}$), to ensure energy conservation.

However, this simple numerical strategy is known to be prone to numerical instabilities (personal communications with M.

Lafaysse and G. Mazzotti, June 2025). This limitation can be easily analyzed from an entropy point of view. While the use of an entropy consistent snowpack model ensures that the $\mathcal{B}^{n+1}$ term in Eq. (17) is positive, entropy destruction can still occur through the $C^{n+1}$ term. Indeed, as the energy fluxes between components are computed (and fixed) without taking into account the variations of all temperatures in the system, one cannot ensure the consistency between the energy flux directions and the temperatures at the end of the timestep. This can lead to entropy destruction at the interface between the different components. This problem notably arises in the case of small layers in the snowpack discretization (for instance related to freezing rain events), which are prone to large temperature fluctuations due to their smaller thermal capacity. In other words, the coupling of internally entropy-consistent sub-models does not imply that the resulting model is itself entropy-consistent. There are entropy sources related to exchange of the energy between the sub-models, which are directly controlled by none of them, but rather by the method of coupling itself.

On the other hand, the use of a tightly-coupled (with backward Euler temporal integration) scheme, where all temperatures are solved at once, ensures the unconditional consistency between all internal fluxes and the temperature field. In such a case, the discrete entropy source is positive everywhere and all the time, yielding a more robust model which can be run without instabilities. This is notably close to the strategy adopted in the numerical resolution of ISBA-MEB (Boone et al., 2017), where
375

the variation of the snowpack surface temperature is taken into account when solving the canopy energy budget, thanks to a
technique related to Schur-complement decomposition (Zhang, 2006).

4.4 A numerical example

380

To illustrate the analysis carried out for the simulation of the snowpack-canopy continuum, we have implemented two simple
march in time strategies: one using a sequential treatment and one solving the system of equations at once (i.e. monolithic or
tightly-coupled). The thermal capacity of the canopy is fixed to 10⁴ J K⁻¹. The snowpack is 50 cm thick with a refined grid at
the top (1 cm thickness for the uppermost cell), with a density ranging from 300 to 100 kg m⁻³, a constant thermal conductivity
of 0.2 W m⁻¹ K⁻¹, and without heat fluxes at its bottom. The energy fluxes F_ca, F_sa, and F_atm are assumed to be controlled by
constant linear thermal conductances of 1 W m⁻² K⁻¹, the flux F_sc is taken as σ (T_sp surf 4 – T_c 4) (with σ the Stefan-Boltzmann
constant), and shortwave absorption is neglected. Finally, the external atmosphere temperature is fixed to 270.15 K, which thus
corresponds to the equilibrium temperature towards which the system converge over time.

385

Simulation results with a timestep of 3600 s are displayed in Fig. 3. For simplicity, we focus here on the snowpack surface
and the canopy air only, without considering the canopy temperature, as it would not change the discussion of the results.
As seen in the lower panels of the figure, the tightly-coupled scheme exhibits a stable and physically consistent behavior,
slowly converging towards a steady state. This slow convergence towards equilibrium is due to the long relaxation time of the
snowpack, where heat must slowly diffuse over the whole height of the snowpack. A similar stable behavior has been observed
for all investigated timesteps and does not break in the case of a vanishingly small spatial discretization at the surface of the
snowpack. As displayed in the figure, this stable behavior is associated with a strictly positive entropy source, consistent with
the net direction of the energy flux between the snowpack and the canopy air: the colder snowpack surface receives energy from
the warmer canopy air. On the other hand, as seen in the upper panels of Fig. 3, the sequential scheme displays instabilities that
blow off over time. These instabilities are associated with a negative entropy source (i.e. entropy destruction), resulting from
the inconsistency between the orientation of the net snowpack-canopy air energy flux and the temperatures. It is interesting to
note that we have observed such instabilities even for a vanishingly small timestep of 1 ns. Our understanding is that since, the
canopy air has no thermal inertia, its characteristic time is essentially null. Therefore, its temperature can swing instantaneously
and remains prone to instabilities even for very small timesteps. On the other hand, a form of stability (in the sense that the
simulation does not blow off over time, despite the presence of spurious oscillations) is recovered in the case of very long
timesteps, in our case longer than 25 days. This illustrates that the intuitive notion that smaller timesteps lead to increased
stability is not always correct. Also, this shows that the possibility of a sequential resolution to be instable does not imply that
such instabilities are necessarily present. It is possible for a given sequential model to have a minimal timestep for stability that
is below the timestep of operation. However, special care should be given to ensure that this stability limit is not crossed.

390

395

400
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Simulation results of a snowpack-canopy system using a sequential (upper panels) or tightly-coupled (lower panels) scheme for a timestep of 3600 s. In both cases, the snowpack surface and canopy air temperatures are shown in blue and orange, respectively. The associated blue arrows indicate whether, during a given timestep, the snowpack surface receives (if orientated upward) or loses (if orientated downward) energy to the canopy air. Orange arrows hold the same meaning from the point of view of the canopy air. The entropy source related to the snowpack-canopy air energy flux is displayed below in green, with the zero level indicated by a dotted horizontal line.</figcaption></figure>

405 **5 A thermodynamically consistent model of dry snow with water vapor**

In this third part, we propose to expand on the simple heat equation studied in Sect. 3, by considering a thermodynamically
consistent numerical model of dry snow including macroscopic transport of water vapor. The goal of this analysis is to support
the implementation of water vapor transport in models (Jafari et al., 2020; Schürholt et al., 2022; Brondex et al., 2023). As
in the heat equation case, we start by analyzing the model in the continuous domain, starting from CIT frameworks. We then
investigate which numerical schemes preserve the entropy structure.

410
## 5.1 The continuous model

The derivation of a thermodynamically consistent model follows the overall steps presented in the case of the heat equation. Namely, we first need to define the entropy of snow and its relations to the prognostic extensive variables. In other words, we need to derive the Gibbs' equation applying to the ice and water vapor mixture. This Gibbs' equation can then be used in combination with the prognostic equations of energy and mass conservation to derive the entropy budget equation, and hence the entropy source. From it, conditions for its positivity can be analyzed.

The main difficulty, compared to the previous section, is that we are now concerned with a multiphase system described by more than one prognostic variable. The derivation from scratch of a concave entropy function for snow that encompasses the expected thermal and chemical material laws and respects Schwarz's theorem (i.e. Maxwell relations) can be quite challenging. This was notably illustrated by several failed attempts on our side to directly derive such a function relating the entropy of snow to other extensive variables. Rather, we have found that an easier, and more powerful path, is to first thermodynamically describe each homogeneous phase (in our case the solid ice and the gas phase) and then to build the entropy of the multiphase medium using the additivity of entropy.

### 5.1.1 Thermodynamical description of the ice phase

Let us consider a finite volume $V_i$ of ice embedded in a volume $V$ of snow. Thermodynamically, the ice phase is characterized by its internal energy density $u_i$ (J m⁻³ of snow), its volume fraction $\phi_i$ (m³ m⁻³), its density of matter $n_i$ (mol m⁻³), and its entropy density $s_i$ (J K⁻¹ m⁻³). Moreover, for this homogeneous phase, Gibbs' equation writes

$$d s_i = \frac{1}{T_i} d u_i + \frac{p_i}{T_i} d \phi_i - \frac{\mu_i}{T_i} d n_i \qquad (24)$$

with $T_i$, $p_i$, and $\mu_i$ the temperature (K), pressure (Pa), and chemical potential of the ice (J mol⁻¹), respectively (Kardar, 2007). Gibbs' equation directly implies that $s_i = s_i(u_i, v_i, n_i)$. Using the extensivity of $s_i$ we also have that

$$s_i = n_i s_i \left( \frac{u_i}{n_i}, \frac{\phi_i}{n_i}, \frac{n_i}{n_i} \right) = n_i \tilde{s}_i (\tilde{u}_i, \tilde{v}_i) \qquad (25)$$

with $\tilde{s}_i$ the molar entropy of pure ice, $\tilde{u}_i$ its molar internal energy, and $\tilde{v}_i$ its molar volume. Moreover, $\tilde{s}_i$ follows a Gibbs' equation with

$$d \tilde{s}_i = \frac{1}{T_i} d \tilde{u}_i + \frac{p_i}{T_i} d \tilde{v}_i. \qquad (26)$$

We note that Eq. (25) implies that the concavity of $s_i(u_i, v_i, n_i) = n_i \tilde{s}_i (\tilde{u}_i, \tilde{v}_i)$ directly follows that of the entropy density $\tilde{s}_i$ (a demonstration is available in Appendix A). Moreover, the extensivity of $u_i$ yields the Gibbs-Duhem relationship

$$\tilde{s}_i d T_i - \tilde{v}_i d p_i + d \mu_i = 0 \qquad (27)$$

suggesting that in the ice phase, the intensive quantities $T_i$, $p_i$, and $\mu_i$ are directly related.
440 The specific expression of the function $\tilde{s}_i(\tilde{u}_i, \tilde{v}_i)$ defines the material (in this case thermal, mechanical, and chemical) properties of the ice phase. Any choice could be made here, as long as the function $\tilde{s}(\tilde{u}_i, \tilde{v}_i)$ is concave. For the rest of the article, we assume that the ice is an incompressible solid of molar volume $\tilde{v}_i^0$ and with a constant molar thermal capacity $c_i$. Under these assumptions, the molar entropy of the ice phase should be taken as

$$ \tilde{s}_i(\tilde{u}_i, \tilde{v}_i) = \tilde{s}_i^0 + c_i \ln \left( 1 + \frac{\tilde{u}_i - \tilde{u}_i^0}{c_i T_0} \right) + \chi_{\tilde{v}_i^0}(\tilde{v}_i) \quad (28) $$

445 where $\tilde{u}_i^0$ and $\tilde{s}_i^0$ are integration constants, and $\chi_{\tilde{v}_i^0}$ the characteristic function of the set $\{\tilde{v}_i^0\}$, defined such that

$$ \chi_{\tilde{v}_i^0}(\tilde{v}_i) = \begin{cases} 0 & \text{if } \tilde{v}_i = \tilde{v}_i^0, \\ -\infty & \text{otherwise.} \end{cases} \quad (29) $$

The role of the characteristic function in Eq. (29) is to impose the incompressibility condition for ice and its constant molar volume. It also plays a role in defining the pressure of an incompressible solid. Indeed, deriving the $\tilde{s}_i$ with respect to $\tilde{u}_i$ and $\tilde{v}_i$ yields the expressions of the temperature $T_i$ and pressure $p_i$ as

$$ 450 \qquad T_i = \left( \frac{d\tilde{s}_i}{d\tilde{u}_i} \right)^{-1} = T_0 + \frac{\tilde{u}_i - \tilde{u}_i^0}{c_i} \qquad (30) $$

and

$$ \frac{p_i}{T_i} \in \partial_{\tilde{v}_i} \tilde{s}_i = \partial_i \chi_{\tilde{v}_i^0}(\tilde{v}_i) = \begin{cases} \mathbb{R} & \text{if } \tilde{v}_i = \tilde{v}_i^0, \\ \emptyset & \text{otherwise,} \end{cases} \quad (31) $$

with $\partial_i \chi_{\tilde{v}_i^0}$ denoting the superdifferential of the concave and lower semi-continuous function $\chi_{\tilde{v}_i^0}$ (see for instance Clarke, 1990; Rockafellar and Wets, 1998). For interested readers who are not familiar with the notion of sub- or super-differential, as brief introduction to this theoretical framework is provided in Appendix C. As the superdifferential of the characteristic function $\chi_{\tilde{v}_i^0}$ is multivalued, the pressure $p_i$ cannot be deduced from the mere Eq. (31). From a physical point of view, it encodes the fact the pressure of an incompressible solid is not fully determined by its internal properties (i.e. energy density and molar volume), but require extra information (typically an external pressure) to be fixed. Similarly, the chemical potential of the ice phase can be derived as

$$ 460 \qquad \mu_i = \tilde{u}_i + p_i \tilde{v}_i - T_i \tilde{s}_i \in -T_i \partial_{n_i} s_i, \qquad (32) $$

showing that there is a direct link between $T_i$, $p_i$, and $\mu_i$, in accordance with Gibbs-Duhem. Specifically, for incompressible ice, it implies that the chemical potential or pressure cannot be determined without fixing one or the other.

### 5.1.2 Thermodynamical description of the gas phase

The thermodynamical description of the gas phase shares many features with that of the solid ice phase. For simplicity, we consider the gas phase to be composed of water vapor only, neglecting the presence of dry air. As for the ice phase, the internal
energy density $u_v$, volume fraction $\phi_v$, matter density $n_v$, and entropy density $s_v$ are related through Gibbs' equation

$$ds_v = \frac{1}{T_v} du_v + \frac{p_v}{T_v} d\phi_v - \frac{\mu_v}{T_v} dn_v, \qquad (33)$$

where we have introduced the temperature, pressure, and chemical potential of the water vapor. Using the extensivity of $u_v$, the Gibbs-Duhem relationship reads

$$\tilde{s}_v dT_v - \tilde{v}_v dp_v + d\mu_v = 0. \qquad (34)$$

Similarly, using the extensivity of $s_v$ we also have

$$s_v = n_v s_v \left( \frac{u_v}{n_v}, \frac{\phi_v}{n_v}, \frac{n_v}{T_v} \right) = n_v \tilde{s}_v (\tilde{u}_v, \tilde{v}_v), \qquad (35)$$

with $\tilde{s}_v$, $\tilde{u}_v$, and $\tilde{v}_v$ being molar quantities, related through

$$d\tilde{s}_v = \frac{1}{T_v} d\tilde{u}_v + \frac{p_v}{T_v} d\tilde{v}_v. \qquad (36)$$

To specify the function $\tilde{s}_v$, we assume water vapor to be an ideal gas with a constant molar isochoric thermal capacity $c_v$. Under these conditions we have

$$\tilde{s}_v (\tilde{u}_v, \tilde{v}_v) = \tilde{s}_v^0 + c_v \ln \left( 1 + \frac{\tilde{u}_v - \tilde{u}_v^0}{c_v T_0} \right) + R \ln \left( \frac{\tilde{v}_v}{\tilde{v}_v^0} \right) \qquad (37)$$

where $\tilde{s}_v^0$, $\tilde{u}_v^0$, and $\tilde{v}_v^0$ are constants of integrations to be adjusted in what follows. By deriving $\tilde{s}_v$ with respect to $\tilde{u}_v$ and $\tilde{v}_v$ we indeed recover that

$$T_v = T_0 + \frac{\tilde{u}_v - \tilde{u}_v^0}{c_v} \qquad (38)$$

and

$$p_v \tilde{v}_v = R T_v. \qquad (39)$$

The chemical potential of water vapor is then given by

$$\mu_v = T_v \frac{\partial s_v}{\partial n_v} = \tilde{u}_v + p_v \tilde{v}_v - T_v \tilde{s}_v. \qquad (40)$$

### 5.1.3 Defining phase equilibrium

The formulations for the ice and water vapor entropy proposed above involve integration constants. The setting of these integration constants defines the internal energy, enthalpy, and entropy differences between the two phases. In the process, this specifies the conditions of equilibrium between the two phases.

For this, we first set the reference temperature to $T_0 = 273.15$ K, without loss of the generality. We then assume (still without loss of generality) that the molar internal energy and entropy of ice vanishes at $T_0$, yielding that $\tilde{u}_i^0 = 0$ and $\tilde{s}_i^0 = 0$ in Eq. (29).
This fully constrains the integration constants in the material description of the ice phase.

For the gas phase, we note $p_0 = 611$ Pa the saturation pressure of water vapor over ice at $T_0$ (Lide, 2005). We set the reference molar volume in Eq. (37) to the molar volume of the saturated water vapor at $T_0$, i.e. $\tilde{v}_0 = \frac{RT_0}{p_0}$. In these conditions, the integration constants $\tilde{u}_v^0$ and $\tilde{u}_v^0$ in Eq. (37) correspond to the internal energy of sublimation and entropy of sublimation at $T_0$ and $p_0$, respectively. These are directly related to the molar enthalpy of sublimation $\Delta\tilde{h}_{i\to v}$ at $T_0$ and $p_0$, with $\tilde{u}_v^0 = \Delta\tilde{h}_{i\to v} - p_0(\tilde{v}_v^0 - \tilde{v}_i^0)$ and $\tilde{s}_v^0 = \frac{\Delta\tilde{h}_{i\to v}}{T_0}$. This naturally enforces the equality of the chemical potentials, and hence phase equilibrium, in the reference conditions.

We note that the dependence of the saturation pressure of water vapor on temperature directly follows the definitions of the entropy proposed above, and does not require being further assumed. Indeed, following the classical derivation of Clausius equation (Kardar, 2007), the application of Gibbs-Duhem in both phases reads $\tilde{s}_i dT_i - \tilde{v}_i dp_i + d\mu_i = 0$ and $\tilde{s}_v dT_v - \tilde{v}_v dp_v + d\mu_v = 0$. At equilibrium, both phases have the same temperature $T$, pressure $p_v^{\text{sat}}$, and chemical potential, yielding

$$ \tilde{s}_i dT - \tilde{v}_i dp_v^{\text{sat}} = \tilde{s}_v dT - \tilde{v}_v dp_v^{\text{sat}}. \quad (41) $$

This readily yields Clausius equation, i.e.

$$ \frac{dp_v^{\text{sat}}}{dT} = \frac{\Delta\tilde{s}_{i\to v}}{\Delta\tilde{v}_{i\to v}} = \frac{\Delta\tilde{h}_{i\to v}}{T\Delta\tilde{v}_{i\to v}}. \quad (42) $$

Note that strictly speaking $\Delta\tilde{h}_{i\to v}$ is not a constant (it varies with temperature and pressure) and therefore, Clausius equation cannot be easily integrated to yield an analytical expression of $p_v^{\text{sat}}(T)$.

### 5.1.4 Thermodynamical description of snow

Once the ice and gas phases and their chemical equilibrium have been described, we can determine an expression for the entropy of the ice and water vapor mixture. Based on the additivity of entropy, the general expression for the expression of the ice + water vapor entropy $\mathcal{S}$ is

$$ \mathcal{S}(u_i, n_i, \phi_i, u_v, n_v, \phi_v) = s_i(u_i, n_i, \phi_i) + s_v(u_v, n_v, \phi_v). \quad (43) $$

To simplify this expression, we assume that (i) the volume of snow is fully filled, i.e. that $\phi_i + \phi_v = 1$ and that (ii) the ice and gas phase are in thermal equilibrium. Following standard statistical physics, the condition of thermal equilibrium is sought as the maximization of the entropy under constant total internal energy $u = u_i + u_v$ (Kardar, 2007). Hence, we have

$$ s(u, n_i, n_v) = \sup_{\substack{u=u_i+u_v \\ \phi_i+\phi_v=1}} [s_i(u_i, n_i, \phi_i) + s_v(u_v, n_v, \phi_v)], \quad (44) $$

with the supremum taken with respect to $u_i$, $\phi_i$, $u_v$, and $\phi_v$. Note that we do not have to enforce the molar volume of ice as a constraint in this maximization, as it is already included through $\chi_{\tilde{v}_i^0}$ in $s_i$. The above maximization problem is associated with
$$ \mathcal{L} = s_i(u_i, n_i, \phi_i) + s_v(u_v, n_v, \phi_v) + \tau(u - u_i - u_v) + \pi(1 - \phi_i - \phi_v) \quad (45) $$

where $\tau$ and $\pi$ are two Lagrange multipliers. The mixture entropy for snow defined by Eq. (44) is then, by construction, a jointly concave function of $(u, n_i, n_v)$, cf. Appendix B1. The maximization of $\mathcal{L}$ with respect to $u_i$ and $u_v$ leads to $T_i = T_v = 1/\tau = T$, i.e. thermal equilibrium as expected. Furthermore, the maximization of $\mathcal{L}$ with respect to $\phi_i$ and $\phi_v$ yields that $p_i = p_v = p = T\pi$, i.e. mechanical equilibrium. Moreover, we have that

$$ \begin{aligned} \partial_u s &= \tau = 1/T \\ \partial_{n_i} s &= \partial_{n_i} s_i = -\frac{\mu_i}{T} \\ \partial_{n_v} s &= \partial_{n_v} s_v = -\frac{\mu_v}{T} \end{aligned} \quad (46) $$

yielding the Gibbs equation for the snow material

$$ ds = \frac{1}{T} du - \frac{\mu_i}{T} dn_i - \frac{\mu_v}{T} dn_v. \quad (47) $$

Thanks to the material laws chosen for the ice and gas phases, namely Eqs. (29) and (37), it is possible to express the analytical formula for the snow entropy, that is

$$ s(u, n_i, n_v) = (n_i c_i + n_v c_v) \ln \left( 1 + \frac{u - n_v \Delta \tilde{u}_v^0}{T_0 (n_i c_i + n_v c_v)} \right) + n_v \left( \Delta s^0 + R \ln \left( \frac{1 - n_i \tilde{v}_i^0}{n_v \tilde{v}_v^0} \right) \right). \quad (48) $$

Such an entropy $s$ is concave in the broad sense, but is nowhere strictly concave owing to Appendix B2, leading to difficulties further discussed in Sect. 6.2.

Expression (48) could be derived with respect to $u$, $n_i$, $n_v$ to obtain the temperature and chemical potentials. There is however a less cumbersome way to derive the expression for $T$, $\mu_i$, and $\mu_v$. Using the additivity of internal energy we have $u = u_i + u_v = n_i c_i (T - T_0) + n_v (c_v (T - T_0) + \Delta \tilde{u}_v^0)$, which directly gives that

$$ T = T_0 + \frac{u - n_v \Delta \tilde{u}_v^0}{n_i c_i + n_v c_v} \quad (49) $$

showing that the snow + water vapor system as an effective volumetric thermal capacity of $n_i c_i + n_v c_v$, and that the reference temperature $T_0$ is reached for $u = n_v \Delta \tilde{u}_v^0$. The chemical potentials can be obtained through $\mu = \tilde{u} + p\tilde{v} - T\tilde{s}$, which applies for both the ice and the water vapor. In the gaseous phase, this yields

$$ \mu_v = c_v (T - T_0) + \Delta \tilde{u}_v^0 + RT - T(c_v + R) \ln \left( \frac{T}{T_0} \right) - T \Delta \tilde{s}_v^0 + RT \ln \left( \frac{p}{p_0} \right) \quad (50) $$

and

$$ \mu_i = c_i (T - T_0) + p \tilde{v}_i^0 - T c_i \ln \left( \frac{T}{T_0} \right), \quad (51) $$
where we used the fact the $p_i = p_v = p$ in accordance with mechanical equilibrium. Finally, by searching the conditions under which $\mu_v = \mu_i$, we note that the saturation pressure of water vapor $p_v^{\text{sat}}$ obeys the implicit equation

$$RT \ln \left( \frac{p_v^{\text{sat}}}{p_0} \right) - p_v^{\text{sat}} \tilde{v}_i^0 = (c_i - c_v)(T - T_0) + (c_v + R - c_i) T \ln \left( \frac{T}{T_0} \right) - RT + T \Delta \tilde{s}_v^0 - \Delta \tilde{u}_v^0 \quad (52)$$

which can be solved with the use of Lambert's W function. Indeed, following Eq. (52) we have

$$\ln \left( \frac{p_v^{\text{sat}} \tilde{v}_i^0 T_0}{p_0 \tilde{v}_v^0 T} \right) - \frac{p_v^{\text{sat}} \tilde{v}_i^0 T_0}{p_0 \tilde{v}_v^0 T} = \frac{A}{RT} + \ln \left( \frac{\tilde{v}_i^0 T_0}{\tilde{v}_v^0 T} \right) \quad (53)$$

with $A = (c_i - c_v)(T - T_0) + (c_v + R - c_i) T \ln \left( \frac{T}{T_0} \right) - RT + T \Delta \tilde{s}_v^0 - \Delta \tilde{u}_v^0$. Following Lambert's W function definition, we have

$$p_v^{\text{sat}} = -p_0 \frac{\tilde{v}_v^0 T}{\tilde{v}_i^0 T_0} W \left( -\exp \left( \frac{A}{RT} + \ln \left( \frac{\tilde{v}_i^0 T_0}{\tilde{v}_v^0 T} \right) \right) \right) \quad (54)$$

This expression is in general different from the Clausius-Clapeyron equation, i.e. $p_v^{\text{sat}} = p_0 \exp \left( \frac{\Delta \tilde{h}_v^0}{R} \left( \frac{1}{T_0} - \frac{1}{T} \right) \right)$. The latter can however be recovered under the assumptions that $\tilde{v}_i^0 \ll \tilde{v}_v$ and of constant molar enthalpy of sublimation, i.e. $\tilde{h}_v - \tilde{h}_i = \Delta \tilde{h}_v^0$, which requires $c_i = c_v + R$.

### 5.1.5 Prognostic equations

Under our assumptions (i.e. thermal equilibrium and ice incompressibility), snow can thus be characterized by three variables: its total internal energy $u$, its ice molar density $n_i$, and its vapor molar density $n_v$. The description of temporal evolution of such a snowpack thus requires the prognostic equations governing the evolution of these three variables. These equations simply corresponds to the internal energy, water vapor, and ice budgets. Following Vidal et al. (1994), they write

$$\left\{ \begin{array}{l} \partial_t u + \nabla \cdot (F_u + \tilde{h}_v F_v) = Q \\ \\ \partial_t n_v + \nabla \cdot F_v = Q_{i \to v} \\ \\ \partial_t n_i = -Q_{i \to v} \end{array} \right. \quad (55)$$

where $F_u$ and $F_v$ are the conduction heat flux and vapor flux in the snow, $Q$ energy absorbed from the exterior (for instance the shortwave absorption), and $Q_{i \to v}$ the rate of internal vapor sublimation. Note that this system of equations is consistent with the works considering macroscopic vapor transport in snow (Calonne et al., 2014; Bouvet et al., 2024). The fast-kinetics (i.e. diffusion limited as in Moyne et al., 1988; Hansen and Foslien, 2015) case can be recovered by defining the phase change term $Q_{i \to v}$ such as the time scale of relaxation towards chemical equilibrium is much faster than the other time scales involved in the problem.

### 5.1.6 The internal entropy production

Using the Eq. (47) and multiplying the system of Eqs. (55) by $1/T$, $-\mu_i/T$, and $-\mu_v/T$, we can proceed similarly as in Sect. 3.1.4 to derive the governing budget equation for entropy as

$$\partial_t s + \nabla \cdot (F_s) = \psi_s + \frac{Q}{T} \quad (56)$$
where the entropy flux is given by

$$F_s = \frac{F_u}{T} + \frac{\tilde{h}_v - \mu_v}{T} F_v = \frac{F_u}{T} + \tilde{s}_v F_v \quad (57)$$

and the internal entropy source by

$$\psi_s = F_u \cdot \nabla \left(\frac{1}{T}\right) + F_v \cdot \left(\tilde{h}_v \nabla \frac{1}{T} - \nabla \left(\frac{\mu_v}{T}\right)\right) + Q_{i \to v}(\mu_i - \mu_v). \quad (58)$$

The entropy source term related to the vapor flux $F_v$ can be further simplified (Vidal et al., 1994; Lebon et al., 2008). Indeed, the Gibbs-Helmotz equation yields that

$$\nabla \left( \frac{\mu_v}{T} \right) = - \frac{\tilde{h}_v}{T^2} \nabla T + \frac{[\nabla \mu_v]_T}{T} = \tilde{h}_v \nabla \left( \frac{1}{T} \right) + \frac{[\nabla \mu_v]_T}{T} \quad (59)$$

where $[\nabla \mu_v]_T = \nabla \mu_v - (\partial_T \mu_v) \nabla T$ is the gradient of vapor chemical potential at constant temperature. Application of Gibbs-Duhem formula for the gas phase at constant temperature yields $[\nabla \mu_v]_T = \tilde{v}_v \nabla p_v$. Hence, the entropy production rewrites

$$\psi_s = J_u \cdot \nabla \left(\frac{1}{T}\right) + \frac{\tilde{v}_v}{T} J_v \cdot \nabla p_v + Q_{i \to v}(\mu_i - \mu_v). \quad (60)$$

The second principle states that this entropy production term should be non-negative, and thus imposes constrains on how the heat flux, vapor flux, and sublimation rate reacts to imbalances in temperature, vapor pressure and chemical potential. We note that a sufficient condition for the second principle to be fulfilled is to have $F_u = -\lambda \nabla T$ (i.e. Fourier's law), $F_v = -D \rho_v \nabla \ln(p/p_0)$ (i.e. Darcy's law, which is equivalent to Fick's law in the isothermal case), and for $Q_{i \to v}$ to be a function of the chemical affinity $A = \mu_i - \mu_v$ such that $Q_{i \to v}(A = 0) = 0$, $Q_{i \to v}(A \ge 0) \ge 0$, and $Q_{i \to v}(A \le 0) \le 0$. These conditions simply state that heat flows from the warm to cold regions, vapor diffuses from regions of high pressure to regions of low pressures, and that sublimation (respectively deposition) occurs when the chemical potential of ice is larger (respectively smaller) than that of vapor.

We note that the heat and vapor fluxes could also include cross-terms (i.e. the Soret and Dufour effects; Mortimer and Eyring, 1980; Vidal et al., 1994). The second principle then would impose for these cross-fluxes to be relatively small compared to the principal Fourier and Darcean/Fickian fluxes. Concretely, it would impose upper bounds for the magnitudes of the Soret and Dufour coefficients (Vidal et al., 1994). Finally, we note that cross-terms cannot exist between the sublimation rate $Q_{i \to v}$ and the gradients of temperature or vapor pressure, as they do not share compatible symmetries (Vidal et al., 1994).

## 5.2 Thermodynamically consistent space-time discretization

Starting from thermodynamics first principles, we were able to derive a consistent model for snow viewed as an ensemble of ice and vapor in thermal, mechanical, but not chemical equilibrium, whose evolution in time do not destroy entropy. As with the heat equation, it remains to investigate which numerical schemes preserve this entropy structure, and thus the stability that comes with it.
Bearing in mind the studies carried out in Sects. 3 and 4, we will only consider a FVM spatial discretization with a Backward Euler temporal integration. The discretized internal energy, water vapor, and ice budgets write for each cell $k$ and time step $n \to n+1$

$$ 
\begin{cases}
\Delta z_k \frac{u_k^{n+1} - u_k^n}{\Delta t} + F_{u,k+\frac{1}{2}} - F_{u,k-\frac{1}{2}} + \tilde{h}_{v,k+\frac{1}{2}} F_{v,k+\frac{1}{2}} - \tilde{h}_{v,k-\frac{1}{2}} F_{v,k-\frac{1}{2}} = \Delta z_k Q_k, \\
\Delta z_k \frac{n_{v,k}^{n+1} - n_{v,k}^n}{\Delta t} + F_{v,k+\frac{1}{2}} - F_{v,k-\frac{1}{2}} = \Delta z_k Q_{k,i \to v}, \\
\Delta z_k \frac{n_{i,k}^{n+1} - n_{i,k}^n}{\Delta t} = -\Delta z_k Q_{k,i \to v},
\end{cases}
\quad (61) $$

with $F_{u,k+\frac{1}{2}}$ and $F_{v,k+\frac{1}{2}}$ being the heat and vapor fluxes between cell $k$ and $k+1$ (or potentially boundary heat and vapor fluxes if the cell is on the boundary), $\tilde{h}_{v,k+\frac{1}{2}}$ the molar enthalpy crossing the interface between cell $k$ and $k+1$ (or potentially the molar enthalpy crossing the boundary of the domain), and $Q_{k,i \to v}$ the sublimation rate in the $k^{th}$ cell. The temporal superscripts are not shown to lighten the notation, but as the scheme is implicit in time, all these quantities are understood to be based on time $t^{n+1}$. The discrete entropy balance can be derived by multiplying the rows of the system of Eqs. (61) by $1/T_k$, $-\mu_{v,k}/T_k$, and $-\mu_{i,k}/T_k$ (again, the temporal superscripts are not shown to lighten the notation) and adding them. After summing over all the cells, we obtain the entropy balance for the whole snowpack

$$ \frac{S^{n+1} - S^n}{\Delta t} = \mathcal{A}^{n+1} + \mathcal{B}^{n+1} + \mathcal{C}^{n+1} + \mathcal{D}^{n+1} + \mathcal{E}^{n+1} + \Psi_{\text{ext}}^{n+1} + \mathcal{F}^{n+1} \quad (62) $$

with

$$ \mathcal{A}^{n+1} = \frac{1}{\Delta t} \sum_{\text{cells}} z_k \left( s_k^{n+1} - s_k^n - \frac{1}{T_k} (u_k^{n+1} - u_k^n) + \frac{\mu_{v,k}}{T_k} (n_{v,k}^{n+1} - n_{v,k}^n) + \frac{\mu_{i,k}}{T_k} (n_{i,k}^{n+1} - n_{i,k}^n) \right), \quad (63) $$

$$ \mathcal{B}^{n+1} = \sum_{\text{cells}} z_k Q_{k,i \to v} \frac{\mu_{i,k} - \mu_{v,k}}{T_k}, \quad (64) $$

$$ C^{n+1} = \sum_{\text{internal edges}} F_{u,k+\frac{1}{2}} \left( \frac{1}{T_{k+1}} - \frac{1}{T_k} \right), \quad (65) $$

$$ \mathcal{D}^{n+1} = \sum_{\text{internal edges}} F_{v,k+\frac{1}{2}} \left( \frac{\tilde{h}_{v,k+\frac{1}{2}} - \mu_{v,k+1}}{T_{k+1}} - \frac{\tilde{h}_{v,k+\frac{1}{2}} - \mu_{v,k}}{T_k} \right) = \mathcal{D}_{\text{Darcy}}^{n+1} + \mathcal{D}_{\text{Num}}^{n+1}, \quad (66) $$

with

$$ \mathcal{D}_{\text{Darcy}}^{n+1} = \sum_{\text{internal edges}} F_{v,k+\frac{1}{2}} \frac{\tilde{v}_{v,k+\frac{1}{2}}}{T_{v,k+\frac{1}{2}}} (p_{v,k} - p_{v,k+1}), \quad (67) $$

and

$$ \mathcal{D}_{\text{Num}}^{n+1} = \sum_{\text{internal edges}} F_{v,k+\frac{1}{2}} \left( \tilde{s}_{k+1} - \tilde{s}_k + \frac{1}{T_{k+1}} \left( \tilde{h}_{v,k+\frac{1}{2}} - \tilde{h}_{v,k+1} \right) - \frac{1}{T_k} \left( \tilde{h}_{v,k+\frac{1}{2}} - \tilde{h}_{v,k} \right) - \frac{\tilde{v}_{v,k+\frac{1}{2}}}{T_{v,k+\frac{1}{2}}} (p_{v,k} - p_{v,k+1}) \right), $$
$$ \varepsilon^{n+1} = \sum_{\text{boundaries}} \left( F_{u,\text{ext}} \left( \frac{1}{T_k} - \frac{1}{T_{\text{ext}}} \right) + F_{v,\text{ext}} \left( \frac{\tilde{h}_{v,\text{bc}} - \mu_{v,k}}{T_k} - \frac{\tilde{h}_{v,\text{bc}} - \mu_{v,\text{ext}}}{T_{\text{ext}}} \right) \right) = \varepsilon_{\text{Fourier}}^{n+1} + \varepsilon_{\text{Darcy}}^{n+1} + \varepsilon_{\text{Num}}^{n+1}, \quad (69) $$

with

$$ \varepsilon_{\text{Fourier}}^{n+1} = \sum_{\text{boundaries}} F_{u,\text{ext}} \left( \frac{1}{T_k} - \frac{1}{T_{\text{ext}}} \right), \quad (70) $$

$$ \varepsilon_{\text{Darcy}}^{n+1} = \sum_{\text{boundaries}} F_{v,\text{ext}} \frac{\tilde{v}_{v,\text{bc}}}{T_{v,\text{bc}}} (p_{v,\text{ext}} - p_{v,k}), \quad (71) $$

$$ \varepsilon_{\text{Num}}^{n+1} = \sum_{\text{boundaries}} F_{v,\text{ext}} \left( \tilde{s}_k - \tilde{s}_{\text{ext}} + \frac{1}{T_k} (\tilde{h}_{v,\text{bc}} - \tilde{h}_{v,k}) - \frac{1}{T_{\text{ext}}} (\tilde{h}_{v,\text{bc}} - \tilde{h}_{v,\text{ext}}) - \frac{\tilde{v}_{v,\text{bc}}}{T_{v,\text{bc}}} (p_{v,\text{ext}} - p_{v,k}) \right), \quad (72) $$

$$ \Psi^{n+1} = \sum_{\text{cells}} \Delta z_k \frac{Q_k}{T_k}, \quad (73) $$

and

$$ \mathcal{F}^{n+1} = \sum_{\text{boundaries}} \left( \frac{F_{u,\text{ext}}}{T_{\text{ext}}} + F_{v,\text{ext}} \frac{\tilde{h}_v - \mu_{v,\text{ext}}}{T_{\text{ext}}} \right), \quad (74) $$

where $\tilde{v}_{v,k+1/2}$ and $T_{v,k+1/2}$ are the molar volume and temperature of the advected vapor, $F_{u,ext}$ and $F_{v,ext}$ are the heat and vapor fluxes at a boundary (orientated toward the snowpack), and $\tilde{h}_{v,bc}$, $\tilde{v}_{v,bc}$, and $T_{v,bc}$ are the molar enthalpy, molar volume, and temperature of the water vapor crossing the boundary. In Eq. (62), the terms $\mathcal{A}^{n+1}$, $\mathcal{B}^{n+1}$, $\mathcal{C}^{n+1}$, $\mathcal{D}^{n+1}$, $\varepsilon^{n+1}$ correspond to entropy production terms -and should thus be non-negative- while $\Psi_{\text{ext}}^{n+1}$ and $\mathcal{F}^{n+1}$ correspond to entropy exchange with the exterior, due to shortwave absorption and boundary fluxes respectively.

As in Sects. 3 and 4, $\mathcal{A}^{n+1}$ is a purely numerical entropy source, whose non-negativity follows from the concavity of the snow entropy. Indeed, the concavity of $s(u, n_i, n_v)$ yields

$$ \begin{aligned} \forall (u, n_i, n_v) : \quad & s(u, n_i, n_v) \le s(u_k^{n+1}, n_{v,k}^{n+1}, n_{i,k}^{n+1}) \\ & \qquad + \frac{\partial s}{\partial u} (u_k^{n+1}, n_{v,k}^{n+1}, n_{i,k}^{n+1}) (u - u_k^{n+1}) \\ & \qquad + \frac{\partial s}{\partial n_i} (u_k^{n+1}, n_{v,k}^{n+1}, n_{i,k}^{n+1}) (n_i - n_{i,k}^{n+1}) \\ & \qquad + \frac{\partial s}{\partial n_v} (u_k^{n+1}, n_{v,k}^{n+1}, n_{i,k}^{n+1}) (n_v - n_{v,k}^{n+1}) \\ & \le s(u_k^{n+1}, n_{v,k}^{n+1}, n_{i,k}^{n+1}) + \frac{1}{T_k} (u - u_k^{n+1}) - \frac{\mu_{i,k}}{T_k} (n_i - n_{i,k}^{n+1}) - \frac{\mu_{v,k}}{T_k} (n_v - n_{v,k}^{n+1}) \end{aligned} \quad (75) $$
which directly implies $\mathcal{A}^{n+1} \ge 0$ when evaluated in $(u_k^n, n_{1,k}^n, n_{v,k}^n)$. The terms $\mathcal{B}^{n+1}$ and $\mathcal{C}^{n+1}$ corresponds to the entropy production due to phase changes and internal heat fluxes terms, respectively. Their non-negativity follows from the choice of a Backward Euler temporal integration, which ensure that the sublimation rates $Q_{k,i \to v}$ and heat fluxes $F_{u,k+\frac{1}{2}}$ are consistent with the cell chemical affinities $\mu_{i,k} - \mu_{v,k}$ and with the temperature difference between cells $(\frac{1}{T_{k+1}} - \frac{1}{T_k})$.

A specificity of the term $\mathcal{D}^{n+1}$ is the presence of advected enthalpy with the vapor flux, which requires defining the value of the enthalpy transported across the internal edges. As seen above, $\mathcal{D}^{n+1}$ can be split into two distinct entropy production terms $\mathcal{D}_{\text{Darcy}}^{n+1}$ and $\mathcal{D}_{\text{Num}}^{n+1}$. $\mathcal{D}_{\text{Darcy}}^{n+1}$ corresponds to the discrete version of the entropy production due to the Darcean vapor flux. Its positivity follows the choice of a Backward Euler scheme, that ensures the consistence between the vapor flux and the pressure gradient that drives it. Notably, its sign is independent of the choice of the advected enthalpy $\tilde{h}_{v,k+\frac{1}{2}}$, as $\frac{\tilde{v}_{v,k+\frac{1}{2}}}{T_{v,k+\frac{1}{2}}} \ge 0$. On the other hand, $\mathcal{D}_{\text{Num}}^{n+1}$ is a purely numerical term, whose sign depends on the definition of the advected enthalpy $\tilde{h}_{v,k+\frac{1}{2}}$. In order for $\mathcal{D}_{\text{Num}}^{n+1}$ to be non-negative, the advected enthalpy should be upstreamed, i.e. taken from the cell giving the net flux of vapor. Indeed, let's assume, without loss of generality, that $p_{v,k} > p_{v,k+1}$ and thus that $F_{v,k+\frac{1}{2}} > 0$. The upstream rule yields $\tilde{h}_{v,k+\frac{1}{2}} = \tilde{h}_{v,k}$ and we can rewrite $\mathcal{D}_{\text{Num}}^{n+1}$ as

$$ \mathcal{D}_{\text{Num}}^{n+1} = \sum_{\text{internal edges}} F_{v,k+\frac{1}{2}} \left( \tilde{s}_{v,k+1} - \tilde{s}_{v,k} - \frac{1}{T_{k+1}} (\tilde{h}_{v,k+1} - \tilde{h}_{v,k}) + \frac{\tilde{v}_{k+1}}{T_{k+1}} (p_{v,k+1} - p_{v,k}) \right). \quad (76) $$

The concavity of $\tilde{s}_v$ and the identity $d\tilde{s}_v = \frac{d\tilde{h}_v}{T} - \frac{d\tilde{v}_v}{T} dp$ yields that

$$ \begin{aligned} \forall (\tilde{h}_v, p_v) : \quad & \tilde{s}_v(\tilde{h}_v, p_v) \le \tilde{s}_v(\tilde{h}_{v,k}, p_{v,k}) + \frac{\partial \tilde{s}_v}{\partial \tilde{h}_v}(\tilde{h}_{v,k}, p_{v,k}) (\tilde{h}_v - \tilde{h}_{v,k}) + \frac{\partial \tilde{s}_v}{\partial p}(\tilde{h}_{v,k}, p_{v,k}) (p_v - p_{v,k}) \\ & \le \tilde{s}_v(\tilde{h}_{v,k}, p_{v,k}) + \frac{1}{T_k} (\tilde{h}_v - \tilde{h}_{v,k}) - \frac{\tilde{v}_{v,k}}{T_k} (p_v - p_{v,k}), \end{aligned} \quad (77) $$

which evaluated at $(\tilde{h}_{v,k+1}, p_{v,k+1})$ proves the non-negativity of $\mathcal{D}_{\text{Num}}^{n+1}$.

Finally, the $\mathcal{E}^{n+1}$ term corresponds to the entropy produced at the boundary of the snowpack due to the exchanges with the exterior. It is composed of a heat flux term $\mathcal{E}_{\text{Fourier}}^{n+1}$, a term $\mathcal{E}_{\text{Darcy}}^{n+1}$ corresponding to Darcean vapor fluxes, and a purely numerical term $\mathcal{E}_{\text{Num}}^{n+1}$. For $\mathcal{E}^{n+1}$ to be non-negative, two conditions must be met. First, the boundary heat and vapor fluxes must be consistent with the temperature and pressure differences between the boundary of the snowpack and the exterior to ensure the non-negativity of $\mathcal{E}_{\text{Fourier}}^{n+1}$ and $\mathcal{E}_{\text{Darcy}}^{n+1}$. This is the case with a Backward Euler scheme, provided that the expressions for the boundary condition make physical sense. Second, the molar enthalpy $\tilde{h}_{v,\text{bc}}$ crossing the boundary must be upstreamed, i.e. taken from the exterior when the net vapor flux is orientated toward the snowpack and taken from the snowpack when the net vapor flux is orientated outward. Then, the non-negativity of the purely numerical contribution of $\mathcal{E}_{\text{Num}}^{n+1}$ can be shown following the same demonstration as for $\mathcal{D}_{\text{Num}}^{n+1}$.
Overall, we thus have that there are three criteria for this numerical model to be entropy-producing: (i) variables must be co-localized, (ii) fluxes/phase changes must be consistent with the gradient/differences that drive them, and (iii) advected vapor enthalpy must be upstreamed.

## 5.3 A numerical example

To illustrate the numerical behavior of a dry snow model with vapor diffusion, we have implemented a simple toy model respecting the criteria derived in the previous section, ensuring the respect of the second principle. Notably, all processes are solved at once using a Backward Euler temporal integration. This model is used to model the evolution of a shallow 40 cm snowpack with a dense crust in the center, subjected to a constant thermal gradient of 100 K m⁻¹ by imposing constant temperatures at the top and bottom (directly inspired by Schürholt et al., 2022). Corresponding vapor fluxes are imposed at the top and bottom to respect vapor flux continuity at the boundary of the snowpack. The vapor diffusivity and thermal conductivity of vapor explicitly depends on the ice volume fraction, following Calonne et al. (2011) and Calonne et al. (2014), respectively. The vapor diffusivities and thermal conductivities of internal interfaces are taken as the weighted harmonic average of the two cells around it (Kadioglu et al., 2008). Finally, the dynamics of the ice-vapor equilibrium was set as $Q_{i→v} = K_{depo}(\rho_v - \rho_v^{sat})$, with $\rho_v$ the vapor molar density in the pore space, $\rho_v^{sat}$ the saturated vapor molar density in the pore space, and $K_{depo}$ a constant set to 1 s⁻¹. In this configuration, the phase change dynamics is faster than the other processes (heat conduction and vapor diffusivity) in the snowpack, and thus rather corresponds to the diffusion-limited case for the vapor (similar to Moyne et al., 1988; Hansen and Foslien, 2015).

Results with a fine 0.8 mm mean cell thickness (corresponding to 500 cells over 40 cm) and with a large timestep of 6 days are displayed in Fig. 4 and for a total simulation durations of 600 days. The large time scales used in this simulation (near weekly timestep over 1.5 years of simulation) have been chosen to correspond to the time scale over which density evolves in response to vapor diffusion and phase changes (which is the motivation from the introduction of water vapor in snowpack models; Krinner et al., 2018; Domine et al., 2019). In agreement with Schürholt et al. (2022) and Brondex et al. (2023), Fig. 4 shows that the presence of a thermal gradient results in a downward movement of the density crust, with density oscillations forming above it. We stress that these oscillations are not numerical artifacts, but an actual physical pattern that develops due to the dependence of the vapor diffusivity and thermal conductivity to the ice volume fraction. Overall, this simulation highlights that the use of a second principle-complying model allows us to perform simulations with very large time steps, tailored to focus on the slow processes of interest, without numerical instabilities arising from the faster processes. In other words, the compliance to the second principle makes such a model robust against the numerical stiffness due to the large time scale differences between the various phenomena at play (Fazio, 2001). This possibility of performing simulations with large time steps is notably valuable for applications involving ice sheet snow and firn, which require simulations lasting for decades to millennia (Hoang et al., 2025).
However, despite its overall robustness, we observed clear numerical artifacts that developed in the case of a coarse and non-uniform spatial resolution (while numerical instabilities are usually expected to occur with fine meshes). The top left panel of Figure 5 shows the results of the same simulation as described above, but in the case of a coarser mesh (30 cells). As seen in the figure, an instability develops, which is likely related to the oscillations in the wake of the moving crust that tend to blow up with a coarse mesh. It can result in unphysical behavior, specifically here the apparition of negative ice volume fractions. We stress that these oscillations are in principle compatible with the second principle, as there are not associated with entropy destruction (all entropy sources being positive by construction of the model). As with the (physical) oscillations developing in the vicinity of the crust, these numerical artifacts only occur in the case of a dependence of the vapor diffusivity and thermal conductivity on the ice volume fraction. Also, the numerical artifacts cannot be resolved by simply choosing a shorter timestep. Therefore, they appear to be linked to the spatial discretization rather than the temporal one. After investigation, we observed that a strategy to avoid the occurrence of these artifacts is to use the downstream values for the interfaces' vapor diffusivity and/or thermal conductivity. As seen in the other panels of Fig. 5, this modification is sufficient to contain the oscillations developing in the vicinity of the density crust and to avoid unphysical ice volume fraction values. The strongest reduction of oscillations is obtained with the combination of both downstreaming the vapor diffusivity and thermal conductivity, where they almost entirely disappear. However, at the moment, this strategy is a simple heuristic, and the reasons why downstreaming the vapor diffusivity and/or thermal conductivity enhances stability remain unclear. This is all the more puzzling as the vapor diffusivity and thermal conductivity have opposite dependencies to the ice volume fraction.

# 6 Discussion

## 6.1 Extension to other cases

In this article, we considered three different cases involving snowpack modelling with different degrees of complexity. We thus limited this work to the study of snowpacks governed by heat fluxes and potentially vapor fluxes with additional phase changes. This is, however, only a fraction of the potential physical mechanisms at play in snow, and which could be accounted for in numerical models. Among the missing mechanisms/processes, we can mention the presence of liquid water, thus with phase changes and liquid water percolation under gravity and capillary effects, the mechanical compaction of the snowpack under its own weight, and the presence of water isotopologues. One can thus wonder about the extensibility of the approach and results presented in this work.

We believe that the methodology used in this work is quite general and should be directly transferable in the case of added mechanisms. We recall that it follows these general steps:

1. Physically describe each phase $\alpha$ of the snow material with an entropy function $s_{\alpha} = (\{x_{\alpha}\})$, where $\{x_{\alpha}\}$ is the set of the densities of extensive quantities characterizing the state of the phase. Typically, these are the internal energy density, volume fraction, and matter densities.

2. Derive the general entropy function of snow by adding the entropy contribution of each phase. Note that while this was not
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Simulation results showing the density evolution of a 40 cm snowpack subjected to a constant temperature gradient of 100 K m⁻¹ over 600 days. Dashed black line: initial density condition, with a density crust at mid-height. Solid color line: density profile over time, with a downward motion of the density crust and the apparition of density oscillations above.</figcaption></figure>

considered in this article, the introduction of capillary effects appears to be achievable at this stage by introducing an entropy of interaction between the phases (equivalent to defining an entropy of the interface between phases; Coussy, 2004).

3- Constraints, such as total volume filling or thermal equilibrium, can then be added by maximizing the entropy of snow under those constraints. This will reduce the number of independent prognostic variables (for instance down to total interval energy, ice content, and vapor content in the case studied in Sect. 5). This step will also define the Gibbs' equation of snow.

4- Define the equations of the prognostic variables (for instance total energy, ice content, and vapor content budgets).

5- Derive the (potentially discrete) entropy budget and its source, by combing Gibbs' equation and the (potentially discrete) prognostic equations.

6- Infer the conditions for the positive of the entropy source.

These are consistent with the general framework of CIT, which has been used to derive the entropy production of various physical systems (Vidal et al., 1994; Lebon et al., 2008). To the best of our knowledge, snow should be no exception to this, even if the intricacy of physical mechanisms (porous medium with heat fluxes, phase changes, vapor and liquid water fluxes, capillary effects, viscous compaction, etc) may render its description more complicated.

Moreover, the few criteria that we derived (co-localization, flux/potential consistency, and upstreaming of advected quantities)
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Illustration of potential numerical artifacts occurring in the case of a coarse and non-uniform spatial resolution for a 600 days long simulation. In each panel, the dashed gray line represents a reference solution, obtained with a fine resolution (500 cells of a mean 0.8 mm thickness) and the blue line a solution obtained with a coarser grid (30 cells of a mean 13 mm thickness). Columns: solutions obtained with a harmonic average (left) or downstream (right) scheme for the interface vapor diffusivities. Rows: solutions obtained with a harmonic average (top) or downstream (bottom) scheme for the interface thermal conductivities.</figcaption></figure>

are fully consistent with the standard linear stability and discrete maximum principle analyses (Butcher, 2008), as well as with heuristics developed for modelling (essentially that tighter coupling provides more robust models Keyes et al., 2013; Fourteau et al., 2024). Therefore, the adherence to co-localization, flux/potential consistency, and upstreaming of advected quantities appear to be general enough to ensure the respect of the second principle in most (if not all) numerical models.

## 6.2 Thermodynamical consistency is not sufficient for numerical stability

As discussed in Sect. 2.2, the compliance to the second principle can be used to formally demonstrate Lyapunov stability for some specific boundary conditions, but its generalization to arbitrary boundary conditions remains an open question. Yet, the local entropy production is commonly accepted to be satisfactory condition for ensuring the stability of numerical schemes, especially for hyperbolic problems (see for instance the seminal contribution by Godunov (1959) or the monograph by Godlewski and Raviart (1996)).
In the numerical example involving a snowpack with vapor diffusion and phase changes, the lack of strict concavity for the entropy highlighted in Appendix B2 revealed that numerical instabilities/artifacts might still exist despite the strict adherence to the second principle. Similar difficulties arise for instance in Faccanoni et al. (2025). In our setting, some specific down-streaming of the permeabilities and of the thermal conductivities is required to avoid spurious oscillations, cf. Fig. 5. We do not have so far a clear understanding on the reasons why such downwind choices provide enhanced stability. Relating this enhanced stability to additional entropy criteria as suggested by Bruining et al. (2003) and Eymard and Tillier (2007) should be addressed in future works.

## 6.3 Should models respect the second principle?

The theory and results presented in this work naturally begs the question of whether models should be expected to respect the second principle, in the way they are expected to follow the first principle, i.e. energy conservation. On one hand, the second principle clearly favors numerical stability, and its compliance in models can only increase their overall robustness. On the other hand, it appears that respecting the second principle is not directly equivalent with numerical stability. It is neither sufficient (models respecting it can still suffer from numerical artifacts if the entropy function is not strictly concave) nor necessary (small entropy inconsistencies do not always translate into clear numerical artifacts that can unequivocally be identified as unphysical). This is to be considered alongside the fact that implementing strictly entropy-consistent models might be cumbersome, as it requires specific material laws (as exemplified by the expression for the saturation pressure in Eq. 52) and tight-coupling between processes and subcomponents. As a consequence, we think that the compliance of models to the second principle of thermodynamics should rather be viewed as an overall guide and a tool, rather than a strict rule to adhere to. Specifically, in cases where numerical models are known to be unstable, it might be fruitful to study whether this instability can be linked to a thermodynamical inconsistency and which design practice can be followed to resolve it. Also, thermodynamics considerations can be used to design efficient coupling strategies, that are less complex than direct tight-coupling, but still ensure entropy non-destruction (e.g. Liu et al., 2021). This might be harnessed in the future to design simple but robust coupling strategies between pre-existing models.

## 7 Conclusions

This article proposes to apply the so-called entropy methods to the case of snowpack modelling. Essentially, entropy methods state that the compliance of a physical system, or physical model, to the second principle of thermodynamics (i.e. the non-destruction of entropy), favors stability. Thus, the goal of this article is to see under which conditions numerical snowpack models and simulation chains in which snowpack models intervene respect the second principle of thermodynamics. For this, we considered three different cases of increasing complexity, namely a simple dry snowpack controlled by the heat equation only, a coupled snowpack-canopy system, and a dry snowpack with macroscopic vapor fluxes and internal phase changes. In each case, we started by ensuring that the physical equations (before any spatial or numerical discretization) used to describe
the system were consistent with the second principle. This step was achieved following the framework of Classical Irreversible Thermodynamics. Then, we studied under which numerical schemes the corresponding numerical models still respected a discrete version of the second principle.

Our conclusion, based on the considered cases, is that they appear to be three general criteria to ensure that the second principle holds in numerical models: (I) the variables must be co-localized, i.e. the temperature, pressure, and chemical potentials in a given layer must be a function of the energy and matter content in that cell only, (II) advected quantities, such as the enthalpy advected with vapor fluxes, must be upstreamed, and (III) the generalized fluxes (i.e. energy fluxes, vapor fluxes, and chemical reaction rates) must be consistent with their driving forces (i.e. temperature, vapor pressure, and chemical potential differences) at the end of the timestep. Criterion (I) can easily be fulfilled by adopting the finite-volume method for the spatial discretization, as already done in most snowpack models, or by mass lumping techniques in the case of the finite-elements method. Similarly, the relatively common choice of a backward Euler time discretization in snowpack models plays an important role in fulfilling criterion (III). That being said, this might not be sufficient, and specific care must also be given to couplings, might they be the couplings between physical processes within the snowpack or the couplings between the snowpack and the rest of the simulation chain in the case of coupled models simulations. For instance, as shown in this article, the combination of two entropy-consistent sub-models do not necessarily yield an entropy-consistent model, as entropy inconsistencies might appear at the interface of the two sub-models, compromising the numerical stability of the ensemble.

As achieving proper thermodynamical consistency might be cumbersome, it might be more interesting to use the second principle as a designing tool when building numerical snow models, rather than a strict rule to adhere to in any circumstances. In other words, thermodynamical consistency should likely be respected when easily doable, and its violations should ideally be motivated by the knowledge that the corresponding entropy sink will remain small enough not to be able to compromise the stability of the numerical system. Similarly, the second principle can play the role of a diagnostic tool in the case of known numerical instabilities in a model or in a modeling chain. Such analysis can indicate the origin of the numerical problem (for instance the inconsistency between a flux and a potential field) and suggest an appropriate fix (for instance the tight coupling of two processes or subcomponents).

## Appendix A: Concavity of volumetric entropy density

The purpose of the Appendix is to show that the concavity of the volumetric entropy $s_\alpha$ of a phase $\alpha \in \{i, v\}$ (i.e. the entropy carried by the phase $\alpha$ per volume of snow), which is a 1-homogeneous (or extensive) function of $u_\alpha$, $n_\alpha$ and $\phi_\alpha$, is equivalent to the concavity of the intrinsic molar entropy of the phase $\tilde{s}_\alpha$ with respect to $\tilde{u}_\alpha = \frac{u_\alpha}{n_\alpha}$ and $\tilde{v}_\alpha = \frac{\phi_\alpha}{n_\alpha}$. Both are related by the formula $s_\alpha(u_\alpha, n_\alpha, \phi_\alpha) = n_\alpha \tilde{s}_\alpha(\tilde{u}_\alpha, \tilde{v}_\alpha)$ thanks to the 1-homogeneity of $s_\alpha$.

It is clear that if $s_\alpha$ is concave, so is $\tilde{s}_\alpha$ as the latter can be thought of as the restriction of $s_\alpha$ to the linear space $n_\alpha = 1$. To establish the reciprocal, let us assume that $\tilde{s}_\alpha$ is concave and set $f^0 = (u_\alpha^0, n_\alpha^0, \phi_\alpha^0)$ and $f^1 = (u_\alpha^1, n_\alpha^1, \phi_\alpha^1)$ two configurations for the energy, number of moles and volume fraction. Then for $\theta \in [0, 1]$, we denote $f^\theta = (u_\alpha^\theta, n_\alpha^\theta, \phi_\alpha^\theta)$ such that

$$f^\theta = \theta f^1 + (1-\theta)f^0.$$
Then the 1-homogeneity of $s_\alpha$ gives

$$ 
\begin{aligned}
s(u_\alpha^\theta, n_\alpha^\theta, \phi_\alpha^\theta) &= n_\alpha^\theta \tilde{s}_\alpha \left( \frac{u_\alpha^\theta}{n_\alpha^\theta}, \frac{\phi_\alpha^\theta}{n_\alpha^\theta} \right) = n_\alpha^\theta \tilde{s}_\alpha \left( \frac{\theta u_\alpha^1 + (1-\theta)u_\alpha^0}{n_\alpha^\theta}, \frac{\theta \phi_\alpha^1 + (1-\theta)\phi_\alpha^0}{n_\alpha^\theta} \right) \\
&= n_\alpha^\theta \tilde{s}_\alpha \left( \frac{\theta n_\alpha^1}{n_\alpha^\theta} \tilde{u}_\alpha^1 + \frac{(1-\theta)n_\alpha^0}{n_\alpha^\theta} \tilde{u}_\alpha^0, \frac{\theta n_\alpha^1}{n_\alpha^\theta} \tilde{\phi}_\alpha^1 + \frac{(1-\theta)n_\alpha^0}{n_\alpha^\theta} \tilde{\phi}_\alpha^0 \right).
\end{aligned}
$$

835 Noting that $\frac{\theta n_\alpha^1}{n_\alpha^\theta} \in [0,1]$ and that $\frac{\theta n_\alpha^1}{n_\alpha^\theta} = 1 - \frac{(1-\theta)n_\alpha^0}{n_\alpha^\theta}$, the concavity of $\tilde{s}_\alpha$ yields that

$$ s(u_\alpha^\theta, n_\alpha^\theta, \phi_\alpha^\theta) \geq \theta n_\alpha^1 \tilde{s}_\alpha(\tilde{u}_\alpha^1, \tilde{v}_\alpha^1) + (1-\theta)n_\alpha^0 \tilde{s}_\alpha(\tilde{v}_\alpha^0, \tilde{v}_\alpha^0) = \theta s_\alpha(u_\alpha^1, n_\alpha^1, \phi_\alpha^1) + (1-\theta)s_\alpha(u_\alpha^0, n_\alpha^0, \phi_\alpha^0), $$
completing the proof of the claim.

# Appendix B: Degenerate concavity of the mixture entropy of snow

## B1 The mixture entropy of snow is concave in general...

840 The general mixture entropy of snow $s(u, n_i, n_v)$, defined by Eq. (44), is by construction a jointly concave function of $u, n_i$, and $n_v$ as soon as the pure phase entropy functions $s_i$ and $s_v$ are concave as well. To establish this property, let us reformulate the problem as a saddle point problem:

$$ s(u, n_i, n_v) = \sup_{u_i, u_v, \phi_i, \phi_v} \inf_{\tau, \pi} \mathcal{L}((u, u_i, u_v, n_i, n_v, \phi_i, \phi_v); (\tau, \pi)), $$

where the Lagrangian

$$ \mathcal{L}((u, u_i, u_v, n_i, n_v, \phi_i, \phi_v); (\tau, \pi)) = s_i(u_i, n_i, \phi_i) + s_v(u_v, n_v, \phi_v) + \tau(u - u_i - u_v) + \pi(1 - \pi_i - \phi_v) $$

845 is jointly concave with respect to its 7 first arguments owing to the concavity of $s_i$ and $s_v$, and linear (thus convex) with respect to $\tau$ and $\pi$. Then owing to Theorem 11.39 in Rockafellar and Wets (1998), strong duality holds, so that one can swap the inf and the sup which turn out to be min and max, leading to

$$ s(u, n_i, n_v) = \min_{\tau, \pi} \max_{u_i, u_v, \phi_i, \phi_v} \mathcal{L}((u, u_i, u_v, n_i, n_v, \phi_i, \phi_v); (\tau, \pi)). \quad (B1) $$

850 Let $(\tau, \pi)$ be given, and let $(u^0, n_i^0, n_v^0)$ and $(u^1, n_i^1, n_v^1)$ be two different choices for the internal energy, the concentration of ice and vapor in the snow. Then for $\ell \in \{0,1\}$, we denote by $(u_i^\ell, u_v^\ell, \phi_i^\ell, \phi_v^\ell)$ the optimizer in Eq. (B1), that is such that

$$ \Lambda((u^\ell, n_i^\ell, n_v^\ell), (\tau, \pi)) = \max_{u_i, u_v, \phi_i, \phi_v} \mathcal{L}((u^\ell, u_i, u_v, n_i^\ell, n_v^\ell, \phi_i, \phi_v); (\tau, \pi)) = \mathcal{L}((u^\ell, u_i^\ell, u_v^\ell, n_i^\ell, n_v^\ell, \phi_i^\ell, \phi_v^\ell); (\tau, \pi)) \quad (B2) $$

Let us show that the above function $\Lambda$ is concave with respect to its first 3 arguments. To this end, we denote by

$$ f^\theta = \theta f^1 + (1-\theta)f^0. $$
855 for $\theta \in [0, 1]$ and $f \in \{u, u_i, u_v, n_i, n_v, \phi_i, \phi_v\}$. Note that $(u_i^\theta, u_v^\theta, \phi_i^\theta, \phi_v^\theta)$ is not in general the optimizer in Eq. (B1) corresponding to $(u^\theta, n_i^\theta, n_v^\theta)$ if $0 < \theta < 1$. Yet, the following inequality readily follows from Eq. (B2):

$$ \Lambda((u^\theta, n_i^\theta, n_v^\theta), (\tau, \pi)) \geq \mathcal{L}((u^\theta, u_i^\theta, u_v^\theta, n_i^\theta, n_v^\theta, \phi_i^\theta, \phi_v^\theta); (\tau, \pi)). $$

As the Lagrangian is concave with respect to its first variables, one deduces from the previous inequality and from (B2) that

$$ \Lambda((u^\theta, n_i^\theta, n_v^\theta), (\tau, \pi)) \geq \theta \Lambda((u^1, n_i^1, n_v^1), (\tau, \pi)) + (1 - \theta) \Lambda((u^0, n_i^0, n_v^0), (\tau, \pi)). $$

860 This shows that $\Lambda$ is concave with respect to its first 3 arguments.

Rewriting $s(u, n_i, n_v) = \inf_{\tau, \pi} \Lambda((u, n_i, n_v); (\tau, \pi))$, one directly gets that $s$ is the infimum of concave functions, and is thus concave (Rockafellar and Wets, 1998, Proposition 2.9(b)).

## B2 ... but is degenerate in our setting

We know from Appendix B1 that the explicit entropy function $s$ defined by Eq. (48) is concave. As a consequence, all the eigenvalues of the Hessian matrix

$$ D^2 s(u, n_i, n_v) = \begin{pmatrix} \partial_u \left(\frac{1}{T}\right) & \partial_{n_i} \left(\frac{1}{T}\right) & \partial_{n_v} \left(\frac{1}{T}\right) \\ \partial_u \left(-\frac{\mu_i}{T}\right) & \partial_{n_i} \left(-\frac{\mu_i}{T}\right) & \partial_{n_v} \left(-\frac{\mu_i}{T}\right) \\ \partial_u \left(-\frac{\mu_v}{T}\right) & \partial_{n_i} \left(-\frac{\mu_v}{T}\right) & \partial_{n_v} \left(-\frac{\mu_v}{T}\right) \end{pmatrix} $$

are non-positive. Noticing that the three quantities $\frac{1}{T}$, $-\frac{\mu_i}{T}$, and $-\frac{\mu_v}{T}$ can be described thanks to the mere two quantities $T$ and $\rho_v = \frac{n_v}{1-n_i \bar{v}_i^0}$ shows that the mapping

$$ (u, n_i, n_v) \mapsto \left( \frac{1}{T}, -\frac{\mu_i}{T}, -\frac{\mu_v}{T} \right) $$

870 cannot be locally invertible in the neighborhood of any point $(u, n_i, n_v)$. As a consequence,

$$ \det D^2 s(u, n_i, n_v) = 0 \quad \text{for all } (u, n_i, n_v). $$

In particular, $s$ is nowhere strictly concave.

## Appendix C: The derivative of the characteristic function in the framework of superderivatives

Superderivatives are a generalization of the notion of derivative to concave non-smooth functions (note that the literature usually deals with the notion of subderivative of convex functions, but the translation to concave functions is straightforward; Rockafellar and Wets, 1998, Chapter 8). In this framework, the superdifferential $\partial f(x)$ of a concave function $f: \mathbb{R}^d \to \mathbb{R} \cup \{-\infty\}$ at a given point $x$ at which $f(x)$ is finite corresponds to the set of tangent slopes that fit above the graph of the function, i.e.

$$ \partial f(x) = \{p \in \mathbb{R}^d \text{ s.t. } f(y) \leq f(x) + p \cdot (y - x) \text{ for all } y \in \mathbb{R}^d\}. $$
880 For points $x$ where $f$ is differentiable, the set $\partial f(x)$ reduces to the singleton $\{\nabla f(x)\}$. The elements of $\partial f(x)$ are referred to as superderivatives of $f$ at point $x$.

A simple example of superderivatives is given for the function $f(x) = -|x|$ in Fig. C1. At the corner point $x = 0$, every slope in the interval $[-1, 1]$ fits above the curve, so that $\partial f(0) = [-1, 1]$. Hence, for $f(x) = |x|$ depicted in Fig. C1, one ends up with

$$ \partial f(x) = \begin{cases} \{1\} & \text{if } x < 0, \\ [-1, 1] & \text{if } x = 0, \\ \{-1\} & \text{if } x < 0. \end{cases} $$

885 The same approach can be used to define the superdifferential of the characteristic function

$$ \chi_{x_0} : x \mapsto \chi_{x_0}(x) = \begin{cases} 0 & \text{if } x = x_0, \\ -\infty & \text{otherwise,} \end{cases} $$

used in Sect. 5 to enforce the incompressibility constraint. As illustrated in Fig. C2 shows that $\chi_{x_0}$ is indeed a concave function and that its superdifferential is given by

$$ \partial \chi_{x_0}(x) = \begin{cases} \mathbb{R} & \text{if } x = x_0, \\ \emptyset & \text{otherwise.} \end{cases} $$

<figure><img src="image_7.png" /><figcaption>**Figure C1.** Example of superderivatives of the convex function $y = -|x|$. The function is non-smooth in the point $x = 0$, and every slope in the interval $[-1, 1]$ fits below the curve. In $x = -2$, the function is smooth and the superderivative reduces to the singleton $\{1\}$, consistently with the standard definition of the derivative.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure C2.** Visualization of the superderivative of the characteristic function $\chi_{x_0}$ in $x_0$.</figcaption></figure>

*Author contributions.* The research was designed by CC and KF. The research was performed by CC, KF, and KJ with help from MD. The manuscript was written by KF with help from all co-authors.

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

*Acknowledgements.* This work was financially supported by the S-NOW project, funded by the Institut des Mathématiques pour la Planète Terre. Kévin Fourteau and Marie Dumont received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (IVORI; grant no. 949516). Clément Cancès also acknowledges support from the French National Research Agency through Labex CEMPI (ANR-11-LABX-0007-01) and project MATHSOUT of the PEPR Mathematics in Interaction (ANR-23-EXMA-0010). The authors further thank Noel J. Walkington for stimulating discussions on the thermodynamically consistent numerical modeling of complex porous media flows.
Bartelt, P. and Buser, O.: The principle of minimum entropy production and snow structure, Journal of Glaciology, 50, 342–352,
https://doi.org/10.3189/172756504781829945, 2004.

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

905 Bartelt, P., Buser, O., and Sokratov, S. A.: A nonequilibrium treatment of heat and mass transfer in alpine snowcovers, Cold Regions Science and Technology, 39, 219–242, https://doi.org/10.1016/j.coldregions.2004.04.005, snow And Avalanches: Papers Presented At The European Geophysical Union Conference, Nice, April 2003. Dedicated To The Avalanche Dynamics Pioneer Dr. B. Salm., 2004.

Bartelt, P., Buser, O., and Kern, M.: Dissipated work, stability and the internal flow structure of granular snow avalanches, Journal of Glaciology, 51, 125–138, https://doi.org/10.3189/172756505781829638, 2005.

910 Bartelt, P., Buser, O., and Platzer, K.: Fluctuation-dissipation relations for granular snow avalanches, Journal of Glaciology, 52, 631–643,
https://doi.org/10.3189/172756506781828476, 2006.

Bessemoulin-Chatard, M. and Chainais-Hillairet, C.: Exponential decay of a finite volume scheme to the thermal equilibrium for drift-diffusion systems, J. Numer. Math., 25, 147–168, https://hal.archives-ouvertes.fr/hal-01250709, 2017.

Bodineau, T., Lebowitz, J., Mouhot, C., and Villani, C.: Lyapunov functionals for boundary-driven nonlinear drift-diffusion equations, Nonlinearity, 27, 2111, https://doi.org/10.1088/0951-7715/27/9/2111, 2014.

Boone, A., Samuelsson, P., Gollvik, S., Napoly, A., Jarlan, L., Brun, E., and Decharme, B.: The interactions between soil-biosphere-atmosphere land surface model with a multi-energy balance (ISBA-MEB) option in SURFEXv8 - Part 1: Model description, Geoscientific Model Development, 10, 843–872, https://doi.org/10.5194/gmd-10-843-2017, 2017.

Born, A., Imhof, M. A., and Stocker, T. F.: An efficient surface energy-mass balance model for snow and ice, The Cryosphere, 13, 1529–1546,
https://doi.org/10.5194/tc-13-1529-2019, 2019.

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

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

Bruining, J., Marchesin, D., and Van Duijn, C. J.: Steam injection into water-saturated porous rock, Comput. Appl. Math., 22, 359–395, 2003.

Butcher, J. C.: Numerical methods for ordinary differential equations, John Wiley & Sons, Ltd, The Atrium, Southern Gate, Chichester, West Sussex PO19 8SQ, England, 2008.

930 Calonne, N., Flin, F., Morin, S., Lesaffre, B., du Roscoat, S. R., and Geindreau, C.: Numerical and experimental investigations of the effective thermal conductivity of snow, Geophysical Research Letters, 38, L23 501, 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, https://doi.org/10.1021/jp5052535, 2014.

Cancès, C.: Energy stable numerical methods for porous media flow type problems, Oil & Gas Science and Technology-Rev. IFPEN, 73, 1–18, 2018.
Cancès, C. and Guichard, C.: Numerical Analysis of a Robust Free Energy Diminishing Finite Volume Scheme for Parabolic Equations with Gradient Structure, Found. Comput. Math., 17, 1525–1584, https://doi.org/10.1007/s10208-016-9328-6, 2017.

Cancès, C., Monsaingeon, L., and Natale, A.: Discretizing the Fokker-Planck equation with second-order accuracy: a dissipation driven approach, hAL: hal-04720129, 2024.

940
Cancès, C. and Gaudeul, B.: A Convergent Entropy Diminishing Finite Volume Scheme for a Cross-Diffusion System, SIAM Journal on Numerical Analysis, 58, 2684–2710, https://doi.org/10.1137/20M1316093, 2020.

Cancès, C. and Guichard, C.: Convergence of a nonlinear entropy diminishing control volume finite element scheme for solving anisotropic degenerate parabolic equations, Mathematics of Computation, 85, 549–580, https://doi.org/10.1090/mcom/2997, 2016.

Carrillo, J. A., Jüngel, A., Markowich, P. A., Toscani, G., and Unterreiter, A.: Entropy dissipation methods for degenerate parabolic problems and generalized Sobolev inequalities, Monatsh. Math., 133, 1–82, https://doi.org/10.1007/s006050170032, 2001.

945

Chainais-Hillairet, C.: Entropy method and asymptotic behaviours of finite volume schemes, in: Finite Volumes for Complex Applications VII-Methods and Theoretical Aspects: FVCA 7, Berlin, June 2014, pp. 17–35, Springer, https://doi.org/10.1007/978-3-319-05684-5_2, 2014.

950
Chainais-Hillairet, C. and Herda, M.: Large-time behavior of a family of finite volume schemes for boundary-driven convection-diffusion equations, IMA J. Numer. Anal., 40, 2473–2504, https://arxiv.org/abs/1810.01087, 2020.

Chatzipantelidis, P., Horváth, Z., and Thomée, V.: On preservation of positivity in some finite element methods for the heat equation, Comput. Methods Appl. Math., 15, 417–437, https://doi.org/10.1515/cmam-2015-0018, 2015.

955
Ciarlet, P. and Raviart, P.-A.: Maximum principle and uniform convergence for the finite element method, Computer Methods in Applied Mechanics and Engineering, 2, 17–31, https://doi.org/10.1016/0045-7825(73)90019-4, 1973.

955

Clarke, F. H.: Optimization and nonsmooth analysis, vol. 5 of *Classics in Applied Mathematics*, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edn., https://doi.org/10.1137/1.9781611971309, 1990.

Coussy, O.: Poromechanics, John Wiley & Sons, Ltd, https://doi.org/10.1002/0470092718, 2004.

Dafermos, C. M.: Hyperbolic Conservation Laws in Continuum Physics, Springer, 2010.

Decharme, B., Brun, E., Boone, A., Delire, C., Le Moigne, P., and Morin, S.: Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model, The Cryosphere, 10, 853–877, https://doi.org/10.5194/tc-10-853-2016, 2016.

960

Di Vita, A.: On the Thermodynamics of a Room, a Heater and a Window, Physical Science International Journal, 4, 166–180, https://journalpsij.com/index.php/PSIJ/article/view/139, last accessed: 09/07/2025, 2013.

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, Journal of Advances in Modeling Earth Systems, 11, 34–44, https://doi.org/10.1029/2018MS001445, 2019.

965

Droniou, J.: Finite volume schemes for diffusion equations: introduction to and review of modern methods, Math. Models Methods Appl. Sci., 24, 1575–1620, 2014.

Essery, R.: A factorial snowpack model (FSM 1.0), Geoscientific Model Development, 8, 3867–3876, https://doi.org/10.5194/gmd-8-3867-2015, 2015.

970

Essery, R., Mazzotti, G., Barr, S., Jonas, T., Quaife, T., and Rutter, N.: A Flexible Snow Model (FSM 2.1.1) including a forest canopy, Geoscientific Model Development, 18, 3583–3605, https://doi.org/10.5194/gmd-18-3583-2025, 2025.
Estilow, T. W., Young, A. H., and Robinson, D. A.: A long-term Northern Hemisphere snow cover extent data record for climate studies and monitoring, Earth Sys. Sci. Data, 7, 137–142, https://doi.org/10.5194/essd-7-137-2015, 2015.

975 Eymard, R. and Tillier, E.: Mathematical and numerical study of a system of conservation laws, J. Evol. Equ., 7, 197–239, 2007.

Eymard, R., Gallouët, T., and Herbin, R.: Finite volume methods, Ciarlet, P. G. (ed.) et al., in Handbook of numerical analysis. North-Holland, Amsterdam, pp. 713–1020, 2000.

Faccanioni, G., Galusinski, C., and Grec, B.: A conservative two-phase flow model with a nonlinear degenerate diffusion, https://hal.science/hal-05021438, hAL: hal-05021438, 2025.

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

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

Glansdorff, P. and Prigogine, I.: Non-equilibrium stability theory, Physica, 46, 344–366, https://doi.org/10.1016/0031-8914(70)90010-8, 1970.

Godlewski, E. and Raviart, P.-A.: Numerical approximation of hyperbolic systems of conservation laws, vol. 118 of *Applied Mathematical Sciences*, Springer-Verlag, New York, 1996.

Godunov, S. K.: Finite Difference Method for Numerical Computation of Discontinuous Solutions of the Equations of Fluid Dynamics, Mat. Sb., 47 (89), 271, 1959.

990 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.

Hoang, T.-K.-D., Quiquet, A., Dumas, C., Born, A., and Roche, D. M.: Using a multi-layer snow model for transient paleo-studies: surface mass balance evolution during the Last Interglacial, Climate of the Past, 21, 27–51, https://doi.org/10.5194/cp-21-27-2025, 2025.

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

Jüngel, A.: The boundedness-by-entropy method for cross-diffusion systems, Nonlinearity, 28, 1963–2001, 2015.

Jüngel, A.: Entropy Methods for Diffusive Partial Differential Equations, Springer Nature, https://doi.org/10.1007/978-3-319-34219-1, 2016.

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

Kardar, M.: Statistical Physics of Particles, Cambridge University Press, https://doi.org/10.1017/CBO9780511815898, 2007.

Keizer, J. and Fox, R. F.: Qualms Regarding the Range of Validity of the Glansdorff-Prigogine Criterion for Stability of Non-Equilibrium States, Proceedings of the National Academy of Sciences, 71, 192–196, https://doi.org/10.1073/pnas.71.1.192, 1974.

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

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, Geoscientific Model Development, 11, 5027–5049, https://doi.org/10.5194/gmd-11-5027-2018, 2018.

Lebon, G., Jou, D., and Casas-Vázquez, J.: Understanding Non-equilibrium Thermodynamics, Springer-Verlag, https://doi.org/10.1007/978-3-540-74252-4, 2008.

1020

Lenaerts, J. T. M., Medley, B., van den Broeke, M. R., and Wouters, B.: Observing and Modeling Ice Sheet Surface Mass Balance, Rev. Geophys., 57, 376–420, https://doi.org/10.1029/2018RG000622, 2019.

Lide, D. R.: CRC handbook of chemistry and physics, 86 edn., Taylor and Francis, Boca Raton, FL, 2005.

Liu, C., Wang, C., and Wang, Y.: A structure-preserving, operator splitting scheme for reaction-diffusion equations with detailed balance, Journal of Computational Physics, 436, 110 253, https://doi.org/10.1016/j.jcp.2021.110253, 2021.

1025

Lyapunov, A.: The general problem of the stability of motion, Int. J. Control, 55, 531–534, https://doi.org/10.1080/00207179208934253, 1992.

Maes, C. and Netočný, K.: Revisiting the Glansdorff–Prigogine Criterion for Stability Within Irreversible Thermodynamics, Journal of Statistical Physics, 159, 1286–1299, https://doi.org/10.1007/s10955-015-1239-4, 2015.

1030

Marty, C., Schlögl, S., Bavay, M., and Lehning, M.: How much can we save? Impact of different emission scenarios on future snow cover in the Alps, The Cryosphere, 11, 517–529, https://doi.org/10.5194/tc-11-517-2017, 2017.

Mortimer, R. G. and Eyring, H.: Elementary transition state theory of the Soret and Dufour effects, Proc. Nat. Acad. Sci., 77-4, https://doi.org/10.1073/pnas.77.4.1728, 1980.

Mott, R., Winstral, A., Cluzet, B., Helbig, N., Magnusson, J., Mazzotti, G., Quéno, L., Schirmer, M., Webster, C., and Jonas, T.: Operational snow-hydrological modeling for Switzerland, Front. Earth Sci., 11, https://doi.org/10.3389/feart.2023.1228158, 2023.

1035

Moyne, C., Batsale, J.-C., and Degiovanni, A.: Approche expérimentale et théorique de la conductivité thermique des milieux poreux humides—II. Théorie, International Journal of Heat and Mass Transfer, 31, 2319–2330, https://doi.org/10.1016/0017-9310(88)90163-9, 1988.

Rank, E., Katz, C., and Werner, H.: On the importance of the discrete maximum principle in transient analysis using finite element methods, Int. J. Numerical Methods Engineering, 19, 1771–1782, https://doi.org/10.1002/nme.1620191205, 1983.

1040

Rockafellar, R. and Wets, R. J.-B.: Variational Analysis, Springer Berlin Heidelberg, Berlin, Heidelberg, https://doi.org/10.1007/978-3-642-02431-3_8, 1998.

Royer, A., Picard, G., Vargel, C., Langlois, A., Gouttevin, I., and Dumont, M.: Improved Simulation of Arctic Circumpolar Land Area Snow Properties and Soil Temperatures, Front. Earth Sci., 9, https://doi.org/10.3389/feart.2021.685140, 2021.

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

1045

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, 2022.
Sturm, M. and Massom, R. A.: Snow in the sea ice system: friend or foe?, chap. 3, pp. 65–109, John Wiley & Sons, Ltd, https://doi.org/10.1002/9781118778371.ch3, 2017.

1050 Thackeray, C. W., Derksen, C., Fletcher, C. G., and Hall, A.: Snow and climate: Feedbacks, drivers, and indices of change, Current Clim. Change Reports, 5, 322–333, https://doi.org/10.1007/s40641-019-00143-w, 2019.

Vavrus, S.: The role of terrestrial snow cover in the climate system, Clim. Dyn., 29, 73–88, https://doi.org/10.1007/s00382-007-0226-0, 2007.

Vidal, C., Dewel, G., and Borckmans, P.: Au-dela de l'équilibre, Hermann, 1994.

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

Webster, M., Gerland, S., Holland, M., Hunke, E., Kwok, R., Lecomte, O., Massom, R., Perovich, D., and Sturm, M.: Snow in the changing sea-ice systems, Nature Clim. Change, 8, 946–953, https://doi.org/10.1038/s41558-018-0286-7, 2018.

Zhang, F.: The Schur Complement and Its Applications, Springer New York, https://doi.org/10.1007/b105056, 2006.