Research article 07 Dec 2021
Research article  07 Dec 2021
Elements of future snowpack modeling – Part 2: A modular and extendable Eulerian–Lagrangian numerical scheme for coupled transport, phase changes and settling processes
 ^{1}AICES Graduate School, RWTH Aachen University, Schinkelstr. 2a, 52062 Aachen, Germany
 ^{2}WSL Institute for Snow and Avalanche Research SLF, Flüelastr. 11, 7260 Davos, Switzerland
 ^{3}Computational Geoscience, University of Göttingen, Goldschmidtstr. 1, 37077 Göttingen, Germany
 ^{4}Methods for Modelbased Development in Computational Engineering, RWTH Aachen University, Eilfschornsteinstraße 18, 52062 Aachen, Germany
 ^{1}AICES Graduate School, RWTH Aachen University, Schinkelstr. 2a, 52062 Aachen, Germany
 ^{2}WSL Institute for Snow and Avalanche Research SLF, Flüelastr. 11, 7260 Davos, Switzerland
 ^{3}Computational Geoscience, University of Göttingen, Goldschmidtstr. 1, 37077 Göttingen, Germany
 ^{4}Methods for Modelbased Development in Computational Engineering, RWTH Aachen University, Eilfschornsteinstraße 18, 52062 Aachen, Germany
