Comment on tc-2021-72 Anonymous Referee # 1 Referee comment on " Elements of future snowpack modeling – part 1 : A physical instability arising from the non-linear coupling of transport and phase changes

Schurholt et al. show that coupled equations for heat transport, vapour diffusion and ice mass conservation in snow permit wave solutions in density. The linear stability analysis is nice work that, together with numerical solutions of the nonlinear equations, demonstrates that these are true mathematical solutions and not numerical artefacts. The setting is limited to be somewhat short of a full snow thermodynamics model, and the question of how mm-scale waves in solutions of the continuous equations relate to a bicontinuous material with mm-scale structure remains open.


Introduction
Neglecting vapor transport in the overall mass balance of a snowpack is considered as a serious uncertainty in snow modeling.
As hypothesized in recent work on shallow tundra snowpacks (Barrere et al., 2017;Domine et al., 2016) persistent temperature gradients throughout the season may contribute to the depletion of snow density at the bottom of the snowpack due to persistent upward vapor fluxes. This problem is relevant for applications e.g. in permafrost where temperature gradients may induce the 20 drying of soils with a feedback on snow metamorphism in the adjacent bottom layer (Domine et al., 2016), in turn affecting the ground thermal regime. Other applications of vapor transport comprise post-depositional re-distribution of stable water monly avoided by using the Lagrangian frame of reference of the settling equation (Brun et al., 1989;Lehning et al., 2002).
This complicates the assessment of phase change models in their originally published Eulerian form. Stand-alone numerical investigations of specific model components are therefore justified and necessary to understand non-linear effects and numerical 60 requirements for the desgin of efficient solvers in future snowpack models.
It is the aim of the present paper to advance the understanding of coupled heat and mass transport in snow by a careful numerical analysis of existing homogenization schemes. To overcome the limited flexibility in existing snowpack models we have implemented a stand-alone solver for the PDEs using the Finite Element (FE) framework FEniCS (Alnaes et al., 2015). The python-based computing platform was previously used for other problems in cryospheric sciences (Cummings, 2016). We focus 65 on the non-linear feedback of heat and vapor transport on an evolving ice phase. We consider the simplest setting and neglect convection in the gas phase and settling in the solid phase to provide a sound reference for future extensions. The question how these effects can be coupled to settling will be addressed in a companion paper (Simson et al., 2021). Our approach will reveal the rich mathematical complexity that is hidden in published models and will provide evidence that the origin of this complexity is the density dependence of the effective transport coefficients in the presence of phase changes. We will show that 70 this coupling causes the formation of wave patterns in the ice volume fraction profile as a true mathematical feature of the nonlinear PDE system. This is confirmed by an analytical, linear stability analysis which relates unstable behavior to the density dependence of the effective (heat and mass) diffusion constants. The results suggest that previously observed oscillations in the numerical treatment (Adams and Brown, 1990;Jafari et al., 2020) were not numerical problems but may rather have signalled emergent physics. With this work we seek to contribute to an understanding if and how these features should be taken into 75 account in future work.
The paper is organized as follows. In Sec. 2 we state the governing partial differential equations from the homogenization schemes (Calonne et al., 2014) and (Hansen and Foslien, 2015). In Sec. 3 we outline the Finite Element solution of the weak formulation of the problem and its implementation in FEniCS. In Sec. 4 we first provide an inter-comparison of the two models in their original test scenarios and the scenario of a thin "Gaussian crust" as a smooth density heterogeneity. In Sec. 5 80 we characterize the observed migration of the ice phase under vapor re-crystallization by comparing the full model with an approximate advection equation for the ice phase which is derived under simplified assumptions. In Sec. 6 we detail the wave patterns that emerge in the numerical solution of the Gaussian crust. The analysis of the mesh-resolution, integration time step and the residuals indicates that the wave patterns are intrinsic features of the mathematical model and not numerical artefacts. This is confirmed by the linear stability analysis in Sec. 6.2. After a few sensitivity tests we will discuss practical consequences 85 of our study for vapor transport modeling in snow models (Sec. 8) and provide summarizing conclusions.

Homogenized heat and mass transport
As a theoretical starting point we focus on two recently published homogenized formulations for an evolving vapor phase, namely (Calonne et al., 2014) and (Hansen and Foslien, 2015). Despite differences in their theoretical homogenization ap-proach, asymptotic expansion (Calonne et al., 2014) and mixture theory (Hansen and Foslien, 2015), both models have a very similar mathematical structure that resembles earlier work of (Bader and Weilenmann, 1992).
2.1 Vapor scheme from (Calonne et al., 2014) The two-scale expansion for (vapor) mass and energy (Calonne et al., 2014) leads to the following set of equations for the vapor density ρ v and the temperature T Here φ i is the ice volume fraction, D eff the effective diffusion coefficient, ρ i is the density of ice density (assumed to be constant), s the ice-air interface area per unit volume, v n the volume averaged interface surface normal velocity indicating ice growth or decay, (ρC) eff denotes the effective volumetric heat capacity (times snow density), k eff the effective thermal conductivity and L the latent heat of sublimation.

