the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 12 Feb 2026
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.
Competing interests: At least one of the (co-)authors is a member of the editorial board of The Cryosphere.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.- Preprint
(1016 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
- CC1: 'Comment on egusphere-2026-510', Pascal Marquet, 23 Feb 2026
-
RC1: 'Comment on egusphere-2026-510', Anonymous Referee #1, 25 Mar 2026
This research aims to improve numerical stability by incorporating constraints derived from the second law of thermodynamics into the model and by developing an algorithm that prevents violations of physical laws, such as decreases in entropy during numerical calculations. Numerical snowpack models involve complex processes, including phenomena that require fine spatial resolution and others that evolve over short timescales, making it difficult to adopt longer time steps. In this respect, the present study has the potential to improve numerical snowpack models in terms of numerical stability while also reducing computational time.
My suggestion is to include a discussion of computational cost. For example, when the same time step is used, does the use of the Tightly Coupled scheme and the Backword Euler method increase the computational cost? If so, to what extent. This discussion helps to explain the advantage to incorporate these schemes in terms of computational cost.
This study addresses several important processes, including coupling with a canopy model and water vapor transport. I expect that, in the future, the framework can be extended to include liquid water movement as well and be further integrated into a numerical snowpack model.
Minor comments
L134-135 The oscillation shown in Fig. 1 is problematic for numerical calculations. In contrast, the entropy-based method appears to provide greater numerical stability and relax the time-step constraint, which may contribute to improving numerical snowpack models. It would also be valuable to include a discussion somewhere in the manuscript on how much this method could reduce computational cost.
L388-398 Figure 3 compares sequential and tightly-coupled simulations to demonstrate whether oscillations occur. Could the authors provide a version of the figure with the same y-axis scale? If they-axis scales for temperature and entropy differ between the upper and lower figures, it becomes difficult to judge whether there is a difference between the two cases in inital stage. For example, has entropy reduction and temperature oscillation already started in the initial stage?
L398- 401 I understand that the thermal inertia of canopy air is very small. However, if the calculation remains unstable no matter how small the time step is because the thermal inertia is set to zero, this may suggest an issue with the calculation setup. I don't think it is necessary to repeat the calculations with non-zero thermal interia included, but it would be better if the authors could comment on the magnitude of the thermal inertia if it is assigned to the canopy air.
L401-403 The fact that the calculations remain stable with very long time steps suggests that it may be possible to calculate changes in external conditions, such as input meteorological data updated hourly, using the same time step. I also wonder whether rapidly varying processes, such as liquid water movement, might likewise be handled with such long time steps. This point may be better discussed in the Discussion section.
L497 Is “uv0 and uv0” in Eq. (37) a typo? Should it be “uv0 and sv0”?
L702-703 Are there any actual observations of snow density oscillations on crusts that may be due to similar mechanisms, including findings from previous studies?
L711-727 Since this study aims to achieve stable calculations, it is important to summarize how instability can be avoided when a coarse grid provide a new source of instability. The manuscript states that the solution is to use downstream values, but the reason for this solution is unclear. My intuition as a non-specialist is that using downstream values may reduce the amount of heat and water vapor transported, thereby suppressing the oscillation. I therefore wonder whether this is only one specific case in which the using downstream value avoided the oscillation, or whether the using downstream value consistently prevents oscillation across multiple calculation patterns. If it is the former, its usefulness as a general solution may be limited. This point should be discussed more, as it is fundamental to the study’s objective of achieving stable computation.
L730-736 Liquid water infiltration is a major cause of rapid behavior and increased computational cost. If this method can be used to allow larger time steps for such processes, I would expect it to lead to substantial improvement in numerical snowpack models.
Citation: https://doi.org/10.5194/egusphere-2026-510-RC1 -
RC2: 'Comment on egusphere-2026-510', Anonymous Referee #2, 04 Apr 2026
The work presented in this article is substantial and ambitious. The proposed entropy-based framework for snowpack modelling, together with the coupling strategy and the number of governing equations involved, reflects a considerable theoretical and numerical effort. In particular, the attempt to enforce entropy consistency within the thermal exchange equations of a snowpack system is technically demanding and represents a noteworthy contribution. More broadly, the manuscript addresses an important and original topic, namely the use of the second law of thermodynamics as a criterion for model consistency and numerical robustness.
At the same time, in its current form, the manuscript reads primarily as a methodological and conceptual framework, with a strong didactic component, rather than as a modelling development whose practical advantages are already fully demonstrated. While the thermodynamic analysis is interesting and carefully developed, the manuscript would benefit from a clearer demonstration of the concrete advantages associated with the introduction of this additional level of complexity. In particular, it would be helpful to better highlight the gains obtained in terms of numerical stability, predictive capability, or computational performance relative to more conventional snowpack formulations based primarily on the first principle.
In my view, the work has the potential to be further strengthened by extending the analysis beyond thermodynamic consistency alone. For example, an exergy-based assessment could offer complementary physical insight by helping to identify the snowpack layers where energy quality degradation and irreversible losses are most significant. More generally, the manuscript would benefit from a clearer quantitative demonstration of the practical added value of the proposed framework, especially in comparison with a conventional snowpack model. At present, the motivation for introducing this added complexity would be more compelling if the manuscript more explicitly showed whether the entropy-based formulation leads to a measurable improvement over a standard first-law-based approach.
Major comments
1. Possible added value of exergy analysis:
The manuscript would benefit from a brief discussion of whether an exergy-based framework could complement the present entropy analysis. Whereas entropy production is well suited to assessing thermodynamic consistency and identifying dissipation mechanisms, exergy could provide an additional and highly relevant physical insight by quantifying where the quality of energy is degraded, where the largest losses of useful energy occur, and which regions of the snowpack are the most dissipative.
This may be particularly valuable in zones characterized by strong thermal gradients. For instance, by associating each heat transfer term with its corresponding Carnot efficiency, it would be possible to identify the regions where exergy destruction is maximal. A similar approach is commonly adopted in thermodynamic cycle analysis, such as in the Rankine cycle, where the boiler is typically identified as the most critical component due to the large irreversibilities occurring at high temperature. By contrast, although the condenser often accounts for a larger energy dissipation from a first-law perspective, the quality of that energy is significantly lower because it is released at low temperature, and therefore the associated exergy remains relatively small. Applying the same concept in the present context could provide a more insightful assessment of the snowpack by identifying the most critical layers, namely those characterized by the highest-quality thermal energy and the largest thermodynamic losses. This could in turn guide the adoption of locally refined meshes or other targeted computational strategies, ultimately leading to a model that can more effectively exploit the second principle to achieve significant improvements in both physical insight and numerical performance.
Even if this extension lies beyond the scope of the present study, a brief discussion of this possible direction would significantly strengthen the manuscript and provide a valuable perspective for future model development.2. Entropy source in Figure 3:
I have some concerns regarding the interpretation of the entropy-related quantity shown in Figure 3. If this curve represents the internal entropy production (or entropy source) associated with the snowpack–canopy energy exchange, then it should remain non-negative according to the second principle of thermodynamics. In that case, negative values would indeed indicate either a thermodynamic inconsistency, a sign-convention issue, or a possible problem in the formulation. By contrast, if the plotted quantity instead represents the time variation of the system entropy or another term of the entropy balance, then negative values may be physically admissible, particularly in the case of an open system. I therefore suggest that the authors clearly clarify what quantity is actually plotted in Figure 3 and explicitly explain how it relates to the entropy balance introduced earlier in the manuscript.3. Need for a baseline comparison model
A quantitative comparison with a baseline first-law-only snowpack model would be particularly valuable. In addition to the general discussion introduced above, I encourage the authors to provide explicit benchmark results comparing the two approaches in terms of timestep sensitivity, computational cost, convergence behavior, and the resulting thermodynamic fields. This would allow the reader to more clearly assess whether the proposed second-law-based formulation provides a measurable practical advantage.4. Row 705:
Around line 705, the authors state that the use of a second principle–complying model allows simulations to be performed with very large time steps without numerical instabilities arising from the faster processes. I believe this statement requires a stronger justification and a more detailed discussion. In particular, I would like to better understand the numerical and physical reasons behind this claimed improvement in stability. Based on the results presented, especially in the first example, some instability or at least questionable behavior appears to be visible, which makes this conclusion less immediately convincing. For this reason, I encourage the authors to provide additional evidence, such as a comparison with a classical first-principle snowpack model, timestep sensitivity analyses, or a clearer explanation of how the second-law formulation specifically suppresses the instabilities associated with the faster processes.5. Figures 4 and 5
I would encourage the authors to also display the entropy generation field in Figures 4 and 5. At present, only snow density–related results are presented. Displaying the entropy generation would significantly strengthen the argument of the paper, especially considering that the proposed model is specifically designed to account for the second principle of thermodynamics. If negative values of the entropy source also appear in these cases, this may suggest potential issues in the coupled formulation or in the numerical implementation More generally, I also recommend including a baseline reference case based on a conventional first-law snowpack formulation, in order to provide a clearer comparison with the proposed approach.6. Third example missing (around line 730)
The manuscript refers to three different cases; however, only two numerical examples are explicitly presented. I would strongly encourage including the third case as well, particularly the simplest conduction-only example, as this would provide a useful validation benchmark. From my point of view, it is important to first present simpler and more accessible examples before moving to the more complex and lengthy formulations involving multiple coupled terms and governing equations. This would help the reader better understand the physical meaning of the proposed entropy-based approach and more clearly assess its added value.7. Justification of complexity
At present, the manuscript introduces a large number of additional equations and thermodynamic terms. However, the numerical examples do not yet clearly demonstrate a substantial improvement over more classical modelling approaches. For this reason, I encourage the authors to better justify why this additional level of complexity is necessary and to clearly highlight the measurable advantages it provides, for instance in terms of accuracy, physical consistency, or numerical robustness.
I also think the discussion of numerical stability could be strengthened. While the manuscript suggests that second-law-based consistency may improve robustness, tighter coupling and additional source terms can also make the numerical system more difficult to solve in practice if not handled carefully. It would be valuable to discuss whether alternative numerical strategies for improving stability have been considered, such as under-relaxation techniques, iterative damping methods, adaptive timestep control, or other numerical stabilization approaches. From the perspective of numerical stability, it may therefore be worth discussing whether other approaches could provide similar or improved robustness while introducing a lower level of model complexity.The manuscript is scientifically substantial and generally well structured, but in its current form it is rather long and at times difficult to read. The progression from the heat-equation example to the coupled snow–canopy case and finally to the vapor-transport formulation is logical, but the presentation becomes increasingly dense, especially in the later sections where many equations, variables, and thermodynamic relations are introduced in a short space. I therefore encourage the authors to improve readability by shortening or streamlining some derivations, clarifying the main take-home messages at the beginning and end of each section, and carefully revising the language throughout. At present, the paper reads more like an extended theoretical development than a focused journal article, and a tighter presentation would make its main contributions easier to appreciate.
Minor language revisions (some)
• Helmoltz -> Helmholtz
• modelling / modeling: choose one dialect and apply it everywhere.
• thermodynamical / thermodynamic: both exist, but thermodynamic is usually cleaner
• orientated -> oriented in most cases.
• non-linear / nonlinear: choose one and keep it consistent.Citation: https://doi.org/10.5194/egusphere-2026-510-RC2
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 186 | 125 | 21 | 332 | 40 | 52 |
- HTML: 186
- PDF: 125
- XML: 21
- Total: 332
- BibTeX: 40
- EndNote: 52
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Dear Co-Authors
The entropy equation is studied in the paper for the entropy of a snow-pack defined as a mixture of water vapor and ice (Ih), with ice-vapor phase changes.
The entropy of this mixture is assumed to be (p.8, line 217) the sum: S^n = Sum_{k=1}^N (\Delta z_k) s^n_k
The entropy of ice (Ih) depends on a reference value s_i^0, with (p.18, line 446): s_i = s_i^0 + (...)
Similarly the entropy of water vapour depends on a reference value s_v^0, with (p.19, line 479): s_v = s_v^0 + (...)
These reference values are determined in the Section 5.1.3 ("Defining phase equilibrium", p.19-20), with (p.20, lines 491-492): "... we first set the reference temperature to T0=273.15K, without loss of the generality. We then assume (still without loss of generality) that the molar internal energy and entropy of ice vanishes at T0, yielding that u_i^0=0 and s_i^0=0 ..." Accordingly (p.20, lines 499-500): "... s_v^0 = (s_i^0=0) + [ Delta-h(i->v) ] / T0. This naturally enforces the equality of the chemical potentials, and hence phase equilibrium, in the reference conditions."
Indeed, for such a simple system made of two phases of the same substance (water ice and water vapour) the reference entropy can be set to zero for one or another of the phase, with the other computed from the other by adding or removing the quantity [ Delta-h(i->v) ] / T0, where Delta-h(l->v) is the latent heat of sublimation.
However, for a more general mixture made of several substances, we known since Gibbs (1875-1878) that the reference values cannot be determined arbitrarily, and we know since Planck (1914, 1917) that these reference must be determined from the third-law of thermodynamics (with concrete calorimetric and/or statistical values computed by Sackur, Tetrode, Gordon, Barnes, Giauque, ...).
These third-law values should be key features to be determined in next more realistic studies including the Deuterium oxide (D2O-vapour/liq/ice) species, which cannot be linked with the water (H2O) species. These generalisations should be considered after the Section 6 ("Discussion" p.28-20) 6.1 "Extension to other cases"), where (p.28, line 732-736): "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." With the aim (p.28, line 743) to: "Derive the general entropy function of snow by adding the entropy contribution of each phase."
I provide in the following examples of reference entropies for D2O species from recent and old Thermodynamic Tables that I have used to compute and study the absolute entropies for the moist-air atmosphere (since Marquet, QJRMS, 2011) and the sea-salt oceans (since Marquet, 2026-a,b). See my publications in https://sites.google.com/view/pascal-marquet/
Best regards,
Dr.Hab. Pascal Marquet
pascalmarquet@yahoo.com
%=============================================
% "Atkins' Physical chemistry (12th edition)"
% Oxford University Press (2023) 927~Pp.
%---------------------------------------------
% Table 2C.7 (p.892) 25°C (1 atm)
% D2O M(g/mol) S(J/K/mol) Cp(J/K/mol)
% vap 20.028 198.34 34.27
% liq 20.028 75.94 84.35
% D2O M(g/mol) S(J/K/kg) Cp(J/K/kg)
% vap 20.028 9903.1 1711.1
% liq 20.028 3791.7 4211.6
%---------------------------------------------
% Table 2C.7 (p.893) 25°C (1 atm)
% H2O M(g/mol) S(J/K/mol) Cp(J/K/mol)
% vap 18.015 188.83 33.58
% liq 18.015 69.91 75.291
% ice 18.015 37.99
% H2O M(g/mol) S(J/K/kg) Cp(J/K/kg)
% vap 20.028 10482 1864
% liq 20.028 3881 4179.4
% ice 18.015 2109
%=============================================
%=======================================================
% Wagman et al. (1982)
% "The NBS tables of chemical thermodynamic properties"
% J. Phys. Chem. Ref. Data, Vol.11, Supplement Number-2,
% 392-Pp. https://srd.nist.gov/JPCRD/jpcrdS2Vol11.pdf
%-------------------------------------------------------
% Table-2H (p.2-38) 298.15K 25°C and 1bar (0.1MPa)
% D2O M(g/mol) S(J/K/mol) Cp(J/K/mol)
% vap 20.0276 198.339 34.27
% liq 20.0276 75.94 84.35
%-------------------------------------------------------
% Table-2H (p.2-38) 298.15K 25°C and 1bar (0.1MPa)
% H2O M(g/mol) S(J/K/mol) Cp(J/K/mol)
% vap 18.0154 188.825 33.577
% liq 18.0154 69.91 75.291
%=======================================================
%=============================================
% Long and Kemp (1936)
% "The Entropy of Deuterium Oxide and
% the Third Law of Thermodynamics.
% Heat Capacity of Deuterium Oxide
% from 15 to 298°K.
% The Melting Point and Heat of Fusion"
%---------------------------------------------
% Journal of the American Chemical Society
% Vol-58, Number-10, October 9, 1936
%---------------------------------------------
% (1 atm) S0
% D2O gas at 273.10K: 45.89 pm 0.1 cal/K/mol
% D2O liq at 298.10K: 17.27 pm 0.05 cal/K/mol
% D2O ice at 276.92K: 10.38 cal/K/mol
%---------------------------------------------
% (1 atm) S0
% D2O gas at 273.10K: 192.0 pm 0.4 J/K/mol
% D2O liq at 298.10K: 72.3 pm 0.2 J/K/mol
% D2O ice at 276.92K: 43.43 J/K/mol
%=============================================
Heavy Water thermodynamic:
https://inis.iaea.org/records/bgt3d-ng660
Hill_MacMillan_Vee_Heavy_Water_1981.pdf (HMMV-81)
%=============================================
https://webbook.nist.gov/cgi/cbook.cgi?ID=C7789200&Mask=6F
Deuterium oxide (NIST page): D2O
Gas (1 bar) : Svap(1000hPa)=198.34 J/K/mol (Chase 1998)
---------------------------------------------
D2O / Chase (1998, Ideal gas, p.1045):
Mr = 20.020004-g/mol = 0.020020004-kg/mol
R/Mr = 8.3145/0.020020004 = 415.3 J/K/kg
---------------------------------------------
Chase (1998) D2O Gas (1 bar 25°C) : Svap(1000hPa)=198.339 J/K/mol
Chase (1998) D2O Gas (1 bar 25°C) : Svap(1000hPa)=9907.04 J/K/kg
---------------------------------------------
HMMV-81 / D2O / 3.82°C / Table-1, p.1 / T-fus = 276.97-K
(HMMV-81, 3.82°C) P-sat = 0.6601-kPa = 660.1 Pa
(HMMV-81, 3.82°C) Delta(S)-(liq/vap) = 8390.2-J/K/kg
(HMMV-81, 3.82°C) Delta(H)-(liq/vap) = 2323.7-kJ/kg
---------------------------------------------
HMMV-81 / D2O / 25°C / Table-1, p.1 :
(HMMV-81, 25°C) P-sat = 2.737-kPa = 2737 Pa = 27.37 hPa
(HMMV-81, 25°C) Delta(S)-(liq/vap) = 7611.6-J/K/kg = Sv-Sl
(HMMV-81, 25°C) Delta(H)-(liq/vap) = 2269.4-kJ/kg = L-vap
(HMMV-81, 25°C) Delta(S)-(liq/vap) = 152.38-J/K/mol = Sv-Sl
(HMMV-81, 25°C) Delta(H)-(liq/vap) = 45433.4-J/mol = L-vap
- - - - - - - - - - - - - - - - - - - - - - - -
0) Delta(H)-(liq/vap) / 298.15 = Delta(S)-(liq/vap)
45433.4 / 298.15 = 152.384 : OK
- - - - - - - - - - - - - - - - - - - - - - - -
1) S-vap(27.37 hPa) = S-vap(1000hPa) - R*ln(27.37/1000)
S-vap(27.37 hPa) = 198.339 - 8.314*ln(27.37/1000)
S-vap(27.37 hPa) = 228.26 J/K/mol
- - - - - - - - - - - - - - - - - - - - - - - -
2) S-liq = S-vap(27.37 hPa) - Delta(S)-(l/v)
= 228.26 - 152.38
S-liq (25°C) = 75.88 J/K/mol
- - - - - - - - - - - - - - - - - - - - - - - -
3) D2O(vap) 100 200 298.15 300
Cp (J/K/kg) 33.299 33.449 34.255 34.278
=============================================================
=============================================================
(normal) Water (NIST page): H2O
- - - - - - - - - - - - - - - - - - - - - - - -
https://webbook.nist.gov/cgi/cbook.cgi?ID=C7732185&Mask=1E9F
- - - - - - - - - - - - - - - - - - - - - - - -
Gas (1 bar 25°C) : S0=188.885 J/K/mol (Cox, Wagmann et al, 1984)
Gas (1 bar 25°C) : S0=188.84 J/K/mol (Chase 1998)
Gas (1 bar 25°C) : S0=10482.2 J/K/kg (Chase 1998)
---------------------------------------------
H2O / Chase (1998, Ideal gas 25°C, p.1324):
Mr=18.01528 g/mol
R/Mr = 8.3145/0.01801528 = 461.5 J/K/kg
H2O / Gas (1 bar) : S0=188.834 +/-0.042 J/K/mol
H2O / Liquid (1 bar) : S0=69.950 +/-0.079 J/K/mol
- - - - - - - - - - - - - - - - - - - - - - -
25°C P-sat = 3169 Pa = 31.69 hPa
L-vap = 43988 J/mol = 2441.7 kJ/kg
S-vap = 69.950 + 43988/298.15 + 8.314*ln(31.69/1000)
188.79 = 69.950 + 147.536 - 28.70
close to 188.84 indeed
=============================================================