Correspondence: Anna Simson (simson@aices.rwthaachen.de)
Hide author detailsCorrespondence: Anna Simson (simson@aices.rwthaachen.de)
A coupled treatment of transport processes, phase changes and mechanical settling is the core of any detailed snowpack model. A key concept underlying the majority of these models is the notion of layers as deforming material elements that carry the information on their physical state. Thereby an explicit numerical solution of the ice mass continuity equation can be circumvented, although with the downside of virtual no flexibility in implementing different coupling schemes for densification, phase changes and transport. As a remedy we consistently recast the numerical core of a snowpack model into an extendable Eulerian–Lagrangian framework for solving the coupled nonlinear processes. In the proposed scheme, we explicitly solve the most general form of the ice mass balance using the method of characteristics, a Lagrangian method. The underlying coordinate transformation is employed to state a finitedifference formulation for the superimposed (vapor and heat) transport equations which are treated in their Eulerian form on a moving, spatially nonuniform grid that includes the snow surface as a free upper boundary. This formulation allows us to unify the different existing viewpoints of densification in snow or firn models in a flexible way and yields a stable coupling of the advectiondominated mechanical settling with the remaining equations. The flexibility of the scheme is demonstrated within several numerical experiments using a modular solver strategy. We focus on emerging heterogeneities in (twolayer) snowpacks, the coupling of (solid–vapor) phase changes with settling at layer interfaces and the impact of switching to a nonlinear mechanical constitutive law. Lastly, we discuss the potential of the scheme for extensions like a dynamical equation for the surface mass balance or the coupling to liquid water flow.
The snow density is probably the most important prognostic variable of any snowpack model as, e.g., reflected by a focus on snow water equivalent in past snow model intercomparison projects (Krinner et al., 2018, and references therein). That said, it actually comes as a surprise when not even the detailed snowpack models, e.g., Crocus and SNOWPACK, explicitly state an ice mass conservation equation in their technical documentation (Brun et al., 1989, 1992; Lehning et al., 2002; Bartelt and Lehning, 2002). Only a more detailed inspection reveals how mass conservation is accounted for, namely rather indirectly by stating a settling law for individual layers and resorting to a “Lagrangian coordinate system that moves with the ice matrix” (Bartelt and Lehning, 2002) to translate the icephase deformation into a thickness evolution of the layers (Brun et al., 1989; Vionnet et al., 2012). While this procedure has been well established for a long time, it is without numerical ambiguities only in the absence of phase changes. In addition, this nonexplicit nature of the most important conservation law in snow makes it virtually impossible to isolate and advance the numerical core of a snowpack model as an encapsulated numerical scheme comprising all involved coupled nonlinear partial differential equations.
This nonexplicit treatment of snow density or ice mass continuity in snowpack models, e.g., SNOWPACK and Crocus, has to be contrasted to other existing work on densification, comprising both standalone numerical snow studies (Meyer et al., 2020) and the vast body of work on firn densification (Lundin et al., 2017). All of the latter models are built around an explicit formulation of the ice mass continuity equation. This conceptual difference renders a general comparison of firn and snow densification mechanisms (Lundin et al., 2017) difficult. For model intercomparisons in the future it is thus desirable to have a numerical core that is able to digest arbitrary snow or firn densification physics with a flexible but rigorous coupling to superimposed nonlinear transport and phase change processes.
Any holistic snowpack model has to account for transport of heat, vapor and liquid water and its induced phase change processes, as well as mechanical settling and apparent metamorphic processes acting on the snow's microstructure. A widespread body of literature exists that proposes different modeling approaches and computational tools for the various flavors and perspectives of this multiphysics coupled situation, e.g., Krinner et al. (2018, and references therein). For the general timescales of interest (diurnal up to seasonal), it is common practice to employ a continuum assumption and to model the snowpack's state as a mixture of ice, vapor, water and air, as initially described in Bader and Weilenmann (1992). Detailed snowpack models, such as SNOWPACK (Lehning et al., 2002) and Crocus (Vionnet et al., 2012), are built upon this type of mixture theory approach and used for a wide range of purposes.
While heat transport, mechanical settling and processes due to the presence of liquid water have been incorporated into SNOWPACK (Bartelt and Lehning, 2002) and Crocus (Vionnet et al., 2012) for a long time, effects due to vapor transport have mostly recently been investigated in separate studies. Temperature gradients between the ground and atmosphere imply upward vapor fluxes in snowpacks. Stronger temperature gradients (due to either a smaller snowpack height or colder surface temperature) in arctic conditions yield higher vapor fluxes (Domine et al., 2019) compared with alpine snowpacks. Depth hoar layers with reduced density and thermal conductivity form at the snowpack's bottom. In alpine snowpacks similar hoar layers develop within the snowpack that may cause avalanches due to their low mechanical stability (Schweizer et al., 2003). Upscaled and homogenized continuum mechanical process models that account for vapor transport have been put forward, for instance by Hansen and Foslien (2015) and Calonne et al. (2014). Both couple the snowpack's evolving temperature profiles to a nonlinear reaction–diffusion type of equation for vapor transport and phase change. While they provide different flavors of how to set up the underlying mathematical model, both approaches are formulated for idealized conditions and investigate vapor diffusion in the absence of settling and therefore neglect its feedback on the apparent snow density. These modelbased investigations and also fieldbased observations in arctic snowpacks on top of permafrost (Domine et al., 2016, 2019) have demonstrated the significance of vaporrelated processes in snow. Hence, it is of great interest to investigate further how vapor interacts with apparent mechanical processes within the snowpack.
Incorporating vapor transport directly into a fully coupled snowpack model is however challenging, e.g., due to the fact that the associated characteristic timescales are small, and expected effects on the snowpack are localized (Schürholt et al., 2021). To resolve these processes on small timescales and at specific locations requires a much higher spatiotemporal resolution than is typically provided by existing operational schemes. In its original version, SNOWPACK for instance uses time steps on the order of 15 min or longer (Bartelt and Lehning, 2002) to facilitate seasonal simulation times. For Crocus time steps are on the order of 15 min (ViallonGalinier et al., 2020) to 1 h (Vionnet et al., 2012). The recent work of Jafari et al. (2020) provides a first attempt to account for vapor transport within a coupled snowpack model. In their paper, they accounted for diffusive vapor transport and phase change following Hansen and Foslien (2015) and analyzed its feedback on the snow density. In order to resolve diffusive processes, simulations were conducted at much shorter time steps of 1 min and a finer spatial resolution of approximately 0.1 cm. For comparison, a typical layer thickness in SNOWPACK is 2 cm (Wever et al., 2016) and the minimum layer thickness in Crocus is 0.5 cm (Brun et al., 1989). While the work of Jafari et al. (2020) demonstrates the general feasibility of vaporcoupled snowpack models, the exact nature of how vapor transport and phase changes interfere with stressinduced settling remains to be investigated in depth.
It is well known that any numerical strategy that aims at simulating simultaneous settlinginduced deformation of the snowpack and (arbitrary) diffusive transport requires a special computational treatment to couple both. Diffusive transport and reactive phase change are best modeled by taking an Eulerian perspective, hence on a static mesh. In contrast there exist a number of different techniques to incorporate the settlinginduced deformation. One option is to use a timedependent coordinate transformation by Morland (1982), who developed a fixed domain transformation to solve onephase diffusion problems with a moving free surface on a finite, timeinvariant computational domain. An alternative approach was put forward by Wingham (2000), who used a different spatiotemporal coordinate transformation for firn densification. Both transformation strategies effectively eliminate the vertical motion (or gradients of it) from the computational update procedure. And exactly the same is (implicitly) performed in the present treatment of densification in snowpack models (Bartelt and Lehning, 2002; Vionnet et al., 2012), where the coordinate transformation embodied in the deformation of the underlying computational grid through the update of layer positions and/or thicknesses is (implicitly) exploited for the ice mass conservation. However, the present descriptions do not take full advantage of a clear and explicit separation into a Lagrangian deformation module that accounts for mechanical settling and an Eulerian transport and phase change module. The benefit of this hybrid computational strategy is that it is easy to understand, computationally feasible, provides a modular error control and increases the interpretability by disentangling numerical artifacts from features of the underlying nonlinear process models. Hybrid numerical schemes that combine an Eulerian process model with a Lagrangiantype spatiotemporal mesh adaptation are not new. These schemes have been used in other disciplines, e.g., for phase change problems (Lacroix and Garon, 1992); as σ coordinates in oceanography, where the ocean's surface and bottom are projected onto coordinates σ = 0 and σ = −1 that follow the ocean floor's topography (Mellor and Blumberg, 1985); or for shallowflow models (Kowalski and Torrilhon, 2019).
The aim of our work is twofold: first, we will describe our numerical strategy for a phasechanging snowpack. The numerical scheme is hybrid, in the sense that it clearly discriminates between a solution of the mechanical settling operator by means of a Lagrangian approach and a solution to the transport and phase change operator by means of an Eulerian approach. To some degree, the numerical model description must be understood as a rigorous reformulation of the numerical schemes from existing computational snowpack models. Yet, in addition to existing schemes we (a) explicitly separate the Eulerian and Lagrangian part of the solver to facilitate a later modular adaption, (b) provide a full finitedifference formulation including correction terms due to the deforming (nonuniform) mesh that are typically omitted, and (c) discuss options to increase the approximation accuracy of the various parts of the numerical scheme. Second, we demonstrate the computational potential by applying and analyzing simulation results for an idealized twolayer, drysnow situation. We consider a model cascade of different process building blocks, which in their most comprehensive version, correspond to fully coupled heat and vapor transport alongside phase changes and settling.
With this work we seek to contribute to anticipated future developments of snow or firn models or likewise extensions of existing ones that aim at flexibility and modularity while providing a simple, mathematically rigorous numerical approximation for a stable and robust integration of generic multiphysics process equations. By modularity and extendability we understand the possibility of considering or neglecting specific process modules and parametrizations in a straightforward way. This modularity would enable us to (a) investigate competing nonlinear effects systematically from a cascade of process models, (b) assess the quality of the numerical approximation independently and (c) conduct a standardized model selection based on welldefined benchmarks.
The paper is structured as follows. In Sect. 2, we recall the drysnow model equations comprising the relevant transport, phase change and mechanical aspects. In Sect. 3, we introduce the Eulerian–Lagrangian numerical scheme and its solution using the method of characteristics. In Sect. 4, we present and discuss results from a number of simulation scenarios, including verification scenarios that consider transport, phase changes and mechanics in the absence of any interaction, as well as coupled scenarios that focus on their interplay. We furthermore investigate the impact of different viscosity parametrizations and assess the behavior when switching to a Glen type of nonlinear constitutive closure. Finally, we compare our results to a conventional layerbased treatment. In Sect. 5, we summarize and discuss our findings, and in Sect. 6 we draw conclusions regarding future snowpack modeling.
(Calonne et al., 2014)2.1 General situation
As a common starting point, snow models take a macroscale perspective that volume averages (Bartelt and Lehning, 2002; Bader and Weilenmann, 1992; Hansen and Foslien, 2015) or homogenizes (Calonne et al., 2014) the snowpack's microstructural state into macroscale variables. If not stated otherwise, we implicitly assume all state variables to be macroscale variables. State variables, model parameters and constants used in this paper are summarized in Table 1.
In the most general case, snow is a mixture of ice, air, vapor and water, and the snow density is given as a mixture of the respective pure densities (Bader and Weilenmann, 1992; Morland et al., 1990). The amount of ice in one reference volume of snow is ϕ_{i}ρ_{i}V, in which ϕ_{i} denotes the ice's volume fraction, ρ_{i} its pure density and V the volume of the reference volume. For dry snow, further contributions due to water vapor can be neglected, and the snow density can be approximated as ϕ_{i} ρ_{i}. The structure and volume fraction of the ice can change over time either due to straininduced settling processes or due to transient phase changes, such as sublimation and deposition or melting and freezing. Our paper focuses on the derivation of a hybrid Eulerian–Lagrangian framework to solve settling, transport and phase changes with an assessment of the computational building blocks. To this end we restrict ourselves to dry snow and allow for one secondary phase (vapor) in an Eulerian treatment coupled to the Lagrangian treatment of the ice phase. With respect to computational model development, we regard the drysnow situation as the more challenging (yet less investigated) one compared to the wetsnow situation, mostly due to a broader spectrum of characteristic spatial and temporal scales involved (see more detailed discussion in Sect. 6).
Note, that water transport and solid–liquid phase change can in principle be integrated following a similar strategy to that presented in this paper. The following section introduces the (macroscale) snowpack model where the subsection structure reflects the laterdescribed modular structure of the numerical core.
2.2 Ice mass balance
The ice volume fraction ${\mathit{\varphi}}_{\mathrm{i}}={\mathit{\varphi}}_{\mathrm{i}}(z,t)$ within a spatiotemporally evolving snowpack of varying snow height H(t) is governed by the ice mass balance and reads
with velocity field v, source term c and ice density ρ_{i} (Hansen and Foslien, 2015; Bader and Weilenmann, 1992). Note that mechanical settling is neglected in Part 1 of this companion paper. The corresponding ice mass balance (Eq. 7 in Part 1) does thus not include the velocity field v.
In a 1D situation that focuses on an evolving vertical snow column, we have the vertical position z as the only relevant spatial coordinate ($z\in [\mathrm{0},H(t\left)\right]$). The velocity field v reduces to vertical velocity $v=v(z,t)$, which depends on time and the position within the column. It is negative for snow height decrease and positive for snow height increase. Vertical motion results either from mechanical settling, hence a consolidation or compaction of the snowpack, or alternatively as a continuity response to changes in ice volume from sublimation, deposition, melting and freezing via the source term c. The continuity response leads to a minor vertical decrease/increase in snow height. Though effects due to consolidation of snow may be significantly more pronounced than those due to phase change processes in the pore space, the latter needs to be accounted for to acknowledge mass conservation of the complete system. At this point in time, we do not consider any additional increases in snow height due to precipitation, yet we discuss how this can be included in the future in Sect. 5.
The source term $c=c(z,t)$ varies with time and position in the column and stands for a gain or loss of ice mass from phase change (Bader and Weilenmann, 1992) per unit volume and unit time. As we constrain this paper to the dry situation we will henceforth refer to c as the deposition rate. c is positive (production) if new ice is built, namely vapor deposits, and it is negative (loss) if ice is lost, namely sublimates. Finally, ρ_{i} denotes the constant pure density of ice and serves as a scaling factor. The ice mass balance (Eq. 1) couples mechanical settling and phase change processes. Considering the equation in its full form is essential for our goal to model and eventually analyze the interplay between these processes. The structure of the ice mass balance resembles an advection–reaction equation that can conveniently be solved by means of Lagrangiantype computational methods, such as the method of characteristics (see Sect. 3). Yet in order to do so, we need to provide a closure for both vertical velocity v and deposition rate c.
2.3 A closure for the velocity field
Velocity v represents mechanical deformation in the snowpack. Its idealized relation to the strain rate is given by
Note that this is simplified with respect to more general, tensorial formulations of 1D consolidation theories; see for instance Audet and Fowler (1992). Yet even the idealized formulation Eq. (2) will be sufficient for our purposes, as it resembles the approach implicitly chosen in snowpack models (Bartelt and Lehning, 2002; Vionnet et al., 2012).
In general one would expect that porous snow inherits the nonlinear constitutive behavior of ice (Kirchner et al., 2001), which leads to
which is a variant of Glen's law. Here, η denotes the compactive viscosity of snow and σ denotes the stress. The choice of the Glen exponent m in earlier work depends on both the physical regime and the computational feasibility. The linear form of Glen's law (m=1) is chosen in Vionnet et al. (2012) and Bartelt and Lehning (2002). For the sake of comparability we thus mainly use a linear version of Glen's law; hence m=1. Our framework, however, also copes with the nonlinear relation, such as m=3, and we later include a comparative example.
The compactive viscosity η depends on the snow's microstructure and is challenging to determine from experiments (Wiese and Schneebeli, 2017). It is typically provided as a parametrized closure for a specific physical situation and strongly correlates with the choice for the Glen exponent m. This fact clearly constrains its universal applicability and makes any transfer of a validated snowpack model to other physical situations challenging. In this article, we will consider both constantviscosity scenarios, as well as an additional scenario with a varying viscosity assuming an empirical viscosity closure from Vionnet et al. (2012)
with the state variables temperature T and ice volume fraction ϕ_{i}; the constants ice density ρ_{i} and phase change temperature T_{ph} = 273 K; and further constants η_{0} = 7.62237 × 10^{6} kg s^{−1}, a_{η} = 0.1 K^{−1}, b_{η} = 0.023 m^{3} kg^{−1} and c_{η} = 250 kg m^{−2}. Finally, f reflects properties of the snow microstructure, i.e., the angularity and the size of the grains, and it is assumed to be 1 in our case. The constantviscosity value applied to a linear Glen's law ${\mathit{\eta}}_{\text{const},m=\mathrm{1}}$ is derived with intermediate values for the ice volume fraction and temperature of the respective initial conditions. These values are plugged into the empirical closure Eq. (4) to solve for viscosity. The same procedure cannot be applied to derive a constantviscosity value for the nonlinear version of Glen's law ${\mathit{\eta}}_{\text{const},m=\mathrm{3}}$ since the viscosity closure (Eq. 4) was initially calibrated to the linear form of Glen's law (m=1). Instead we choose a snow deformation rate from the literature ($\dot{\mathit{\u03f5}}$ = 10^{−6} s^{−1}; Johnson, 2011) and determine the maximum stress value from the initial snow density. These strain rate and stress values are then inserted into the constitutive relation (Eq. 3), which is finally solved for viscosity. To avoid infinite ice volume growth above physical values (ϕ_{i} > 1), the viscosity must tend to infinity for ϕ_{i}→1. Therefore, the constantviscosity values are restricted to ice volumes below 0.95 by multiplication with an icevolumefractiondependent power law (Appendix Eq. A6). This power law yields ∼ 1 for ϕ_{i}≤0.95 and exponentially increases for higher ice volumes. Multiplied with the constantviscosity values, viscosity remains constant below ϕ_{i} < 0.95 and exponentially increases above it, which stops further densification and settling. This procedure does not intend to reproduce the correct physics for lowporosity ice, although it mathematically leads to a similar crossover behavior.
In the absence of strong horizontal deformation and deviatoric stress components, it is reasonable to assume a stressfree condition at the snow's surface and a hydrostatic stress condition in its interior:
g is the gravitational acceleration, and ρ_{snow} refers to the snow's density, which is clearly dominated by the ice fraction via ρ_{snow}≈ϕ_{i}(z)ρ_{i}. It varies with the position z in the snow column due to a vertically varying ice volume fraction ϕ_{i}(z). Integration of Eq. (5) and combination with Eqs. (2) and (3) yields an expression for the velocity gradient:
ζ is the integration variable. A second integration along the vertical axis finally yields an expression for the velocity at position z in the snow column:
in terms of total height H(t) and the ice volume fraction ϕ_{i}(z,t) and with $v(z=\mathrm{0},t)\equiv \mathrm{0}$. This definition of the vertical velocity yields a process that complies with the obvious physical constraints: (a) the velocity vanishes at the bottom of the snow column, hence preventing artificial penetration into the ground. This is similar to displacement requirements in SNOWPACK (Bartelt and Lehning, 2002). (b) The vertical velocity accumulates with height, which prevents any artificial disaggregation of the snowpack. (c) The vertical velocity relaxes towards zero as the ice volume fraction tends towards its maximum volume fraction ${\mathit{\varphi}}_{\mathrm{i}}<{\mathit{\varphi}}_{\mathrm{i},max}<\mathrm{1}$. In the remainder of this paper, we will use Eq. (7) to account for the mechanical settling of the snowpack.
2.4 Transport and phase changes
The ice deposition rate c as relevant to solve Eq. (1) typically depends on a cascade of coupled heat and mass transport for the involved phases of ice, water and vapor. In this article, we will consider a process model proposed by Hansen and Foslien (2015) that reflects a drysnow condition in which void space is filled by vapor only. Note, however, that this coupled process model could readily be substituted or extended by another one, e.g., from Calonne et al. (2014), Jafari et al. (2020) or Schürholt et al. (2021).
Next, we state the essential aspects and process equations of the model proposed in Hansen and Foslien (2015) and describe how it can be used to recover the ice deposition rate.
Assuming a drysnow condition, the ice production is solely determined by mass transport between vapor and ice. The vapor mass balance reads
in which ρ_{v} denotes the vapor density and D_{eff} the effective vapor diffusion coefficient. Vapor production corresponds to a negative ice deposition rate −c that represents sublimation. Following Hansen and Foslien (2015), vapor density in the pore space can be assumed to be at saturation density ${\mathit{\rho}}_{\mathrm{v}}^{\text{eq}}$, so ${\mathit{\rho}}_{\mathrm{v}}\equiv {\mathit{\rho}}_{\mathrm{v}}^{\text{eq}}$. The latter is well investigated, and empirical relations exist that specify its temperature dependency ${\mathit{\rho}}_{\mathrm{v}}^{\text{eq}}\left(T\right)$. In this work, we will employ an empirical relation from Libbrecht (1999). The full expression can be read in Appendix A1. Due to the closure for vapor density ${\mathit{\rho}}_{\mathrm{v}}^{\text{eq}}$, the vapor mass balance (Eq. 8) can be rewritten using the temperature dependence of the equilibrium vapor density:
Assuming the snow to be in thermal equilibrium at the microscale, we can likewise write the energy balance in terms of the temperature, which reads
The parameters (ρC)_{eff} and k_{eff} stand for the effective heat capacity of snow and effective thermal conductivity, respectively. Both parameters depend on the ice volume fraction, and their definition is stated in Appendix A2. The righthand side of the heat equation (Eq. 10) accounts for latent heat release, which is coupled to phase change processes.
The system of the two equations, Eqs. (9) and (10), and the two unknowns, temperature T and deposition rate c, is solved by replacing c in Eq. (10) with Eq. (9), which yields a nonlinear equation for temperature:
The spatiotemporal temperature evolution is then used to recover the ice deposition rate c from either Eq. (9) or Eq. (10).
The complete process model is now given by the ice mass balance Eq. (1), its mechanically induced vertical velocity Eq. (7), and the coupled system for temperature Eq. (11) and ice deposition rate determined by either Eq. (9) or Eq. (10). Each of the equations will be solved in a separate module. The ice mass balance in conjunction with the vertical velocity has the form of a nonlinear advection equation, whereas the remaining equations are of parabolic nature, which is reflected in our general approach to solve the system.
3.1 General approach to the computational strategy
Based on the distinction into diffusion and advectiondominated processes, we propose a twostep solution scheme:
Step 1. This step accounts for the mesh deformation and solves the advectiondominated mechanical settling, i.e., the ice mass balance Eq. (1), by means of a Lagrangian approach that tracks the movement of the coordinates including changes from metamorphism.
Step 2. This step determines the spatiotemporal evolution of temperature and deposition rate fields as introduced in Sect. 2.4 based on an Eulerian approach that solves the diffusiondominated transport and phase changes via a finitedifference implementation on a deformed (unstructured) mesh.
Note that here we employed a finitedifference method because it provides a feasible algorithm that is applicable to the scenarios considered in the paper using a 1D snow column. It also naturally integrates with the Lagrangian part of the solution (Step 1), as we can reuse the same mesh. In principle, it is also possible to couple the twostep approach with a finiteelement solution for the temperature and deposition rate, for instance when aiming for a 2D or 3D model in a complex geometry that incorporates realistic mountain slope topographies. When using a finiteelement solver, however, we have to keep in mind that the deposition rate and temperature fields need to be reconstructed from the solution at each time step. Especially when wanting to use higherorder elements, this might limit computational feasibility.
Our solution scheme alternates both steps via straightforward firstorder operator splitting. This is found to work well for our simulation scenarios yet could be readily exchanged with a higherorder splitting scheme, e.g., a secondorder Strang splitting (LeVeque, 2002), if required.
The computational model is implemented in Python, and it is modular and extendable, in the sense that each module can be separately activated and deactivated. This not only simplifies the verification of individual process building blocks but also allows an indepth investigation of the various coupling effects and the model's nonlinear feedback. Alternative formulations e.g., of the parametrized velocity field are implemented and can easily be exchanged. Finally, the modular structure facilitates the implementation of additional closure relations or the integration of entire new process modules.
3.2 Computational grid
In this paper, we consider a 1D snow column, which is discretized into nz+1 spatial mesh nodes denoted by z_{k} with $k\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{\dots},nz\mathit{\}}$. We applied 101 computational nodes (nz=100) except for some simulations that required a higher resolution of 251 nodes (nz=250). The mesh is nonuniform in general, meaning that the distance between neighboring nodes ${z}_{k+\mathrm{1}}{z}_{k}$ varies throughout the snow column and with time. Note that the z axis is oriented opposing gravitational acceleration, such that z_{0} denotes the position of the ground and z_{nz} the position of the snowpack's free surface. Time increments are denoted by t_{n} with $n\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{\dots},nt\mathit{\}}$ and nt being the maximum number of time steps in a complete simulation run. For each of the field variables subscript k denotes the vertical coordinate and superscript n denotes the time step; hence $T({z}_{k},{t}_{n})={T}_{k}^{n}$.
3.3 Lagrangian solution of the ice mass balance
When the snowpack is subject to vertical motion, e.g., settling, its physical height decreases; hence its vertical extent shrinks. One option to reflect this in a computational method is to adjust the spatial node coordinates accordingly. The challenging fact in our situation is that the vertical motion within the snow column (nonlinear advection) is coupled to phase changes, i.e., a change in the ice volume fraction via the source term in the ice mass balance (Eq. 1). The method of characteristics is a suitable method to solve such a nonlinear advection equation with source terms. It can be interpreted as a simultaneous motion tracking of snow reference volumes, referred to as the integration along socalled characteristics, while also accounting for its metamorphism along the trajectory. By construction, the method correctly tracks the snowpack's moving free surface. Due to the fact that the snow column's evolution is determined with respect to a reference volume that moves vertically at speed v in the snowpack, the method of characteristics is called a Lagrangian approach.
In order to derive the specific update rule for the ice mass balance Eq. (1), we first apply the product rule to its initial Eulerian version
and then reformulate the equation in a Lagrangian reference frame, hence with respect to nodes moving at the vertical velocity v. Changing to the moving reference frame effectively compensates for the advection term in Eq. (12) and yields
Equation (14) accounts for the settling of material particles within the snowpack. We will use it to update the coordinates of the mesh nodes directly, which results in a continuous mesh deformation as illustrated in Fig. 1. Equation (13) captures the evolution of the ice volume fraction along the trajectory of a moving ice volume within the snowpack. It accounts for volume changes due to (a) mass production and loss in response to phase changes and (b) vertical variation in the vertical velocity. Further details and generalizations of the method of characteristics can be found in Farlow (1993).
Equations (13) and (14) can be solved analytically for a constant vertical velocity and deposition rate. In our case however, the velocity closure is provided by Eq. (7) and the deposition rate results from solving yet another process model (Eqs. 10 and 11), which requires a numerical solution. Since we expect the response of the ice volume fraction to be slow (with respect to other processes in the system), we will rely on a firstorder explicit Euler time integration scheme:
In order to update the mesh coordinates according to Eq. (16) for the vertical velocity closure derived before, we need to numerically approximate Eq. (7) at each node z_{k}, which results in
with $\mathrm{\Delta}{z}_{j}^{n}:={z}_{j+\mathrm{1}}^{n}{z}_{j}^{n}$ where $j\in [\mathrm{0},nz)$, m being the Glen exponent, η being viscosity and σ_{j} denoting the stress exerted by the overburdened snow mass
where g is gravitational acceleration. Note that the stress at the uppermost node k=nz is zero, so velocity v(z_{nz}) is only controlled by the movement below and is thus equivalent to the velocity at the next lower node (v(z_{nz−1})). The forward Euler scheme of Eqs. (15) and (16) via the method of characteristics combined with the velocity update (Eq. 17) essentially resembles the treatment of mass conservation as it is, for instance, presently carried out in SNOWPACK. However, the explicit formulation and numerical treatment of Eqs. (15) and (16) allows us to also employ other (e.g., higherorder, implicit) solution schemes for both equations if this is required to capture detailed aspects of the spatiotemporal coupling of phase changes (c) and settling (via ∂_{z}v) (cf. also the discussion in Sect. 5). To solve Eq. (15), we directly discretize the velocity's spatial derivative ∂_{z}v, which corresponds to the strain rate ${\dot{\mathit{\u03f5}}}_{k}^{n}=\frac{\mathrm{1}}{\mathit{\eta}}{\mathit{\sigma}}_{k}^{m}$ given via Eq. (3). This is beneficial, as it avoids numerically approximating the velocity gradient. The complete numerical update of ice volume fraction ϕ_{i} and mesh coordinates z can now concisely be written as
Similarly to existing layerbased schemes (see for instance Sect. 3.4. in Bartelt and Lehning, 2002, or its recent extension in Jafari et al., 2020), the method of characteristics provides information on the settling of layers within the snowpack. Yet, in addition, it serves as a basis for a fully modular and flexible computational strategy that (a) by construction accounts for the twoway feedback between the ice volume fraction and mass production or decay rates resulting from phase changes as a response to transport processes within the snowpack and (b) allows for a flexible adoption and extension of the process model (used to determine c) and the velocity closure. The latter could for instance serve as a pathway to integrate a datadriven velocity closure (or assimilation) from measurements. Such flexibility in numerical tools will be important in the future to conduct model comparisons, such as presented in Schürholt et al. (2021) within holistic snowpack models, or even a formalized Bayesian model selection that allows for inferring the most plausible process model out of a pool of candidate models given certain data. A remaining difficulty now is to provide a (Eulerian) numerical scheme for diffusive processes that can operate on a spatially varying unstructured mesh.
3.4 Eulerian solution of transport and phase changes on a moving mesh
The process model accounting for vapor transport and heat transport (Eqs. 9 and 11) has to be solved with respect to a moving computational mesh according to Eq. (16). Both equations have the same generic structure; namely
with $\mathit{\alpha}={\mathit{\alpha}}_{\mathrm{T}}=(\mathit{\rho}C{)}_{\text{eff}}+(\mathrm{1}{\mathit{\varphi}}_{\mathrm{i}})\frac{\mathrm{d}{\mathit{\rho}}_{\mathrm{v}}^{\text{eq}}\left(T\right)}{\mathrm{d}T}L$, $\mathit{\beta}={\mathit{\beta}}_{\mathrm{T}}={k}_{\text{eff}}+L{D}_{\text{eff}}\frac{\mathrm{d}{\mathit{\rho}}_{\mathrm{v}}^{\text{eq}}\left(T\right)}{\mathrm{d}T}$ and $\mathit{\gamma}={\mathit{\gamma}}_{\mathrm{T}}=\mathrm{0}$ for the heat equation (Eq. 11) and $\mathit{\alpha}={\mathit{\alpha}}_{c}=(\mathrm{1}{\mathit{\varphi}}_{\mathrm{i}})\frac{\mathrm{d}{\mathit{\rho}}_{\mathrm{v}}^{\text{eq}}\left(T\right)}{\mathrm{d}T}$, $\mathit{\beta}={\mathit{\beta}}_{c}={D}_{\text{eff}}\frac{\mathrm{d}{\mathit{\rho}}_{\mathrm{v}}^{\text{eq}}\left(T\right)}{\mathrm{d}T}$ and $\mathit{\gamma}={\mathit{\gamma}}_{c}=\phantom{\rule{0.125em}{0ex}}c$ for the vapor transport equation (Eq. 9).
An implicit firstorder finitedifference approximation of Eqs. (9) and (11) for a spatially varying mesh of increments $\mathrm{\Delta}{z}_{k}^{n}$ results in
Note that parameters α_{f} and β_{f} for $f\in \mathit{\{}T,c\mathit{\}}$ also vary in space and time and will be (explicitly) evaluated based on the snowpack's state at time n. The terms E_{c} and E_{T} are higherorder mesh errors for the vapor and temperature equations. These higherorder mesh errors account for the necessary correction due to nonuniformity of the mesh and are controlled by the temperature gradient; they vanish for equidistant meshes or constant temperatures. The complete form of the higherorder mesh errors is given in Appendix B, and their effect on the accuracy of the simulation is discussed in Sect. 4.3.
The complete numerical update can be concisely written in matrix form, which matches with the way it is implemented in the software:
First, Eq. (24) is solved for temperature T^{n+1}. Next, the updated temperature is used to solve Eq. (25) for the deposition rate c^{n+1}. The complete matrix definitions are given in Appendix C. Note that, formally, it would be possible to add up matrices A_{T} and E_{T} as well as A_{c} and E_{c}. We decided to keep them in this particular form to stress the similarity of this formulation with a standard finitedifference approximation on an equidistant mesh, in which we are left with B_{T} and B_{c} and E_{T} and E_{c} vanish.
3.5 Iterative coupling of Eulerian and Lagrangian solutions
The derived numerical update routines for temperature, deposition rate, vertical velocity and ice volume fraction comprise the four main modules that are sequentially called to update the respective state variables for one time step. A schematic illustration is given in Fig. 2. The equations for heat and vapor transport have already been implemented by Calonne et al. (2014) and Hansen and Foslien (2015). A feedback on the ice volume fraction in the absence of a vertical velocity has been investigated in Part 1 of the companion paper (Schürholt et al., 2021). The modules for vertical velocity and the coupled update of ice volume fraction and mesh coordinates, through the method of characteristics, are novel in our approach. Our implementation is modular in the sense that it allows for a coupling with other process models that comply with a nonuniform mesh.
The time step size for the next time step n+1 is dynamically updated in the computational scheme. Since diffusive processes are dominant, we utilize the mesh Fourier number based on the diffusivity $\frac{{\mathit{\beta}}_{\mathrm{T}}}{{\mathit{\alpha}}_{\mathrm{T}}}$ of heat of the current time step n:
Since this choice for the time step computation did not yield instabilities, we excluded the vapor's diffusivity for the time step computation. Note that in response to settling processes, the mesh sizes vary and decrease (see Fig. 1) with time, and so does the time step.
In the following, we describe how the modularity of the model is applied and used to assess the individual effect of the different process building blocks by a strategical activation and deactivation of the modules.
3.6 Application of the model
We applied the developed numerical scheme to perform several simulations with varying combinations of activated and deactivated advection and diffusiontype process building blocks, e.g., transport and phase changes, such as those also considered in Schürholt et al. (2021), vs. transport and phase changes in the presence of settling and their corresponding coupled scenarios. Furthermore, this scheme allowed the numerical verification of separate building blocks. While the scenarios are still idealized, they demonstrate the robustness of the Eulerian–Lagrangian scheme against the selection of varying subsets of model components. Table 2 provides an overview of the various combinations we considered as they have been introduced in Sect. 3. Note that we use the terms vertical velocity and settling velocity interchangeably.
Firstly, we focus on the effects due to pure mechanical settling on the snowpack (Case 1). Next, we consider isolated heat transport (Case 2) as well as its interplay with settling processes (Case 3). Similarly, we consider coupled heat and vapor transport first in the absence of settling (Case 4) and later with settling (Case 5). For Case 5, we evaluate the effect of included or excluded higherorder mesh errors E_{T} and E_{c} (see Sect. 3.4) on the temperature profiles. Case 1, Case 3 and Case 5 consider the constant viscosity for a linear Glen's law (m=1) ${\mathit{\eta}}_{\text{const},m=\mathrm{1}}$, as introduced in Sect. 2.3. Furthermore, we investigate the impact of an empirical, temperaturecontrolled and icevolumefractioncontrolled viscosity closure (Eq. 4), first on settling only (Case 6) and then on the fully coupled processes (Case 7). Next, we show that our general approach can be combined with the nonlinear Glen's law (Eq. 3) by using m=3 (Case 8) and an accordingly adjusted constant viscosity ${\mathit{\eta}}_{\text{const},m=\mathrm{3}}$. For a detailed explanation of the general derivation of the viscosity values see Sect. 2.3. Lastly, we compare our new modeling approach to that of layerbased schemes (cf. Sect. 1).
3.7 Computational setup, initial and boundary conditions
Initial condition. The initial ice volume fraction ϕ_{i} reflects a layered situation as depicted in Fig. 3, with two snow layers of equal thickness. The bottom layer has an initial snow density of 150 kg m^{−3}, and the upper layer's density is 75 kg m^{−3}. The transition from the upper layer to the lower layer is linearly smoothed out over 2 cm, which for a grid defined according to Sect. 3.2 corresponds to 5 computational nodes for the coarser and 11 computational nodes for the finer discretization. The snow densities are in the range of “damped new snow” and “new snow”, respectively (Paterson, 1994). Snow densities in this range are expected, e.g., for new snow in the European Alps (Helfricht et al., 2018) or in the Rocky Mountains (Judson and Doesken, 2000). We choose this layered snowpack to ensure an extreme and very active snow regime with a strong dynamical coupling of the processes. Temperature is initially constant at 263 K throughout the whole snowpack. The deposition rate is directly deduced from temperature (see Eq. 23) and therefore requires neither initial nor boundary conditions. From the initial condition we derived the constantviscosity values ${\mathit{\eta}}_{\text{const},m=\mathrm{1}}\approx $ 9.1 × 10^{7} Pa s and ${\mathit{\eta}}_{\text{const},m=\mathrm{3}}\approx $ 16 × 10^{12} Pa s; see also Appendix A3.
Boundary condition. We consider a constant temperature of 273 K at the bottom boundary and a constant temperature of 253 K at the free surface.
Simulation time. We simulate 2 d (48 h), 3 d (62 h) and 4 d (96 h) scenarios.
4.1 Settling (Case 1)
First, we investigate the effects of mechanical settling on the snowpack (Case 1 in Table 2) and in particular the evolution of the vertical velocity (Fig. 4a) and the ice volume fraction (Fig. 4b). The vertical velocity decreases from top to bottom and relaxes during the first 48 h. Vertical velocity varies more in the lower layer compared to in the upper layer within one time step. This pattern remains prominent as time proceeds while the overall velocity variation decreases. This effect is due to the increase in the overburdened snow mass from top to bottom. Settling proceeds the fastest just after the start of the simulation, when the snowpack is at maximum height, and correspondingly its snow density is the lowest. In the course of time the ice volume fraction increases faster in the lower layer than in the upper layer, and it is the highest at the bottom of the snowpack (Fig. 4b). This observation is also visualized in Fig. 5, which depicts profiles of snow density. Furthermore, the extent of the upper layer decreases only slightly, approximately 3.5 cm, over the simulation time, whereas the lower layer reduces to half of its initial height (approximately 12.5 cm). These effects are expected and reflect the correlation between the amount of compaction and the total overburdened mass. The total settlement of 14 cm after 2 d means a 30 % snow height reduction. Bartelt and Christen (2007) simulated an 11.6 to 54.8 cm snow height reduction after 5 d for an initially 90 cm high snowpack of 115 kg m^{3} density, which is a snow height reduction of 12 % to 60 %. Taking into account that snow settles more slowly with increasing density, our results fit to the highest settling rate derived by Bartelt and Christen (2007).
4.2 Heat transport in the absence and presence of settling (Cases 2 and 3)
In this subsection, we first consider isolated heat transport. This simulation scenario refers to Case 2 in Table 2. Temperature (Fig. 6a) and temperature gradients (Fig. 6c) reach a stationary state after approximately 60 h. Heat flux differences between the two layers are clearly visible in the temperature gradient plot. Next, heat transport is superposed by mechanical settling (Fig. 6b and d), representing Case 3. As a result, snow height decreases while the internal temperature profiles evolve. Active mechanical processes yield a steepened temperature gradient and hence a higher value of the heat flux (Fig. 6d). This effect can be attributed to