100
The surface area density s, reflecting the current state of the microstructure and principally evolving in time, is assumed to be constant here in accordance with (Calonne et al., 2014).
The source term −ρ i sv n on the right hand side (r.h.s.) of the vapor equation (1) quantifies phase changes, i.e. the net condensation rate. This implies the form of the energy source term L sv n through latent heat. In the simplest setting (Calonne et al., 2014), the volume averaged interface velocity is given by where β is the inverse growth velocity and ρ eq v (T ) the saturation vapor density. Note that v n strongly depends on the temperature, indicating the temperature feedback on the vapor mass balance. In order to close the system (1)-(3) parametrizations for the effective PDE coefficients must be provided (see below in section 2.3). The feedback of a spatio-temporally evolving ice phase φ i (z, t) has not been considered in (Calonne et al., 2014). 110 2.2 Vapor scheme from (Hansen and Foslien, 2015) A similar coupling scheme of the vapor transport to the energy equation has been put forward in (Hansen and Foslien, 2015).
Using the same notation as above, these equations can be written in the form Here the r.h.s. of both equations is expressed in terms of the condensation rate c.
The comparison of both models reveals that the vapor mass balance of (Hansen and Foslien, 2015) (4) can be obtained from the corresponding vapor mass balance of (Calonne et al., 2014) (1) using the following reasoning: The deviation of the vapor concentration ρ v from its equilibrium value ρ eq v (T ) will be mostly small. By assuming that the deviation ρ v − ρ eq v (T ) diffuses faster than it relaxes locally, ρ v = ρ eq v (T ) is used in the diffusion term in (1), while a deviation from equilibrium is maintained 120 in the reaction term. Splitting the reaction term into an equilibrium contribution ∂ t ρ eq v and a perturbation, and using the known temperature dependence of the saturation vapor density finally yields (4). (Hansen and Foslien, 2015) supports this derivation with a (non-rigorous) scaling argument to avoid two conflicting equations for the temperature T .
Though both presented models are strikingly similar in structure, it should be emphasized that while system (1)-(3) is closed and can be solved for ρ v and T , the same is not true for system (4) -(5) if we were to use the same condensation rate closure (3).

125
Due to the simplifying assumption outlined above, Hansen's model lacks an evolution equation for the (perturbed) vapor mass balance. Since both, vapor mass balance (4), and energy balance (5) are actually formulated in terms of the temperature T , (Hansen and Foslien, 2015) proceeds by consolidating both equations into a single one that no longer contains the condensation The latter now poses a non-linear, yet closed PDE that can be solved for the evolving temperature. Back-substituting the resulting temperature into either (4) or (5) allows to evaluate the condensation rate c, without remaining freedom to impose yet another closure relation. An additional, important implication of the near-equilibrium assumption in (Hansen and Foslien, 2015) is the impossibility to impose arbitrary boundary conditions for the vapor phase. By construction of the model, the boundary conditions are Dirichlet values for the equilibrium vapor concentration consistent with the boundary temperatures.

135
Similar to (Calonne et al., 2014) the feedback of a spatio-temporally evolving ice phase φ i (z, t) was not been considered in (Hansen and Foslien, 2015).

Parametrization of the PDE coefficients
Despite similarities in the forms of the PDEs, both models have used different parametrizations for the transport coefficients D eff , k eff and (ρC) eff and ρ eq v . While the two-scale homogenization presented in (Calonne et al., 2014) contains derived 140 expressions for the effective properties in terms of the microstructure as a byproduct, mixture theory in (Hansen and Foslien, 2015) relies on independently postulated parametrizations. The implications of these differences are presently still under debate (Fourteau et al., 2020).
In general, there is a broad agreement that all effective parameters are primarily influenced by the density or ice volume fraction φ i . The goal of the present paper is not an exhaustive inter-comparison of different formulations of the parametrized 145 coefficients in the analysed snow pack models. We rather want to focus on the intrinsic features and physical implications of the underlying process models. We therefore aim to keep differences due to specific flavours of the coefficients at the minimal level. To this end we use the same parametrization for the equilibrium vapor pressure and the effective heat capacity in both models but use different parametrizations for the effective thermal conductivity/diffusion constant. In the following we refer to these PDE parametrizations as Calonne parametrization and Hansen parametrization. All coefficients are stated explicitly in 150 appendix A.

