Some insights from the second principle of thermodynamics for snowpack modeling

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.