the decrease in snow height while keeping the temperatures at the boundaries fixed and

the permanent change in thermal conductivity and thermal diffusivity due to their dependency on variations in the ice volume (Eq. A3).
The temperature profile will reach the stationary state once the ice volume fraction has reached its maximum value.
4.3 Heat and vapor transport in the absence and presence of settling (Cases 4 and 5)
By using the vapor formulation from Hansen and Foslien (2015), transport of vapor through and phase changes in the snowpack both require an apparent temperature gradient such that the evolution of vapor transport can only be considered in conjunction with heat transport. In Fig. 7a, we compare the deposition rate (negative for sublimation) due to heat and vapor transport only (Case 4 in Table 2) with the deposition rate obtained when considering additional settling processes, representing the fully coupled processes (Case 5 in Table 2). Both profiles are characterized by moderate deposition rates throughout the snow column with a pronounced negative (sublimation) peak at the center of the snow column, which is located in the transition area of the layers. The profile for the fully coupled processes shows a higher sublimation peak (approximately 4 times higher). Figure 7b shows the time evolution of the fully coupled processes (Case 5). In the first hours, sublimation is low in the transition area. After approximately 6 h, the pronounced sublimation rate peak, as already described for (Fig. 7a), develops and increases until the end of the simulation (48 h). The increased sublimation in the layer transition area may be driven by strong vapor density gradients (Fig. 7c) above the transition area that can be inferred from a strong, local temperature gradient (Fig. 7d). This temperature gradient is further enhanced (Fig. 6d) by compaction due to settling for Case 5, which yields even stronger variations in the material properties in the transition area than without compaction and explains the stronger sublimation rates for the fully coupled processes.
Furthermore, both profiles in Fig. 7a show a small peak in the deposition rate (positive x direction) just above the aforementioned sublimation rate peak. This peak is very weak for Case 4 and more prominent for Case 5. This deposition rate peak is highly interesting as it is interpreted as the onset of spatiotemporal oscillations as observed and investigated in greater detail in the companion paper (Schürholt et al., 2021). Schürholt et al. (2021) describe these wiggles as “smooth oscillations” that are “intrinsic features” of the equations. The results in Fig. 7a nicely demonstrate that (a) our Eulerian–Lagrangian scheme can capture this behavior and (b) the instability prevails and even increases in the presence of settling processes. The results suggest that mechanics likely increase local phase change activity in the vicinity of layer boundaries, which potentially has a large effect on weak layer formation.
The deposition rates obtained with our model are between −2 and 2 $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$, which fits to the range of −1.728 to 1.728 $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ presented in Jafari et al. (2020). Sublimation rate peaks on the order of 0.1 to 1.2 $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$ have also been computed with the numerical test cases by Hansen and Foslien (2015). For comparison with experiments, deposition rates can be derived via $\text{SSA}\cdot {v}_{\mathrm{n}}\cdot {\mathit{\rho}}_{\mathrm{i}}$, with v_{n} being the ice crystal's interface growth velocities and SSA the ice's specific surface area (see Calonne et al., 2014, Eq. 21 therein). For a simple characteristic scale analysis, we considered SSA in the range of 0.6 × 10^{4} to 1×10^{4} m^{−1} (Schleef et al., 2014) and an interface growth velocity on the order of 1 × 10^{−9}m s^{−1} (Krol and Löwe, 2016; Calonne et al., 2014). Combining these literature values yields deposition rates on the order of 0.5 × 10^{4} $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{d}}^{\mathrm{1}}$, which is significantly larger than our simulation results. Interface growth velocities on the order of 10^{−13} or 10^{−14} m s^{−1} would match with the simulated magnitudes for the deposition rate.
Lastly, we evaluate the impact of included higherorder mesh errors E_{T} and E_{c} (see Sect. 3.4) on the temperature distribution. We determine the error by computing the temperature deviation between the solution that considers higherorder mesh errors and the solution that does not. The deviation is then quantified in an L1 norm. The error increases with simulation time and is 0.13 K after 24 h, 0.23 K after 36 h and 0.28 K after 48 h. After 48 h the deviation is highest for the computational nodes just above the layer transition, where high temperature gradients are present (see Fig. 6). Note that the error for the deposition rate could be derived similarly. From the temperature error, the deposition rate error can be derived as the deposition rate is directly derived from temperature via the vapor transport equation; we consider one error measure as sufficient to emphasize the impact of mesh errors.
4.4 Settlinginduced evolution of the ice volume fraction in the absence and presence of transport (Cases 1 and 5)
In this section, we compare isolated settling (Case 1 in Table 2) and the fully coupled processes (Case 5 in Table 2) with respect to their impact on the evolving ice volume fraction. Figure 8 shows the corresponding ice volume fraction profiles after 2 d. Both profiles are very similar (Fig. 8a), which suggests that the density evolution is dominated by settling processes and coupled heat and vapor transport play a minor role. When focusing on the upper boundary of the transition area (Fig. 8b), we find however a locally decreased ice volume fraction for the fully coupled processes (Case 5). This suggests a local ice volume decay for active vapor transport and implies phase changes. This observation is consistent with the enhanced sublimation rate observed in Fig. 7 and indicates the formation of a density heterogeneity.
4.5 Heat and vapor transport coupled to settling with a dynamic viscosity (Cases 4, 6, and 7)
Figure 9 shows the evolution of the ice volume fraction and viscosity over 4 d for Case 7 (Table 2), which is the fully coupled processes coupled to dynamic viscosity. The ice volume fraction increases gradually throughout the snow column (Fig. 9a). In Fig. 9b, we see that the viscosity of the upper layer has smaller values, and they also increase more slowly compared to viscosity values in the lower layer. In contrast, viscosity increases by approximately 1 order of magnitude in the lower layer. This derives from the applied viscosity formula that is controlled by the variables of temperature and ice volume fraction. Based on the formula, viscosity varies more with respect to ice volume fraction changes than to temperature changes. The ice volume fraction varies more in the lower layer, which then also yields more variation in viscosity. Additionally, the height of the lower layer decreases less than that of the upper layer. This outcome may be related to the lower layer's higher viscosity (higher resistance to deformation).
Figure 10a compares the ice volume fraction profiles after a 3 d simulation time of Case 7 to those of Case 6 (fully coupled processes, nondynamic viscosity). For the dynamic viscosity, the ice volume fraction is higher in the lower layer and lower in the upper layer compared to that of constant viscosity. This is due to the dynamic viscosity's temperature dependence. Temperatures at the bottom are close to the melting point and yield lower viscosities. Thus, settling proceeds faster and compaction is stronger in the lower part. The opposite is true for the upper layer. Since we used an intermediate value of 263 K to derive the constant viscosity, variations in the center of the snowpack are less pronounced. Figure 10b shows the deposition rate of Case 7 compared with Case 4, which refers to deactivated settling. Similarly to Fig. 7 both deposition rate profiles have a sublimation rate peak in the transition area. Dissimilar to Fig. 7b is that the peaks are approximately at the same normalized snow height and that the peak of the fully coupled processes is not higher than the one of deactivated settling. Instead Case 7 shows less deposition in the vicinity of the transition area above the sublimation rate peak compared to Case 4. This suggests that the sublimation rate peak is less pronounced when coupled to the proposed dynamic viscosity, but settling still has an effect on the deposition rate. Additionally, Case 7 shows the small peak in deposition rate just above the sublimation rate peak, which is similar to Case 4 and discussed in Sect. 4.3.
4.6 Nonlinear Glen's law in a fully coupled drysnowpack model of constant viscosity (Case 8)
In this simulation scenario, we present the results of the fully coupled processes for the nonlinear Glen's law (Eq. 3 with m=3, Case 8). As discussed before (Sect. 2.3), the viscosity closure (whether it is a constant value or an empirical closure) strongly depends on the choice of the Glen parameter m. This requires us to adjust the constantviscosity value accordingly; see Sect. 2.3 for details.
Figure 11a (Case 8 in Table 2) shows vertical velocity profiles and the evolution of the ice volume fraction with time. The vertical velocity is almost constant in the upper layer and then decreases in the lower layer. This effect is similar to the vertical velocity profiles as presented and explained for the linear version of Glen's law (Fig. 4a), but it is more pronounced due to the nonlinearity in the constitutive law. Compared with previous scenarios, the overall vertical velocity is lower. This is probably related to the magnitude of the constant viscosity and cannot be directly related to the nonlinear constitutive law. A further sensitivity study in the future would be most informative.
In Fig. 11b the upper layer's ice volume fraction and thickness remain almost constant with time. In contrast, the lower layer decreases 9 cm in height while the ice volume fraction increases with time from top to bottom.
As shown for Case 7 (Fig. 7a) the deposition rate profile shows a sublimation peak in the layer transition area (Appendix D1) that increases with time. Overall, however, deposition rates tend to be lower compared to those of preceding computations. The reduced phase change activity in the layer transition area can be directly related to smaller vertical variations in the temperature profile. This effect may be due to less variation in the vertical velocities that yield a more uniform deformation and a less pronounced variation in the ice volume fraction across the layer transition area.
4.7 Comparison against layerbased schemes (based on Case 6)
In this section, we compare results of our proposed Eulerian–Lagrangian scheme with conventional layerbased models. We would like to emphasize that a twolayer snowpack model certainly constitutes an extremely simplified case, as layerbased schemes are usually operated with a significantly higher number of snow layers. Yet it is informative to conduct this analysis to point out differences, as these can certainly accumulate during long simulation times.
In layerbased snowpack models state variables are assigned layerwise, and the twolayer snowpack (Fig. 3) would have three computational nodes at the following locations: at the bottom of the lower layer, at the top of the lower layer and at the top of the upper layer. The two nodes located at the top of the lower and upper layers would then represent the physical state of the lower and the upper layer, respectively. Velocity is again derived from stress exerted by the overburden snow mass. Since the upper layer is represented by the computational node at the top, it is unloaded and requires a special treatment for stress. We adopt the approach by Vionnet et al. (2012) and apply a “nonphysical stress” equivalent to half of the layer's own weight, yet we apply it to the uppermost computational node (Sect. 3.4 in Vionnet et al., 2012). Next, vertical velocity is computed likewise with Eq. (7) and viscosity with Eq. (4). We compare both approaches based on Case 6 of Table 2, hence in the presence of mechanical settling and for a dynamic viscosity closure. Since we neglect heat and vapor transport, the viscosity changes over time are solely due to the evolution of the ice volume fraction alone.
In Fig. 12, we see that the layerbased scheme sustains a layerwise vertical velocity (Fig. 12a) and ice volume fraction evolution (Fig. 12c): one value for the velocity and one value for the ice volume fraction describe an entire layer. In contrast, using the generalized Lagrangian approach described in Sect. 3, we yield a sublayer resolution of the vertical velocities (Fig. 12b) and ice volume fractions (Fig. 12d). For both approaches (layerbased and Eulerian–Lagrangian) the vertical velocity is higher in the top part of the snowpack and zero at the bottom. For early times, the layerbased scheme determines a vertical velocity that is 1 order of magnitude higher than values computed with the Eulerian–Lagrangian scheme. This may be related to the comparably high (nonphysical) stress at the top of the upper layer. At the end of the simulation, the snowpack has settled almost twice as much with the layerbased scheme, which highlights the impact of this conceptual difference. This effect may result from an overestimation of velocity with layerbased schemes. Following our proposed method, the ice volume fraction is higher in the lower part of the snowpack and reaches higher values (Fig. 12d). Furthermore, the ice volume fraction at the top of the snowpack does not change during the simulation since there is no stress from overburden mass. In contrast, for the layerbased scheme ice volume fraction grows at this location (Fig. 12c). This is again due to the chosen stress condition at the top. Of course this discrepancy becomes smaller as we increase the number of layers, and this effect may reduce. However this slight offset in the stress condition will always be present and lead to uncertainties. In the proposed computational approach the spatial resolution of processes can be easily changed to assess its impact on snowpack evolution. In a future study, it might be interesting to quantitatively compare results against Jafari et al. (2020), who also rely on a rather fine spatial resolution.
4.8 Thin layers at the top and bottom of the snowpack (Case 5)
For the final scenario we implement a thin layer at the top and at the bottom of the twolayer snowpack. With this test case we want to show our model's feasibility for a potential future comparison with operational snowpack models such as Crocus that sustain layers at the top and bottom of the snowpack (Vionnet et al., 2012).
The initial condition for snow density is in principle equivalent to Fig. 3 except that the upper and lower 2 cm now form a new layer each, of 50 and 200 kg m^{−3}, respectively. The transition to the neighboring layer is linearly smoothed out over 1 cm for both thin layers. Figure 13a shows the ice volume fraction profiles for three times. The profile for 0 h reflects the initial condition for the ice volume fraction. The three layer transitions are discernible as steps in the profiles. Figure 13b depicts the deposition rate evolution. Sublimation is stronger at the layer transitions but shows a different evolution at all transitions. While the uppermost layer transition remains at a constant sublimation rate with time, sublimation at the lowermost layer transition is very strong in the beginning and then decreases. The central layer transition shows that sublimation continuously increases until the end of the simulation time. This is consistent with our observations in Fig. 7b. Figure 7a and b show that after 48 h the lowermost thin layer has been reduced to less than half its initial height while the upper thin layer's thickness has remained almost constant. We suggest that the initially strong sublimation for the lowermost transition area is related to high temperature gradients at the lower boundary at the start of the simulation. Sublimation then reduces due to effects from consolidation. The effects of settling are very small in the vicinity of the uppermost layer transition, and the density difference in the layers is only 25 kg m^{−3}, which explains the constant and intermediate sublimation rate.
In this paper, we described in detail a hybrid Eulerian–Lagrangian computational approach to model the evolution of a dry snowpack. The model accounts for transport of heat and vapor, phase changes (sublimation and deposition), and mechanical settling processes. The ice mass balance is explicitly accounted for in the model formulation. It captures the evolution of the ice volume fraction in response to settling and phase changes. It constitutes an advectiondominated partial differential equation of hyperbolic type and is therefore conveniently solved with the method of characteristics, a popular Lagrangiantype scheme for such processes. Here, Lagrangian refers to the fact that the scheme tracks the motion of small reference volumes within the snow column by adjusting the node positions while at the same time accounting for phase changes within the moving snow. Solving the ice mass balance requires us to specify the vertical velocity as well as the mass production rate (sublimation rate/deposition rate). A closure for the velocity is derived by combining a common mechanical stress–strain relation with Glen's law and numerically approximating the resulting integrals. The deposition rate is due to vapor transport through a varying temperature field and can be determined from a diffusivetype process model that accounts for simultaneous heat and vapor transport. Due to its diffusive type (parabolic), a fixedgrid approach is most appropriate, referred to as an Eulerian approach. More specifically, we solved coupled heat and vapor transport by means of a firstorder implicitintime finitedifference approximation. The Eulerian scheme for the process model's diffusive part complies with the nonuniform mesh that results from the Lagrangian scheme for the ice volume fraction evolution. In order to solve the complete drysnow process model for the coupled evolution of the ice volume fraction, temperature field, vapor field and settling processes, the Eulerian and Lagrangian parts are iteratively applied following a straightforward operator split approach.
We have implemented our proposed numerical scheme as a series of sequential updates within one simulation time step. The implementation follows a modular, extendable approach, such that each process building block can easily be considered or neglected for verification or validation purposes. We applied our numerical core to conduct a series of simulation scenarios comprising isolated processes (pure settling, pure heat transport), twoprocess coupling scenarios (heat transport in the presence of settling, coupled heat and vapor transport) and fully coupled processes (heat and vapor transport in the presence of settling). We furthermore investigated different viscosity closures as well as a linear and a nonlinear version of Glen's law. A twolayer snowpack, consisting of a lower layer of higher density and an upper layer of lower density, served as a test case to demonstrate the feasibility of our approach. We simulated fields and profiles for the temperature, deposition rate, ice volume fraction and vertical velocity with a high spatial (∼ mm to cm) and temporal (∼ s to min) resolution.
We showed that our model implementation facilitates the comparison of various parametrizations and processes. This is enabled by the following:

Our Eulerian–Lagrangian scheme along with its vectorized implementation is flexible and extendable. Alternative model closures, e.g., for the viscosity and the vertical velocity, can easily be integrated. To close for the velocity, we have successfully tested a nonlinear strain rate closure commonly used in firn models (Lundin et al., 2017). The Lagrangian part of the solver (that accounts for the evolution of the ice volume fraction) can be singled out and coupled to an alternative process model, e.g., when accounting for firn conditions instead of dry snow.

The combination of an implicit Eulerian routine for the diffusiondominated operators (that controls the time stepping) and a Lagrangian routine for the advectiondominated operators ran stably and robustly for all considered simulation cases (different viscosity closures, different versions of Glen's law).

The numerical scheme allows for a high spatial resolution that resolves processes on the sublayer level. By construction, it relies on a mechanically consistent vertical velocity. This improves the accuracy since it makes the ad hoc specification of an artificial stress value for the uppermost layer (e.g., Vionnet et al., 2012), as required for conventional layerbased schemes, obsolete.

The modular setup of the software allows a systematic study of various model formulations, in which we selectively considered different combinations of process building blocks without finetuning the stability of the solver. This is important to enable empirical–numerical investigations of the relevance of different process couplings.

The incorporation of the higherorder mesh errors into the vapor and heat transport equations that account for deviations due to the nonuniform mesh increases accuracy especially in areas of large variations in node distance and of high temperature gradients.
In this paper, we mostly relied on numerical approximations that are of either first order (time integration and operator splitting) or second order (diffusion operator); the corresponding numerical solvers can be extended without conceptual difficulty, e.g., changing from a firstorder time integration to a higherorderintime integrator. We commented on this in the relevant section (Sect. 3). Our simulation consistently showed that vapor transport and phase change in the presence of strong temperature gradients can induce a stronger phase change activity and in particular a localized sublimation rate peak above the transition area between two layers. We furthermore showed that this has the potential to result in a localized ice volume fraction reduction above the transition area. This in itself is not new, as a similar behavior has been deduced in Hansen and Foslien (2015) and analyzed in detail in our companion paper (Schürholt et al., 2021). In addition to the existing results, we have shown that the increased phase change activity persists in the presence of settling (even more pronounced), for a constantviscosity closure in combination with a linear as well as a nonlinear version of Glen's law.
In our paper, we deliberately focused on discussing modularity and extendability in the context of snowpack modeling, e.g., by assessing a whole process cascade for one relatively simplified physical setting. In order to discuss these aspects in depth, we restricted ourselves to one relatively simple physical setting. We are well aware that as of today, our proposed numerical approach is not ready for operational use, and that was not our intention. At this point in time, we rather would like to contribute to the discussion on how future snowpack modeling can benefit from a consistently formulated, hybrid Eulerian–Lagrangian solver. Nevertheless, it is important to discuss whether the suggested scheme is amenable to further extensions required for an operational snowpack model.
Most importantly, the proposed scheme would need a generalization for surface mass gain (precipitation) or losses (sublimation). This bears two technical challenges. First, concurrent settling and precipitation result in a nonmonotonic vertical motion of the snowpack's upper surface, for which several techniques have been proposed in the past, e.g., based on a regularization approach (Wingham, 2000) or via kinetic boundary conditions as applied to sedimentation on ocean floors in Audet and Fowler (1992). A straightforward approach based on appending the computational grid sequentially during precipitation events likewise seems computationally feasible. A second challenge associated with the incorporation of precipitation events is the question of how to initialize the complete state (temperature, vapor and deposition rate) in the new snow layers. The latter is less critical in conventional layerbased schemes, as the necessary information reduces to “one value per layer”. While the first challenge mostly consists in overcoming technical subtleties in the actual implementation, the second requires a thoughtful formulation of physically consistent boundary conditions. Neither of the two challenges seems to pose a severe risk.
Another important addition to our proposed snowpack model is the presence of liquid water in the snow. Conceptually, similar modeling approaches could be used to derive a model for wet snow. While including potential phase changes from melting and freezing could be straightforwardly implemented via the source term c, it is the advective transport of liquid water that is more demanding. Liquid water transport is commonly modeled via the Richards equation (Wever et al., 2014) which would benefit from existing hybrid Eulerian–Lagrangian solution strategies, as shown for saturated media without mechanical settling (Huang et al., 1994). Furthermore, a model for wet snow requires a second energy balance to account for the liquid water temperature. Once set up, it can be integrated into our computational workflow (Fig. 2).
Finally, operational models generally include the capability to account for topological change within the snow column, to capture either layer coalescence if two initially separated snow layers merge into one or layer separation if an initially homogeneous layer splits into two. By construction, our computational approach does not require a dedicated treatment for layer coalescence or separation. Both are implicitly accounted for in the continuous description of stratigraphy as long as the feature that is to be resolved is larger than the chosen spatial resolution of the computational grid. In the current version of our model we use a dynamic time step adaptation based on the mesh Fourier number, which leads to a decreasing time step size from 40 s at initiation to 3 s after 48 h. Layer coalescence, e.g., after certain time intervals or whenever a specific time step size is undercut, could facilitate longer simulations runs. Otherwise the resolution can be increased to avoid layer splitups for regions with high gradients. Furthermore, our approach prevents the complete degeneration of layers as the ice volume fraction is constrained by the snow's maximum apparent density per construction of the scheme. Yet, while the theory suggests that layer coalescence and separation are not problematic, there might still be troublesome realistic test cases, especially when thinking about long simulation times. In order to address these and verify robustness, a series of benchmark tests have to be conducted. If necessary, the Eulerian–Lagrangian scheme in its current version can be equipped with occasional remeshing (along with a resampling of field variables) triggered by the degeneration of welldefined mesh quality criteria.
We believe that a flexible and extendable computational approach, such as the one described in this paper, will be key for future snowpack modeling to facilitate the use of standardized benchmark problems (potentially used during a continuous integration) and allow us to systematically assess the model's predictive power, including uncertainty quantification, parameter estimation and model selection.
A1 Vapor saturation density
An empirical expression for the vapor saturation density ${\mathit{\rho}}_{\mathrm{v}}^{\text{eq}}\left(T\right)$ in terms of temperature T is formulated based on the empirical formulation for vapor saturation pressure from Libbrecht (1999) and reads
with coefficients a_{0} = 3.6636 × 10^{12} $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{2}}$, a_{1} = −1.3086 × 10^{8} $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}$, a_{2} = −3.3793 × 10^{6} $\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{2}}$, f = 461.31 $\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{\mathrm{1}}$, T_{m} = 273.15 K and T_{ref} = 6150 K. f is the specific gas constant for water vapor. Note that division by fT accounts for the conversion from pressure (Pa) (as used in Part 1) to density (kg m^{−3}).
A2 Model parameters in the transport and phase change equations
The effective vapor mass diffusion coefficient D_{eff}(ϕ_{i}) in terms of the ice volume fraction ϕ_{i} is taken from Calonne et al. (2014) but is extended by the heaviside function Θ to hinder vapor diffusion for ice volumes above twothirds:
with D_{0} = 2.036 × 10^{−5} m^{2} s^{−1} being the vapor diffusion constant in air.
The effective thermal conductivity k_{eff}(ϕ_{i}) in terms of the ice volume fraction ϕ_{i} is taken from Calonne et al. (2011) and reads
with coefficients a_{0} = 0.024, a_{1} = −1.23 × 10^{−4} and a_{2} = 2.5 × 10^{−6} and ice density ρ_{i}.
The effective heat capacity (ρC)_{eff}(ϕ_{i}) in terms of the ice volume fraction ϕ_{i} is taken from Calonne et al. (2014) and Hansen and Foslien (2015) and reads
with C_{i} being ice heat capacity, C_{a} air heat capacity, ρ_{i} ice density and ρ_{a} air density.
A3 Constant viscosity for the twolayer case
A3.1 Linear Glen's law, ${\mathit{\eta}}_{\text{const},m=\mathrm{1}}$
We derived intermediate ice volume fraction ϕ_{i,const} = 0.1125 and temperature T_{const} = 263 K values from the initial condition of the twolayer case and insert them as constants into Eq. (4).
A3.2 Nonlinear Glen's law ${\mathit{\eta}}_{\text{const},m=\mathrm{3}}$
Equation (4) does not hold for the Glen exponent m=3; therefore we derive an adjusted constant viscosity ${\mathit{\eta}}_{\text{const},m=\mathrm{3}}$ via the constitutive equation (Eq. 3)
with ${\dot{\mathit{\u03f5}}}_{\text{lit}}\equiv $ 10^{−6} s^{−1} being a strain rate value from the literature (Johnson, 2011) and ${\mathit{\sigma}}_{max}\equiv $ 547.71 Pa the maximum stress value obtained from the initial snow density profile of the twolayer case. Equation (A5) is then solved for the constant viscosity ${\mathit{\eta}}_{\text{const},m=\mathrm{3}}$.
A3.3 Restrict infinite ice volume growth
To hinder infinite ice volume growth, the constant viscosity η_{const,m} is combined with a power law that yields exponential growth of viscosity for cells with ϕ_{i}>0.95:
with pl1=690 and pl2=650. The constant viscosity is then multiplied with the power law (η_{const,m}PL(ϕ_{i})) so that computational nodes with ϕ_{i}>0.95 are assigned and viscosity grows exponentially. Note that for better readability the multiplication with the power law is omitted in the equations of this paper.
For the temperature equation (Eq. 22) the higherorder mesh error is
and for the vapor transport equation (Eq. 23) it is
For the matrix equations Eqs. (25) and (24) the higherorder mesh errors are defined as E_{T} and E_{c}.
Matrix A is defined as follows:
For the heat equation (Eq. 24) (A_{T}) the entries are
and for the vapor transport (Eq. 25) (A_{c}) the entries are
with
Matrix B is defined as follows:
For Eq. (24) (B_{T}) the entries are
and for Eq. (25) (B_{c}) they are
Matrix E is defined as follows:
consisting of the terms for the heat equation (E_{T}) (Eq. 24)
and for the vapor equation (E_{c}) (Eq. 25)
D1 Nonlinear Glen's law
The model code is available via https://doi.org/10.5281/zenodo.5588308 (Simson and Kowalski, 2021).
The study concept, underlying model and methodology were devised by AS, HL and JK. Model analysis, software implementation and simulation runs were performed by AS supervised by JK. Test case analysis and discussion, data visualization, and manuscript preparation were carried out by AS with contributions from JK and HL.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We thank the two anonymous reviewers for their helpful suggestions on the manuscript. The authors were supported by the Helmholtz School for Data Science in Life, Earth and Energy (HDSLEE). The work was furthermore supported by the Federal Ministry of Economic Affairs and Energy, on the basis of a decision by the German Bundestag (50 NA 1908).
This paper was edited by Mark Flanner and reviewed by two anonymous referees.
Audet, D. and Fowler, A.: A mathematical model for compaction in sedimentary basins, Geophys. J. Int., 110, 577–590, https://doi.org/10.1111/j.1365246x.1992.tb02093.x, 1992. a, b
Bader, H.P. and Weilenmann, P.: Modeling temperature distribution, energy and mass flow in a (phasechanging) snowpack. I. Model and case studies, Cold Reg. Sci. Technol., 20, 157–181, https://doi.org/10.1016/0165232x(92)90015m, 1992. a, b, c, d, e
Bartelt, P. and Christen, M.: A computational procedure for instationary temperaturedependent snow creep, Springer Berlin Heidelberg, https://doi.org/10.1007/BFb0104195, 2007. a, b
Bartelt, P. and Lehning, M.: A physical SNOWPACK model for the Swiss avalanche warning: Part I: numerical model, Cold Reg. Sci. Technol., 35, 123–145, https://doi.org/10.1016/s0165232x(02)000745, 2002. a, b, c, d, e, f, g, h, i, j
Brun, E., Martin, E., Simon, V., Gendre, C., and Coleou, C.: An Energy and Mass Model of Snow Cover Suitable for Operational Avalanche Forecasting, J. Glaciol., 35, 333–342, https://doi.org/10.3189/s0022143000009254, 1989. a, b, c
Brun, E., David, P., Sudul, M., and Brunot, G.: A numerical model to simulate snowcover stratigraphy for operational avalanche forecasting, J. Glaciol., 38, 13–22, https://doi.org/10.1017/s0022143000009552, 1992. a
Calonne, N., Flin, F., Morin, S., Lesaffre, B., du Roscoat, S. R., and Geindreau, C.: Numerical and experimental investigations of the effective thermal conductivity of snow, Geophys. Res. Lett., 38, L23501, https://doi.org/10.1029/2011GL049234, 2011. a
Calonne, N., Geindreau, C., and Flin, F.: Macroscopic Modeling for Heat and Water Vapor Transfer in Dry Snow by Homogenization, J. Phys. Chem. B, 118, 13393–13403, https://doi.org/10.1021/jp5052535, 2014. a, b, c, d, e, f, g, h, i, j
Domine, F., Barrere, M., and Sarrazin, D.: Seasonal evolution of the effective thermal conductivity of the snow and the soil in high Arctic herb tundra at Bylot Island, Canada, The Cryosphere, 10, 2573–2588, https://doi.org/10.5194/tc1025732016, 2016. a
Domine, F., Picard, G., Morin, S., Barrere, M., Madore, J.B., and Langlois, A.: Major Issues in Simulating Some Arctic Snowpack Properties Using Current Detailed Snow Physics Models: Consequences for the Thermal Regime and Water Budget of Permafrost, J. Adv. Model. Earth Sy., 11, 34–44, https://doi.org/10.1029/2018ms001445, 2019. a, b
Farlow, S. J.: Partial Differential Equations for Scientists and Engineers, Dover Publications, Mineola, New York, USA, reprint of the John Wiley & Sons, New York, USA, 1982 edition, 1993. a
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/tc918572015, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m
Helfricht, K., Hartl, L., Koch, R., Marty, C., and Olefs, M.: Obtaining subdaily new snow density from automated measurements in high mountain regions, Hydrol. Earth Syst. Sci., 22, 2655–2668, https://doi.org/10.5194/hess2226552018, 2018. a
Huang, K., Zhang, R., and van Genuchten, M. T.: An EulerianLagrangian approach with an adaptively corrected method of characteristics to simulate variably saturated water flow, Water Resour. Res., 30, 499–507, https://doi.org/10.1029/93WR02881, 1994. a
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, Front. Earth Sci., 8, 249, https://doi.org/10.3389/feart.2020.00249, 2020. a, b, c, d, e, f
Johnson, J. B.: Snow Deformation, in: Encyclopedia of Snow, Ice and Glaciers, 2011 edition, edited by: Singh, V. P., Singh, P., and Haritashya, UK, Springer, Dordrecht, the Netherlands, 1041–1045, https://doi.org/10.1007/9789048126422_501, 2011. a, b
Judson, A. and Doesken, N.: Density of Freshly Fallen Snow in the Central Rocky Mountains, B. Am. Meteorol. Soc., 81, 1577–1587, https://doi.org/10.1175/15200477(2000)081<1577:DOFFSI>2.3.CO;2, 2000. a
Kirchner, H. K., Michot, G., Narita, H., and Suzuki, T.: Snow as a foam of ice: Plasticity, fracture and the brittletoductile transition, Philos. Mag. A, 81, 2161–2181, https://doi.org/10.1080/01418610108217141, 2001. a
Kowalski, J. and Torrilhon, M.: Moment Approximations and Model Cascades for Shallow Flow, Commun. Comput. Phys., 25, 669–702, https://doi.org/10.4208/cicp.oa20170263, 2019. a
Krinner, G., Derksen, C., Essery, R., Flanner, M., Hagemann, S., Clark, M., Hall, A., Rott, H., BrutelVuilmet, 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.: ESMSnowMIP: assessing snow models and quantifying snowrelated climate feedbacks, Geosci. Model Dev., 11, 5027–5049, https://doi.org/10.5194/gmd1150272018, 2018. a, b
Krol, Q. and Löwe, H.: Analysis of local ice crystal growth in snow, J. Glaciol., 62, 378–390, https://doi.org/10.1017/jog.2016.32, 2016. a
Lacroix, M. and Garon, A.: Numerical Solution of phase change problems: an EulerianLagrangian approach, Numer. Heat Tr. BFund., 19, 57–78, https://doi.org/10.1080/10407799208944922, 1992. a
Lehning, M., Bartelt, P., Brown, B., Fierz, C., and Satyawali, P.: A physical SNOWPACK model for the Swiss avalanche warning: Part II. Snow microstructure, Cold Reg. Sci. Technol., 35, 147–167, https://doi.org/10.1016/s0165232x(02)000733, 2002. a, b
LeVeque, R. J.: Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, https://doi.org/10.1017/CBO9780511791253, 2002. a
Libbrecht, K. G.: Physical properties of ice, available at: http://www.cco.caltech.edu/~atomic/snowcrystals/ice/ice.htm (last access: 29 July 2021), 1999. a, b
Lundin, J. M., Stevens, C. M., Arthern, R., Buizert, C., Orsi, A., Ligtenberg, S. R., Simonsen, S. B., Cummings, E., Essery, R., Leahy, W., Harris, P., Helsen, M. M., and Waddington, E. D.: Firn Model Intercomparison Experiment (FirnMICE), J. Glaciol., 63, 401–422, https://doi.org/10.1017/jog.2016.114, 2017. a, b, c
Mellor, G. L. and Blumberg, A. F.: Modeling Vertical and Horizontal Diffusivities with the Sigma Coordinate System, Mon. Weather Rev., 113, 1379–1383, https://doi.org/10.1175/15200493(1985)113<1379:mvahdw>2.0.co;2, 1985. a
Meyer, C. R., Keegan, K. M., Baker, I., and Hawley, R. L.: A model for Frenchpress experiments of dry snow compaction, The Cryosphere, 14, 1449–1458, https://doi.org/10.5194/tc1414492020, 2020. a
Morland, L., Kelly, R., and Morris, E.: A mixture theory for a phasechanging snowpack, Cold Reg. Sci. Technol., 17, 271–285, https://doi.org/10.1016/s0165232x(05)800060, 1990. a
Morland, L. W.: A fixed domain method for diffusion with a moving boundary, J. Eng. Math., 16, 259–269, https://doi.org/10.1007/bf00042720, 1982. a
Paterson, W.: The Physics of Glaciers, third edition, Elsevier, Oxford, England, https://doi.org/10.1016/C2009014802X, 1994. a
Schleef, S., Löwe, H., and Schneebeli, M.: Influence of stress, temperature and crystal morphology on isothermal densification and specific surface area decrease of new snow, The Cryosphere, 8, 1825–1838, https://doi.org/10.5194/tc818252014, 2014. a
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 Discuss. [preprint], https://doi.org/10.5194/tc202172, in review, 2021. a, b, c, d, e, f, g, h
Schweizer, J., Jamieson, J. B., and Schneebeli, M.: Snow avalanche formation, Rev. Geophys., 41, 1016, https://doi.org/10.1029/2002rg000123, 2003. a
Simson, A. and Kowalski, J.: Eulerian_Lagrangian_snow_solver: final paper submission TC, Zenodo [code], https://doi.org/10.5281/zenodo.5588308, 2021. a
ViallonGalinier, L., Hagenmuller, P., and Lafaysse, M.: Forcing and evaluating detailed snow cover models with stratigraphy observations, Cold Reg. Sci. Technol., 180, 103163, https://doi.org/10.1016/j.coldregions.2020.103163, 2020. a
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. Model Dev., 5, 773–791, https://doi.org/10.5194/gmd57732012, 2012. a, b, c, d, e, f, g, h, i, j, k, l
Wever, N., Fierz, C., Mitterer, C., Hirashima, H., and Lehning, M.: Solving Richards Equation for snow improves snowpack meltwater runoff estimations in detailed multilayer snowpack model, The Cryosphere, 8, 257–274, https://doi.org/10.5194/tc82572014, 2014. a
Wever, N., Valero, C. V., and Fierz, C.: Assessing wet snow avalanche activity using detailed physics based snowpack simulations, Geophys. Res. Lett., 43, 5732–5740, https://doi.org/10.1002/2016GL068428, 2016. a
Wiese, M. and Schneebeli, M.: Earlystage interaction between settlement and temperaturegradient metamorphism, J. Glaciol., 63, 652–662, https://doi.org/10.1017/jog.2017.31, 2017. a
Wingham, D. J.: Small fluctuations in the density and thickness of a dry firn column, J. Glaciol., 46, 399–411, https://doi.org/10.3189/172756500781833089, 2000. a, b
 Abstract
 Introduction
 Physical model
 Computational approach
 Results and discussion
 Summary and conclusions
 Future work and challenges
 Appendix A: Formulas required for the process model
 Appendix B: Higherorder mesh errors to correct for nonuniform mesh
 Appendix C: Matrices from temperature and vapor transport equations
 Appendix D: Additional figures
 Code availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Review statement
 References
 Abstract
 Introduction
 Physical model
 Computational approach
 Results and discussion
 Summary and conclusions
 Future work and challenges
 Appendix A: Formulas required for the process model
 Appendix B: Higherorder mesh errors to correct for nonuniform mesh
 Appendix C: Matrices from temperature and vapor transport equations
 Appendix D: Additional figures
 Code availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Review statement
 References