Feedback from an evolving ice phase
The consistent treatment of phase changes requires a dynamic ice phase that evolves through recrystallization alongside with the vapor phase in a mass-conserving way. This was neither considered in (Calonne et al., 2014) nor in (Hansen and Foslien, 2015). To investigate the feedback of an evolving ice phase on the two models from above we supply (1)-(3) and (4)-(5) with a 155 dynamic ice mass conservation equation (Bader and Weilenmann, 1992;Krol and Löwe, 2018). In the absence of settling, but presence of phase changes the continuity equation reduces to an ordinary differential equation for each location in space to balance the source terms in the vapor equations above. In Calonne's model v n can be evaluated from the closure relation (3), whereas Hansen's c has to be reconstructed from either (4) or (5) as outlined before.

160
In summary, all symbols and parameter values used in this study are provided in Table 1.

Finite Element solution in FEniCS
To minimize the coding overhead and focus on the physical problem while keeping access to advanced numerical adjustments, the coupled PDE model was implemented using the python-based Finite Element framework FEniCS (Alnaes et al., 2015).
It provides a high-level API for the solution of PDEs in their weak formulation, with capabilities for parallelization, and 165 cross-platform portability via docker images. FEniCS' flexibility arises from its interface to the Unified Form Language (UFL) (Alnaes et al., 2015) which supports an intuitive declaration of discrete FE formulations of variational problems.
Below we outline the spatial and temporal discretization for the weak formulation of both systems (1)-(3) and (4)-(5), and describe the solution strategy employed for the coupled PDEs.

170
For the spatial discretization we note that the non-linear PDE systems, of interest can be formally rewritten in the form where F (u) is a nonlinear differential operator and u is the (vector valued) solution that comprises the unknown fields. Following a standard Galerkin approach, the PDE's weak formulation is obtained by multiplying F (u) = 0 with test functions v and integrating over the domain. The Galerkin method then approximates the solution as the sum over a discrete set of basis func-175 tions u i , viz u ≈ u h = i u i v i . This procedure allows to apply the differential operator to the test functions. The test functions have a small support only, which simplifies the integral evaluations significantly. By choosing different test functions v = v i , the PDE system is reduced to finding the roots of a (non-linear) algebraic system. This standard Galerkin procedure is e.g detailed in (Donea and Huerta, 2003). The implementation into the open source finite element software FEniCS is straightforward. Time derivatives are approximated using the theta method (also known as Rothe's method) (Donea and Huerta, 2003). For a PDE of the form where u n denotes the state vector's solution at time t n . Note that for θ = 0 and θ = 1 this reduces to standard first order explicit and implicit Euler methods, respectively, whereas for θ ∈ (0, 1) is constitutes a weighted average of both.

Coupling scheme
Due to the different time scales of the involved equations, a monolithically coupled solution for the vector u = (T, ρ v , φ i ) would be most consistent, yet turns out to be inefficient. The vapor equation has by far the fastest dynamics, followed by the energy.

190
The ice mass balance instead has a much slower dynamics. Since both water vapor and energy equation We apply different theta values for each differential operator. For the diffusion operators in the vapor and energy equations we use θ = 0.5 (Crank-Nicolson) which is known to be stable and converge of second order for linear operators. As the phase 200 change induced source terms are expected to be stiff, a fully implicit, unconditionally stable Backwards Euler scheme with θ = 1.0 is applied.

Model comparison
To evaluate the models, we tested them on three different scenarios comprising specific combinations of initial conditions (IC) and boundary conditions (BC). The first two scenarios discussed in sections 4.1 and 4.2 are taken from (Calonne et al.,205 2014) and (Hansen and Foslien, 2015), respectively, while the last (section 4.3) is the simulation of a Gaussian shaped crust.
The three scenarios test different relevant aspects of snow modeling, namely the transient response to time dependent BCs, piecewise-linear ice profiles covering the entire range of snow densities and lastly a high density layer in a small sample with smooth, high density gradients. For each of the three physical scenarios, always evaluate three model formulations: First, Calonne's equations (1), (2), (7) with Calonne coefficient closure from Appendix A (referred to Cal-Eq/Cal-Par). Second, 210 Hansen's equations (6), (7) with the Hansen coefficient closure from Appendix A (Han-Eq/Han-Par). Third, the mixed case of Hansen's equations (A5), (7) with Calonne coefficient closure (Han-Eq/Cal-Par).

Scenario 1: Homogeneous snow -transient heating at the boundary
The first scenario is taken from (Calonne et al., 2014) who investigated the response of a homogeneous snow layer to transient heating. To this end we use the IC 215 T (z, 0) = T 0 (10) with T 0 = 273K and φ i,0 = 0.3. We employ a a fixed Dirichlet BC at z = 0 and a transient temperature drop at z = H, viz where the transient drop is characterized by the Heaviside step function Θ with parameters τ = 5h and T H = 263K.
For this combination of IC and BC we obtain the results in Figure (1 which coincides with Han-Eq/Han-Par. Only a small difference can be observed close to the right boundary, where the changes in T and ρ v are the fastest and therefore the assumption ρ v = ρ eq v (T ) underlying (Hansen and Foslien, 2015) is violated.

Scenario 2: Heterogeneous snow -fixed temperature boundary conditions
Next we investigate the test case envisaged by (Hansen and Foslien, 2015) that comprises a piece-wise linear density profile as an approximation for a layered snowpack. The IC is given by with T 0 = 273K and T H = 253K. The BC are given by .    typically encountered in tomography experiments (Hammonds et al., 2015). In the following we solely focus on differences between (Cal-Eq/Cal-Par) and (Han-Eq/Han-Par) and investigate this difference at different physical times. 245 We refer to this scenario as a Gaussian crust. As IC we employ with T 0 = 273K, T H = 253K, φ i,0 = 0.3, ∆φ i,0 = 0.2, z 0 = 0.01m, σ 2 = 5 · 10 −7 m. For the Gaussian crust we use the same 250 boundary conditions as in (17).   First, the Gaussian crust shows a quasi-advection towards the warm boundary, despite the absence of an explicit advection 255 term in the ice equation. During this quasi-advection, density gradients steepen on the cold side and flatten on the warm side.
Far away from the crust a linear density gradient emerges in the domain as a consequence of the boundary conditions. The difference between the models lies in the apparent advection velocity. This is consistent with the observed differences in the phase changes since recrystallization rate differ approximately by a factor of two.
considerable later time with smaller magnitude. For the Hansen model these are apparent immediately.
These two observations are further detailed below.

Quasi-advection of the ice phase
The quasi-advection of density heterogeneities in the ice phase despite the absence of an explicit advection term in (7) is a consequence of phase changes (vapor sublimation and re-sublimation) as an intrinsic feature of the coupled system.

265
To understand the origin of this quasi-advection, we approach the coupled system (1), (2) and (7) with several analytical approximations. Justified by the separation of time-scales we first restrict ourselves to stationary heat transfer and neglect the phase changes in the energy equation following arguments given in (Calonne et al., 2014). This is a considerable simplification since the two-way coupling between heat and vapor is eliminated. Consequently, the heat equation ∇k eff ∇T = 0 can be solved for the BC (17) independent of the vapor transport process even for an inhomogeneous ice profile, with its exact solution given The local temperature (and its gradient) is thus a functional of the heterogeneous, non-constant ice volume fraction φ i via the dependence of the effective conductivity k eff on φ i . We express this in the form 275 which defines the functional F T via (21). Motivated by the similarity between the solutions of (Calonne et al., 2014) and (Hansen and Foslien, 2015) found in the previous section, we adopt the simplification from (Hansen and Foslien, 2015) and account for deviations of the vapor concentration from equilibrium only in the phase change (source) term. Again, we use a stationary diffusion equation ∇D eff ∇ρ eq v = ρ i sv n due to the separation of time-scales between vapor and ice. This allows to rewrite the ice mass conservation in a conservative form as a single advection equation with a functional that expanded ρ eq v (T ) by means of the chain rule and substituted in the previously derived temperature expression (22) follow from the BC ρ = ρ eq v of the vapor phase. We stress the significance of this result. First, by using the approximations for the heat and vapor transfer, the three coupled heat and mass conservation equations can be equivalently cast into a single continuity equation for the ice phase. Second, the form of this continuity equation is a non-linear and non-local advection equation, which explains the nature of this quasi-290 advection as a variant of shock formation reminiscent of the non-linear Burgers equation (Du et al., 2012).
To test the derived approximation we have compared the numerical solution of (23) to the full solution of the Calonne model as shown in figure (4) for the Gaussian crust. In general, the agreement between the two models is very close, yielding almost identical results in homogeneous regions, i.e. towards the boundaries and close to the center of the crust. In the high gradient regions both models predict a steepening of the gradient, which is faster in the advection equation (21). Since (21)  Despite the remaining differences, the approximation (21) allows to trace back unambigously the quasi-advection and gradient steepening to the dependence of the effective diffusion constant and the effective thermal conductivity on ice volume fraction: If D eff was linear in φ i and k eff was constant, (21) would take the form of a simple advection equation ∂ t φ i + v∇φ i = 0 300 with constant v. This would imply a shape invariant advection of any initial density profile. The fact that D eff decreases, and k eff increases with ice volume fraction, decreases the ice flux functional G in high density regions over lower density regions.
Both effects contribute in the same direction and explain the emergent asymmetry in the crust. (3), the numerical solution of the advection equation also displays oscillations at the boundary, already after a few hours. This will be further detailed in the following.

Pattern Formation
All numerical solutions so far are subject to oscillations "at some point in time" irrespective of the details. Therefore we will investigate the nature of these oscillations, now only considering the full Calonne model applied to the Gaussian crust.

Oscillatory solutions
The following example shows simulations with varying mesh size (number of elements n e ) and varying time steps dt for the while (5)c shows close-up of the sublimating side of the crust. For very low mesh resolution (n e = 100, corresponding to 0.2mm) and large time step (1 min) the solution oscillates from node-to-node (a) while for sufficiently high resolution the numerical solution converges to a smooth wave pattern independent of temporal and spatial resolution. This is interesting, 315 as it suggests that these waves are true, intrinsic features of the full Calonne model equations, rather than an artefact of the numerical scheme. As an additional check we have analyzed the residuals of the numerical solution (cf. Appendix (B)) which confirms that patterns persist even when enforcing the residuals to approach zero.

Linear stability analysis
To comprehend the oscillatory nature of the solution we analyzed the problem theoretically within perturbation theory. Pattern 320 formation in non-linear PDE systems can be generally understood by investigating the dynamics of perturbations around a known stationary state via linear stability analysis of the dynamical equations. To do so, we start from the full, non-linear, only subject to the simplification that on the r.h.s. the division by ρ eq (T ) in Eq. (3) is subsumed in the prefactor α which is assumed to be constant.
at the bottom z = 0 and top z = H of the domain. To proceed analytically, we need to make an additional assumption to obtain an exact stationary solution of the system. To this end we will assume that the equilibrium vapor concentration is a 335 linear function in T . This is reasonable within small layers where temperature differences are small and a linearization of the equilibrium vapor curve can be always justified. This linearization is written in the form as an expansion around a reference temperature T ref (e.g. the mean temperature in the sample). With this assumption, a stationary solution of Eq.(26) can be obtained by inspection. It is easily verified that 340 satisfies (26) with BC (27) for arbitrary, constant volume fraction φ i,0 .
To carry out the stability analysis we use a vector notation and combine the fields in the vector u := (T, ρ v , φ i ). Then the 345 PDE system (28) can be written in matrix form with state dependent (3 × 3) matrices K(u), C(u) constant matrix R and constant vector f defined by Then, the (stable) stationary state of (29) is denoted by the vector

355
The non-linearities of (30) emerge from the dependence of the matrices K ≡ K(u) and C ≡ C(u) on u. Note that this dependence is solely through the third component of u (i.e. φ i ) which highlights the special role of the density (evolution).
Next we investigate the linear stability of the fixed point (36). To this end we make an ansatz and investigate the dynamics of the deviation from the stationary state u (0) by deriving an equation for u (1) through a pertur-360 bation expansion of (30) in first order. This procedure (carried out in detail in Appendix (C)) yields an equation for the Fourier modesũ (1) (k) of the perturbation for wave number k which can be written in the form with a wave number dependent matrix with constant matrix coefficients K (0) , V (0) , C (0) . The eigenvalues of M k control the behavior of the system in the vicinity of the stationary state. The matrix M k is non-Hermitian due to the "advection" matrix V (0) and the "phase change" matrix R. Thus these two processes can principally cause eigenvalues with positive real parts (causing a growth of perturbations) and non-zero imaginary parts (causing oscillatory behavior). In the case of V (0) = R = 0, Eq. (39) predicts negative eigenvalues for all k > 0 and no pattern formation as it should be for simple diffusion equations.

370
The matrix M k can be diagonalized in closed form using the symbolic algebra software MAPLE and eigenvalues subsequently separated into real and imaginary parts. Evaluating the eigenvalues with the corresponding parameter values from Table (1) shows that one of the eigenvalues has positive real part and non-zero imaginary part indicating a wave instability. This is controlled by the matrix V (0) with the coefficients k 1 and D 1 (cf. Eq. (C14) in appendix (C)) which characterize the sensitivity of the effective coefficients on density. Thus the instability originates from the fact that the effective diffusivity 375 and/or the effective conductivity of snow depend on ice volume fraction (or snow density) and it is triggered by transitions in snow density, i.e. layers.
We can compare the prediction of the stability analysis with the simulation of the full model for the Gaussian crust by analyzing the growth of amplitude of Fourier modes for wave number k through the Power spectral density of the perturbation.
To this end we have computed the third component of the perturbation vector u (1)

Sensitivity studies
Finally we conducted a few sensitivity studies to facilitate a discussion about the relevance of our findings in future snow 390 modeling.
To confirm that observed wave patterns in the Gaussian crust are robust against variations of absolute values in temperature gradients and initial crust density we simulated a denser Gaussian crust with different values for the temperature gradient. Figure (7)(a) shows that the non-linear PDE system essentially obeys time vs temperature gradient scaling. Virtually the same wave solution is obtained at different physical times t which scale as t ∼ (∇T ) −1 . This test also shows that for a higher crust 395 density, the density depletion on the sublimating side is more pronounced.
Second we subjected a smoothly varying density profile to a temperature gradient of 50K/m and conducted a long simulation over 170 days. The results are shown in Fig. (7)(b). This confirms that if the density profile is sufficiently smooth, even an entire season simulation shows a stable, smooth advection of the ice phase.
Finally, we mimic the situation of a snowpack over a dry soil where no vapor can enter the system from below by imposing  From a mere physical point of view, our comparison of (Calonne et al., 2014) and (Hansen and Foslien, 2015) has shown that both schemes yield similar results as long as the same parametrizations are used as closure for the PDE coefficients. The impact 410 of purely diffusive vapor transport on macroscopic density changes is rather small (Fig. (7), in agreement with (Jafari et al., 2020). Accurate estimates will certainly rely on the choice for the parametrization of the transport coefficients which naturally cause differences in the results (Fig. (1)). In homogeneous parts of the snowpack these differences are small, larger differences can be expected at layer transition regions (Fig. (2)) where phase changes sensitively react to the underlying homogenized process model equations.

415
For the present work, however, the precise numbers were of minor importance. The primary goal was an assessment of the non-linearity that is contained in published vapor homogenization schemes and their numerical requirements in view of previously reported numerical issues Jafari et al. (2020); Adams and Brown (1990); Hansen and Foslien (2015). While we have certainly challenged the numerical scheme by predominantly exploring high temperature and high density gradients, the observed time -gradient scaling (Fig. (7)) indicates that the underlying non-linear mechanisms are generally robust against 420 these details.

Oscillatory behavior in the numerical solution
We have discovered that both models (Calonne et al., 2014) and (Hansen and Foslien, 2015) exhibit oscillations in the numerical solution if they are dynamically coupled to an evolving ice phase (Eq. (7)). From our analysis of different combinations of initial and boundary conditions (Sec. 4) we conclude that two type of oscillations can occur, node-to-node oscillations (indicating 425 numerical problems) and smooth oscillations (features of the underlying equations) In view of numerical problems, we have shown that the dynamic coupling of the ice phase to heat and mass transport is equivalent to an advection of the ice phase (Sec 5). The non-linear nature of this advection causes a self-amplified increase of density changes/layer transitions, imposing special requirements on mesh resolution, time stepping and potentially shockcapturing schemes (Shu and Osher, 1988) to avoid oscillations. The fact that node-to-node oscillations occur on the "outflow" 430 boundary of the ice phase supports this. In most of our simulations we have imposed a Dirichlet boundary condition of the vapor phase which is strictly in equilibrium on the boundary. By virtue of Eq. (23) this implies that the snow density cannot change directly at the boundary, while in the interior of the domain the ice is piling up towards the boundary under the advection (Fig. 4) This is reminiscent for outflow boundaries in fluid dynamics for high local Péclet numbers (Donea and Huerta, 2003): The disagreement between the information transported through the domain from the cold boundary and the 435 information prescribed on the warm boundary is numerically resolved by steep node-to-node oscillations (Fig. 3. Similar issues can arise whenever abrupt changes in the advection velocity cannot sufficiently be resolved on a given mesh, as in Figure 5 in the middle of the domain. This interpretation of node-to-node oscillations is supported by the fact that oscillations can be completely suppressed on the lower boundary by choosing a boundary conditions ∇ρ v = ∂/∂T ρ eq v ∇T for the vapor equation, which implies a vanishing gradient for the ice profile at the boundary. Our tests indicated that node-to-node oscillations can 440 be largely reduced by choosing a high resolution in space and time. A more efficient approach might be the utilization of stabilization methods where our approximation (Eq. (23)) may serve to provide a local advection velocity of the ice phase which needs to be supplied e.g. to a SUPG stabilization scheme (Donea and Huerta, 2003).
As a nasty coincidence, after increasing mesh and time resolution to suppress numerical problems, the true solution converges to a smooth, travelling wave pattern (Fig. 5, Fig. 6). As confirmed by the theoretical analysis, these oscillations are 445 true features of the non-linear homogenization equations and triggered by gradients in the density (Sec.6.2). These patterns can only form when the ice equation is dynamically coupled in a mass conserving way to the heat and vapor equation. At the boundary, these physical oscillations are triggered as a consequence of the Dirichlet boundary condition as explained above: The competition of a fixed value of φ i directly on the boundary for z = 0 (implication of the imposed boundary conditions on the vapor phase for Eq. 7) and the increase of φ i in the interior of the domain gives rise to a transition layer at the boundary 450 with a high density gradient in the ice phase that triggers the instability according to the mechanism revealed in Sec. 6.2. This is again supported by the fact that the physical oscillations at the boundary can be suppressed by choosing the BC for the vapor equations such that the gradient in ice volume fraction vanishes. This behavior is reminiscent of the wave instability triggered by a Dirichlet boundary condition in a 1D non-linear reaction-advection-diffusion system (Vidal-Henriquez et al., 2017) which travels up-stream into the domain.

Weak-layer formation by wave instability?
We have shown that spatial heterogeneities in the density/ice volume fraction can amplify under the coupled thermodynamic description in snow. This has been pointed out before Adams and Brown (1990). We have investigated this phenomenon within the idealized scenario of a Gaussian crust where density gradients self-amplify under the non-linear advective dynamics with a subsequent instability and the emergence of wave patterns. The eigenvalues of the linearized PDE system in our linear stability 460 analysis has revealed the mechanism: A non-homogeneous density paired with the density dependence of the effective diffusion coefficients triggers the instability. Pattern formation following an instability is ubiquituous in non-linear (diffusion-reaction) PDEs (Cross and Hohenberg, 1993).
In the present case, the observed instability may have far reaching consequences which comes as a (surprising) side-product of our numerical study: As the instability causes the depletion of density on the sublimating side of a crust (Fig. 7(a)  components used here, predict a considerable reduction of the density in a sub-millimeter layer above the crust (Fig. 7(a). This mechanism is solely a consequence of the continuum description of snow. We stress, that this is different from a previously suggested reasoning (Colbeck and Jamieson, 2001) and it may occur as a superimposed effect on additional microstructural controls through near crust metamorphism (Hammonds et al., 2015). In the future, it will be therefore important to assess 470 how this instability is affected by adding mechanical settling, an evolving microstructure or enhanced mass transport from convection. First, as detailed in (Calonne et al., 2014), the homogenized equations (1), (2) are not valid for arbitrary growth rates. This is a tricky situation since the PDE system exhibits self-propelled dynamics into a state of fast growth for high density gradients, thereby leaving its own range of validity. However, since the the homogenization scheme (Calonne et al., 2014)  to sub-millimeter scales requiring ridiculously small mesh sizes to resolve them, clearly interfering with characteristic length scales of the microstructure. What is a consistent way of dealing with this situation? Resolving them numerically and averaging them out? Suppressing them numerically by artificial diffusion? Or does this behavior signal true multi-scale effects where the assumed "seperation of scales" fundamentally breaks down? Answering these questions appears to be a key demand for future work on homogenization to provide a robust recipe how to use derived equations correctly for applied modeling. Any snowpack contains variations in its properties that may be described either by continuous profiles containing large gradients (as done here) or by a discontinuous stacking of layers. We want to point out here that if a continuous description of density variations were to be replaced by discontinuous layers, the investigated wave problem won't disappear. In a hypothetical discontinuous layer setting, as commonly pursued in snowpack models, the mass continuity of the ice phase at the layer 495 interface would require to derive a dynamic equation (likely non-linear) equation for the migration of the layer interface on which the continuity of temperature and heat and mass fluxes were to be imposed. From the continuous description used here, no firm conclusion can be drawn on the behavior of the interface evolution between discontinuous, homogeneous layers. We hypothesize though, that the wave instability of the continuous PDE formulation may translate into an oscillatory instability for the position of the interface. In view of the mathematical overhead of tracking continuity conditions at moving interfaces, 500 we tend to recommend a continuous description in future snow modeling with numerical schemes that cope with arbitrary gradients in the properties.

Advantages and disadvantages of the numerical framework
We have used FEniCS for the stand-alone implementation to minimize the FEM implementation effort, yet retaining full control on the numerical solution. Overall FEniCS provides a convenient, modular setup for exchanging PDE coefficients (k eff , D eff ) 505 or boundary conditions, e.g. for more sophisticated exchange of vapor with the soil or the atmosphere. Alongside our study, we evaluated the FEniCS framework by comparing numerically different implementations. These experiences are shared for future reference.
We found that integration by parts in the weak formulation is not only necessary to apply Neumann boundary conditions, but also increases the precision, regardless of the order of interpolation polynomials. Operator splitting turned out to be of 510 limited value, the decrease in precision was not outweighed by a decrease in numerical complexity. As expected, a nondimensionalization of the equations and the corresponding re-scaling to values of order unity did not impact the solution, as long as solver convergence settings are adapted. Increasing the polynomial order of the test functions was equivalent to increasing the mesh resolution with the corresponding number of nodes. However, we experienced large errors, if the polynomial order of the variables of the coupled equations was not the same, e.g. solving φ i in first order, ρ v , T in second order caused large errors 515 on the entire domain and also violations of Dirichlet boundary conditions. Adding auxiliary algebraic variables (as done here for the source terms) can improve the solvability of the system, when convergence settings of the coulped system are adapted.
Using so called sub-functions in the formulation of the variational problem for coding modularity can introduce errors. We encountered deviations between a sub-function value and its hard-coded counterpart.
While FEniCS has provided an excellent numerical framework for the present study, a clear drawback of FEniCS would 520 however emerge if mechanical settling was to be considered, as a necessary follow-up extension. In the presence of settling, the ice phase conservation equation is an advection dominated problem which is notoriously difficult to solve on a Eulerian FE mesh without numerical smoothing. Remeshing is currently not supported in FEniCS. To this end we have explored another, fully different numerical route to enable a flexible coupling of transport and phase changes to mechanical settling which is presented in part 2 of this companion paper (Simson et al., 2021).

Conclusions
We have shown that the widely accepted form of homogenized vapor transport equations in snow predicts mathematical features (density waves) with interesting physical implications (weak layers) which constitute a considerable numerical challenge for future snowpack modeling. the dependence of the effective heat and vapor diffusion coefficients (k eff , D eff ) on snow density. Since this dependence is the most fundamental non-linearity in coupled heat and vapor transport in snow, it is unlikely that this effect will luckily disappear when considering more complex, non-linear extensions of the model. The instability is a true feature of published equations and comes at play when the ice phase is dynamically coupled to the vapor phase by phase changes in a mass conserving way.
The instability is triggered by high density gradient regions which either (pre-)exist in snowpacks in the form of layers 535 or which are generated from the self-amplification of density gradients under the coupled dynamics. This amplification is a consequence of the effective advection of the ice phase due to phase changes. This is explained within the derivation of an approximate, non-linear and non-local advection equation. Given the observed (approximate) equivalence between time and imposed temperature gradient, the system always undergoes a self-propelled evolution into its own instability. The instability might be practically irrelevant as long as smooth density profiles are considered. But the instability will certainly become 540 relevant if a snowpack model should be applicable to simulate a sublimating side of a crust as a potential origin of weak layer formation.
We have outlined open questions and limitations of the present study related to the homogenization scheme, the numerical scheme and the concept of discontinuous layers. While the present study required a stand-alone numerical implementation, it seems to be of key importance that future snow models will be flexible enough for conducting advanced numerical studies.

545
Only then re-implementations of stand-alone numerical experiments will become obsolete and the rich, non-linear behavior of snow could be predicted from a snow model alone.
Code and data availability. Will be made available upon acceptance
For the effective heat capacity we used a volume averaged formulation given in (Calonne et al., 2014) with C i , C a , ρ i given in Table 1 555 In the Calonne case we used for the effective diffusion parameters Eq. 12 from Calonne et al. (2011) and the self-consistent estimate used in Calonne et al. (2014) given by with k 0 = 0.024, k 1 = −1.23 10 −4 , k 2 = 2.5 10 −6 and D 0 given in Table 1 560 In the Hansen case we used Eq. (87) and (88) from (Hansen and Foslien, 2015) given by and k a , k i , L, D 0 given in Table 1 565

Appendix B: Analysis of residuals
Due to the lack of an analytical solution for the non-linear problem, we used the nodal residuals as an indicator for the solution error. FEniCS does not provide access to the nodal residuals to verify convergence in the NonlinearVariationalSolver class. To this end we recovered the residuals manually by decoupling the equations and iterating over individually fine-tuned solvers for each equation. The results are shown in Figure B1. While in this solution scheme the residuals are several orders of Figure B1. Results upon lowering the residuals of all equations. While some pattern is still visible, the order of magnitude is small enough to assume small errors on the solution as well.
magnitude smaller than in our regular setup, the onset of the wave instability remains the same ( Figure B2).

Appendix C: Perturbation expansion of the heat-vapor-ice system
Carrying out the perturbation expansion of the PDE system (30) around the stationary state (31) using the ansatz (37) requires to derive two auxiliary relations which are stated below. where we use the following shorthand notation for the matrix coefficients Now we can state the expansion of the PDE system (30) to first order in u as follows: 585 C (0) + (e 3 · u (1) )C (1) ∂ t u (0) + ∂ t u (1) − K (0) + (e 3 · u (1) )K (1) ∇ 2 z u (0) + ∇ 2 z u (1) − (e 3 · ∇ z u (0) ) + (e 3 · ∇ z u (1) ) K (1) + (e 3 · u (1) )K (2) ∇ z u (0) + ∇ z u (1) = Ru (0) + Ru (1) + f (C6) Equating zero and first order terms in u (0) provides an equation that is satisfied by the stationary state as it should be. By 590 collecting terms linear in u (1) the governing equation reads The advection term can be re-written in the form K (1) ∇ z u (0) ⊗ e 3 ∇ z u (1) and after multiplying the corresponding matrices we arrive at the final, linear PDE system with constant coefficients In the latter step we employed the expansions of the diffusion coefficients around the reference volume fraction φ i,0 in the form Generally, it seems feasible (though tedious) to carry out an expansion of (C8) in terms of the eigenfunctions of the differential operator satisfying the BC to incorporate the BC in the stability analysis (since u (0) already satisfies the original Dirichlet conditions, the perturbation u (1) must vanish on the boundaries). However, we limit ourselves here to the simpler case of a stability analysis in an infinite domain and take continuous Fourier transforms of (30) with respect to z. Denoting the Fourier variable by k we end up with the final result stated in Eq. (38) 610 Author contributions. The study was designed by HL and JK. Implementation, simulations, and figures were done by KS with contributions and supervision from HL and JK. Derivation of the quasi-advection and the linear stability analysis was done by HL. The manuscript was written by HL with contributions from KS and JK.
Competing interests. The authors declare that they have no conflict of interest