Multiscale modeling of heat and mass transfer in dry snow: influence of the condensation coefficient and comparison with experiments

Temperature gradient metamorphism in dry snow is driven by heat and water vapor transfer through snow, which includes conduction/diffusion processes in both air and ice phases, as well as sublimation and deposition at the ice–air interface. The latter processes are driven by the condensation coefficient α, a poorly constrained parameter in the literature. In the present paper, we use an upscaling method to derive heat and mass transfer models at the snow layer scale for values of α in the range 10⁻¹⁰ to 1. A transition α value arises, of the order of 10⁻⁴, for typical snow microstructures (characteristic length ∼ 0.5 mm), such that the vapor transport is limited by sublimation–deposition below that value and by diffusion above it. Accordingly, different macroscopic models with specific domains of validity with respect to α values are derived. A comprehensive evaluation of the models is presented by comparison with three experimental datasets, as well as with pore-scale simulations using a simplified microstructure. The models reproduce the two main features of the experiments: the non-linear temperature profiles, with enhanced values in the center of the snow layer, and the mass transfer, with an abrupt basal mass loss. However, both features are underestimated overall by the models when compared to the experimental data. We investigate possible causes of these discrepancies and suggest potential improvements for the modeling of heat and mass transport in dry snow.

1 Introduction

Natural snowpacks are frequently subjected to temperature gradients induced by meteorological conditions. In the case of temperature gradient in dry snow, heat and water vapor are transported through the snowpack by heat conduction through ice and air and by vapor diffusion in air. These phenomena are coupled by the sublimation–deposition processes at the ice–air interfaces. In practice, such transfer processes can be enhanced by natural air convection induced by the temperature gradient (e.g., Jafari et al., 2022) or by forced convection generated by the wind at the snowpack surface (e.g., Albert and McGilvary, 1992; Calonne et al., 2015). For the sake of simplicity, both types of convection are neglected in the following. All those processes lead to changes in the snow microstructure called the temperature gradient metamorphism (TGM), which transforms snow into faceted crystals (FC) in the case of moderate gradients and into depth hoar (DH) for stronger gradients (see Fierz et al., 2009). Those transformations of the microstructure can sometimes come along with a redistribution of mass in the snow layer, a density drop, or even the formation of an air gap at the base of the snowpack, as observed in the Arctic (e.g., Domine et al., 2019) or in some cold room experiments (e.g., Kamata and Sato, 2007; Wiese, 2017; Bouvet et al., 2023). As a result of changes in microstructure and density, the TGM also induces significant changes in the snow physical and mechanical properties, such as thermal conductivity, vapor diffusivity, or elastic properties (e.g., Srivastava et al., 2010; Calonne et al., 2014a; Wautier et al., 2015), affecting the snowpack behavior at larger scale. Hence, an accurate representation of the heat and mass transport processes during the TGM is key to accurately model the snow cover, as required for many applications such as avalanche forecasting or climate studies (Jordan, 1991; Lehning et al., 2002; Vionnet et al., 2012).

Models to describe the heat and mass transfer at the pore scale, referred to as the micro-scale, have been proposed (e.g., Flin et al., 2003; Flin and Brzoska, 2008; Kaempfer and Plapp, 2009; Vetter et al., 2010; Bouvet et al., 2022). Based on explicit representations of the 3D snow microstructure, often from X-ray tomography images, simulations at that scale are usually performed on small snow volumes due to numerical cost limitations. In micro-scale modeling, heat and mass transfer processes are coupled through interface conditions that account for the sublimation and deposition processes at the air–ice interface and involve an interface growth velocity. In snow physics, the interface growth velocity is classically given by the Hertz–Knudsen equation and strongly depends on the condensation coefficient α, also called attachment, sticking, deposition, or kinetic coefficient (e.g., Flin et al., 2003; Libbrecht, 2005; Brzoska et al., 2008; Kaempfer and Plapp, 2009; Furukawa, 2015; Krol and Löwe, 2016; Fourteau et al., 2021a; Granger et al., 2021). This coefficient describes the probability that a water vapor molecule striking the ice surface will be incorporated into it and that it theoretically ranges from 0 to 1 for an infinite flat surface (see, e.g., Libbrecht, 2005; Furukawa, 2015). The analogous coefficient for sublimation can also be defined but is often assumed equal to the condensation coefficient. At present, the condensation coefficient is still poorly understood and quantified, notably because of its complex dependencies on temperature, supersaturation (or temperature gradient), and ice crystalline orientation (see, e.g., Libbrecht, 2021). Estimates of the condensation coefficient can be found in the literature. Typical values obtained from single-ice-crystal growth experiments range from 10⁻⁴ to 10⁻¹ (see, e.g., Libbrecht and Rickerby, 2013). Indirect estimates from snow modeling at the pore scale range from 10⁻⁴ to 10⁻³ (see, e.g., Flin, 2004; Bouvet et al., 2022).

In the last few decades, several models have been presented to describe the heat and mass transfer at the scale of a snow layer, referred to as the macro-scale. At that scale, the snow microstructure is not explicitly represented, and simulations can be carried out on entire snowpacks. The first models assumed saturated vapor conditions in the snow (e.g., de Quervain, 1963; Anderson, 1976; Powers et al., 1985). Later, using a phenomenological approach, Albert and McGilvary (1992) proposed describing the heat and water vapor transfer through a snowpack subjected to an airflow but without restricting the water vapor to its saturation value. The model uses two coupled advection–diffusion equations, including a source term arising from phase change at the pore scale. A similar heat and mass transfer model was analytically obtained by Calonne et al. (2014b, 2015), using an upscaling method. In that case, the macroscopic equivalent modeling was derived from its description at the pore scale using the homogenization of multiple-scale expansions. This theoretical method also provides the exact expression of the effective parameters arising at the macro-scale and the domains of validity of the macroscopic modeling. Two main effective parameters emerge from the model: (i) the effective thermal conductivity k^eff, which depends on the ice and air conductivity and on the snow microstructure, and (ii) the effective diffusion D^eff, which depends on the vapor molecular diffusion coefficient and the snow microstructure. The source term is related to the Hertz–Knudsen equation, which involves the condensation coefficient α. Calonne et al. (2014b) have shown that this model is valid for interface growth velocities below 3 × 10⁻¹¹ m s⁻¹, which typically corresponds to slow kinetics.

Other approaches largely rely on the assumption of saturated vapor conditions, which seems valid for faster kinetics and rather high values of α (e.g., Sturm and Benson, 1997; Kamata and Sato, 2007; Hansen and Foslien, 2015). Hansen and Foslien (2015) developed a heat and mass transfer model using a mixture theory. Assuming that the water vapor is saturated (based on the value of α=0.0144 from Delaney et al. (1964)), the authors derived a unique thermal equation which yields an apparent thermal conductivity that depends on the air and ice conductivities, the water vapor diffusivity, the latent heat of sublimation–deposition, and the temperature derivative of the Clapeyron equation. A similar formulation of the apparent thermal conductivity was also proposed by Yosida et al. (1955). Recently, Fourteau et al. (2021b) investigated the influence of α on the apparent diffusion coefficient in snow. By performing numerical simulations on 3D images, they showed that this apparent diffusion coefficient is equal to D^eff for α values smaller than $\approx {\mathrm{10}}^{-\mathrm{4}}$ and then increases with increasing α until it reaches a plateau for α values larger than 10⁻², with the value at the plateau smaller than the molecular diffusion of water vapor in the air. In a companion paper, Fourteau et al. (2021a) computed from 3D images the apparent thermal conductivity of snow, assuming that the water vapor at the ice–air interface is equal to the water vapor at the saturation point given by the Clapeyron equation. In this case, they showed that the apparent thermal conductivity is enhanced by the sublimation–deposition process arising at the pore scale. Their results are consistent with the model of Moyne et al. (1988) for wet porous media, based on the same hypothesis at the micro-scale and derived using the volume averaging method.

Further uses of the abovementioned models, as well as their implementation in full snow cover models, are limited by some challenges. One is the difficulty in choosing between models, as they differ in many ways. They were derived using different methods; involve different balance equations and effective parameters; and are valid for different, often unclear, domains of validity in terms of α values. This should be clarified, especially by estimating the α values for which the assumption of saturated water vapor is theoretically valid. A second challenge is that none of these models was thoroughly evaluated to assess their performances. This might be partly due to the limited number of suited datasets to compare with. The datasets from the cold-laboratory experiments of Kamata and Sato (2007) and, recently, of Bouvet et al. (2023), however, seem appropriate for such comparisons, as they provide the time series of the vertical profiles of snow density and temperature, as well as the forcing conditions to be reproduced in the simulations.

This paper aims (i) to define the heat and mass transport modeling in dry snow for the full range of α values and (ii) to evaluate the model's ability to reproduce natural snow evolution during the TGM. To this end, first, the homogenization of multiple-scale expansions is applied to derive the macroscopic equivalent modeling of heat and vapor transfer for α values ranging from 10⁻¹⁰ to 1, following Calonne et al. (2014b). The physics considered at the pore scale includes heat conduction, vapor diffusion, and phase change and neglects any transport linked to the curvature effect and convection. The macroscopic models and the involved macroscopic properties are compared to the ones from the literature and are illustrated for simplified snow microstructures. Second, the derived macroscopic models are evaluated using three cold-laboratory experiments of TGM from Kamata and Sato (2007) and Bouvet et al. (2023). The experiments are reproduced with the macroscopic models, and results between observations and simulations are analyzed.

2 Derivation of the macroscopic modeling

2.1 Upscaling method

We apply the homogenization technique of multiple-scale expansions (Bensoussan et al., 1978; Sanchez-Palencia, 1980) to the physics of heat and vapor transport in dry snow. The homogenization method allows us to model the local physical processes in heterogeneous media by an equivalent continuous macroscopic description if the condition of separation of scales is satisfied (Bensoussan et al., 1978; Sanchez-Palencia, 1980; Auriault, 1991; Auriault et al., 2009). This separation of scales can be expressed as $\mathit{\epsilon }=l/L\ll \mathrm{1}$, where l and L are the characteristic lengths of the heterogeneities at the pore scale and of the macroscopic sample or excitation, respectively. This condition implies the existence of a representative elementary volume (REV) of size l for both the material and the excitation. Following the methodology presented by Auriault (1991), the macroscopic equivalent model is obtained from the description of the physics at the pore scale by (i) assuming the medium to be periodic, without the loss of generality, as the condition $\mathit{\epsilon }=l/L\ll \mathrm{1}$ is fulfilled; (ii) writing the description of the physics at the pore scale in a dimensionless form; (iii) evaluating the obtained dimensionless numbers with respect to the coefficient of separation of scale ε; (iv) looking for the unknown fields in the form of asymptotic expansions in powers of ε; and (v) solving the successive boundary value problems that are obtained after introducing these expansions in the pore-scale dimensionless description. The macroscopic equivalent model is obtained from compatibility conditions that are the necessary conditions for the existence of solutions to the boundary value problems.

2.2 Physical processes at the pore scale

As in Calonne et al. (2014b), we assume that a snow layer of the characteristic length L can be represented by a collection of spatially periodic REVs of the characteristic length l, such that the coefficient of separation of scale $\mathit{\epsilon }=l/L\ll \mathrm{1}$.

Figure 1 Physical phenomena under consideration at the representative elementary volume (REV) scale.

In what follows, Ω is the REV domain, Ω_i is the ice domain, and Ω_a is the air domain (Fig. 1). The ice grain interface is denoted Γ, and n_i is the unit's outward vector of Ω_i. The subscripts i or a are related to quantities defined in Ω_i and Ω_a, respectively. As illustrated in Fig. 1, the processes of heat and mass transport in dry snow considered are (i) the heat conduction through ice and air; (ii) the water vapor diffusion in air; and (iii) the sublimation of ice and the deposition of vapor at the ice grain interface, characterized by an interface growth velocity (Libbrecht, 2005; Kaempfer and Plapp, 2009; Barrett et al., 2012), following the Hertz–Knudsen equation. This latter equation, initially derived to describe the condensation–evaporation processes at a liquid–gas interface, is widely used in snow physics and is supported by a great deal of experimental evidence (e.g., Libbrecht, 2005; Kaempfer and Plapp, 2009; Furukawa, 2015; Libbrecht and Rickerby, 2013; Krol and Löwe, 2016. Air convection and snow densification are not taken into account here. Assuming that the properties of air and ice are isotropic, these physical processes at the pore scale are described by the following set of equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{i}}{C}_{\mathrm{i}}{\displaystyle \frac{\partial T_{\mathrm{i}}}{\partial t}}-\mathrm{div}\left(k_{\mathrm{i}}\mathbf{grad}{T}_{\mathrm{i}}\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{i}},\text{(2)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{a}}{\displaystyle \frac{\partial T_{\mathrm{a}}}{\partial t}}-\mathrm{div}\left(k_{\mathrm{a}}\mathbf{grad}{T}_{\mathrm{a}}\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\text{(3)}& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}}{\partial t}}-\mathrm{div}\left(D_{\mathrm{v}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{v}}\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\text{(4)}& {\displaystyle }{T}_{\mathrm{i}}=T_{\mathrm{a}}\mathrm{on}\mathrm{\Gamma },\text{(5)}& {\displaystyle }{k}_{\mathrm{i}}\mathbf{grad}{T}_{\mathrm{i}}\cdot {\mathit{n}}_{\mathrm{i}}-k_{\mathrm{a}}\mathbf{grad}{T}_{\mathrm{a}}\cdot {\mathit{n}}_{\mathrm{i}}=L_{\mathrm{sg}}\mathit{w}\cdot {\mathit{n}}_{\mathrm{i}}\mathrm{on}\mathrm{\Gamma },\text{(6)}& {\displaystyle }{D}_{\mathrm{v}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{v}}\cdot {\mathit{n}}_{\mathrm{i}}=({\mathit{\rho }}_{\mathrm{i}}-{\mathit{\rho }}_{\mathrm{v}})\mathit{w}\cdot {\mathit{n}}_{\mathrm{i}}\simeq {\mathit{\rho }}_{\mathrm{i}}\mathit{w}\cdot {\mathit{n}}_{\mathrm{i}}\mathrm{on}\mathrm{\Gamma },\end{array}$$

where t is the time (s), T is the temperature (K), k is the thermal conductivity (W m⁻¹ K⁻¹), ρ is the density (kg m⁻³), C is the specific heat capacity (J kg⁻¹ K⁻¹), L_sg is the latent heat of sublimation–deposition (J m⁻³), w is the interface growth velocity (m s⁻¹), ρ_v is the partial density of water vapor in air (kg m⁻³), D_v is the water vapor diffusion coefficient in air (m² s⁻¹), and div and grad are the divergence and gradient operators with respect to the physical space variable X, respectively. At the interface, the heat and mass transfer are coupled through the normal interface growth velocity, $w_{\mathrm{n}}=\mathit{w}\cdot {\mathit{n}}_{\mathrm{i}}$, which is given by the Hertz–Knudsen equation,

$$\begin{array}{}\text{(7)}& w_{\mathrm{n}}=\mathit{w}\cdot {\mathit{n}}_{\mathrm{i}}={\displaystyle \frac{\mathrm{1}}{\mathit{\beta }}}\left[{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{v}}-{\mathit{\rho }}_{\mathrm{vs}}\left(T_{\mathrm{a}}\right)}{{\mathit{\rho }}_{\mathrm{vs}}\left(T_{\mathrm{a}}\right)}}-d_{\mathrm{0}}K\right]\mathrm{on}\mathrm{\Gamma },\end{array}$$

such that w_n is positive when the ice grain grows and negative when it sublimates. β is the interface of the kinetic coefficient (s m⁻¹), ρ_vs is the saturation of the water vapor density in air (kg m⁻³), d₀ is the capillary length (m), and K is the interface of the mean curvature (m⁻¹). The interface kinetic coefficient β is linked to the condensation coefficient α by

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\beta }}}=\mathit{\alpha }{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{vs}}\left(T_{\mathrm{a}}\right)}{{\mathit{\rho }}_{\mathrm{i}}}}\sqrt{{\displaystyle \frac{k_{\mathrm{B}}{T}_{\mathrm{a}}}{\mathrm{2}\mathit{\pi }m}}},\end{array}$$

where m is the mass of a water molecule (kg), and k_B is the Boltzmann constant equal to 1.38 × 10⁻²³ J K⁻¹. As already mentioned, the condensation coefficient α characterizes the probability of a water molecule hitting the surface of the solid to be incorporated to the crystal, or vice versa, and ranges from 0 to 1. Although this coefficient depends on several parameters as temperature, supersaturation, and crystalline orientation, we assume that this parameter is constant over the REV at the first order. The saturation vapor density ρ_vs at a given air temperature T_a is given by the Clausius–Clapeyron law as follows:

$$\begin{array}{}\text{(9)}& {\mathit{\rho }}_{\mathrm{vs}}\left(T_{\mathrm{a}}\right)={\mathit{\rho }}_{\mathrm{vs}}^{\mathrm{ref}}\left(T^{\mathrm{ref}}\right)\exp \left[{\displaystyle \frac{L_{\mathrm{sg}}m}{{\mathit{\rho }}_{\mathrm{i}}{k}_{\mathrm{B}}}}\left({\displaystyle \frac{\mathrm{1}}{T^{\mathrm{ref}}}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{a}}}}\right)\right].\end{array}$$

In the current work, we chose the reference values T^ref=263 K, leading to a ${\mathit{\rho }}_{\mathrm{vs}}^{\mathrm{ref}}\left(T^{\mathrm{ref}}\right)$ value of 2.173 × 10⁻³ kg m⁻³. For simplicity, we assume that none of the material properties (ρ, C, k_B, D_v, β, and m) depend on the temperature. Also, the effect of curvature on the ice interface growth is considered insignificant compared to the effect of temperature and is thus neglected. Consequently, using Eq. (8), the Hertz–Knudsen equation can be rewritten as

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{w}_{\mathrm{n}}& =\mathit{w}\cdot {\mathit{n}}_{\mathrm{i}}={\displaystyle \frac{\mathrm{1}}{\mathit{\beta }{\mathit{\rho }}_{\mathrm{vs}}\left(T_{\mathrm{a}}\right)}}\left[{\mathit{\rho }}_{\mathrm{v}}-{\mathit{\rho }}_{\mathrm{vs}}\left(T_{\mathrm{a}}\right)\right]\\ & ={\displaystyle \frac{\mathit{\alpha }}{{\mathit{\rho }}_{\mathrm{i}}}}{w}_{\mathrm{k}}\left(T_{\mathrm{a}}\right)\left[{\mathit{\rho }}_{\mathrm{v}}-{\mathit{\rho }}_{\mathrm{vs}}\left(T_{\mathrm{a}}\right)\right]\mathrm{on}\mathrm{\Gamma },\end{array}\end{array}$$

where $w_{\mathrm{k}}=\sqrt{k_{\mathrm{B}}{T}_{\mathrm{a}}/\mathrm{2}\mathit{\pi }m}$ is defined as a kinetic velocity which depends on the temperature at the ice–air interface. Taking this result into account, Eqs. (5) and (6) can be rewritten as

$$\begin{array}{}\text{(11)}& {\displaystyle }\begin{array}{rl}{k}_{\mathrm{i}}\mathbf{grad}{T}_{\mathrm{i}}& \cdot {\mathit{n}}_{\mathrm{i}}-k_{\mathrm{a}}\mathbf{grad}{T}_{\mathrm{a}}\cdot {\mathit{n}}_{\mathrm{i}}\\ & =L_{\mathrm{sg}}{\displaystyle \frac{\mathit{\alpha }}{{\mathit{\rho }}_{\mathrm{i}}}}{w}_{\mathrm{k}}\left(T_{\mathrm{a}}\right)\left[{\mathit{\rho }}_{\mathrm{v}}-{\mathit{\rho }}_{\mathrm{vs}}\left(T_{\mathrm{a}}\right)\right]\\ & =L_{\mathrm{sg}}{\displaystyle \frac{D_{\mathrm{v}}}{{\mathit{\rho }}_{\mathrm{i}}}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{v}}\cdot {\mathit{n}}_{\mathrm{i}}\mathrm{on}\mathrm{\Gamma },\end{array}\text{(12)}& {\displaystyle }{D}_{\mathrm{v}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{v}}\cdot {\mathit{n}}_{\mathrm{i}}=\mathit{\alpha }{w}_{\mathrm{k}}\left(T_{\mathrm{a}}\right)\left[{\mathit{\rho }}_{\mathrm{v}}-{\mathit{\rho }}_{\mathrm{vs}}\left(T_{\mathrm{a}}\right)\right]\mathrm{on}\mathrm{\Gamma }.\end{array}$$

2.3 Dimensionless pore-scale description

The next step is the normalization of the above pore-scale description in Eqs. (1)–(4) and (11)–(12). For that reason, all the dimensional variables in this description are written such that each variable φ reads $\mathit{\phi }={\mathit{\phi }}_{\mathrm{c}}{\mathit{\phi }}^{*}$, where the subscript c denotes a characteristic quantity (constant), and the asterisk superscript denotes a dimensionless variable. Note that the microscopic length l is chosen as characteristic length, such that l_c=l; i.e., the so-called microscopic point of view is adopted (Auriault, 1991). The formal dimensionless set of equations that describes the physics at the pore scale can thus be written as

$$\begin{array}{}\text{(13)}& {\displaystyle }\left[F_{\mathrm{i}}^T\right]{\mathit{\rho }}_{\mathrm{i}}^{*}{C}_{\mathrm{i}}^{*}{\displaystyle \frac{\partial T_{\mathrm{i}}^{*}}{\partial t^{*}}}-{\mathrm{div}}^{*}\left(k_{\mathrm{i}}^{*}{\mathbf{grad}}^{*}{T}_{\mathrm{i}}^{*}\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{i}},\text{(14)}& {\displaystyle }\left[F_{\mathrm{a}}^T\right]{\mathit{\rho }}_{\mathrm{a}}^{*}{C}_{\mathrm{a}}^{*}{\displaystyle \frac{\partial T_{\mathrm{a}}^{*}}{\partial t^{*}}}-{\mathrm{div}}^{*}\left(k_{\mathrm{a}}^{*}\mathbf{grad}{T}_{\mathrm{a}}^{*}\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\text{(15)}& {\displaystyle }\left[F_{\mathrm{a}}^{\mathit{\rho }}\right]{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}^{*}}{\partial t^{*}}}-{\mathrm{div}}^{*}\left(D_{\mathrm{v}}^{*}{\mathbf{grad}}^{*}{\mathit{\rho }}_{\mathrm{v}}^{*}\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\text{(16)}& {\displaystyle }{T}_{\mathrm{i}}^{*}=T_{\mathrm{a}}^{*}\mathrm{on}\mathrm{\Gamma },\text{(17)}& {\displaystyle }\begin{array}{rl}\left[K\right]& k_{\mathrm{i}}^{*}{\mathbf{grad}}^{*}{T}_{\mathrm{i}}^{*}\cdot {\mathit{n}}_{\mathrm{i}}-k_{\mathrm{a}}^{*}{\mathbf{grad}}^{*}{T}_{\mathrm{a}}^{*}\cdot {\mathit{n}}_{\mathrm{i}}\\ & =\left[H\right]L_{\mathrm{sg}}^{*}{\displaystyle \frac{D_{\mathrm{v}}^{*}}{{\mathit{\rho }}_{\mathrm{i}}^{*}}}{\mathbf{grad}}^{*}{\mathit{\rho }}_{\mathrm{v}}^{*}\cdot {\mathit{n}}_{\mathrm{i}}\mathrm{on}\mathrm{\Gamma },\end{array}\text{(18)}& {\displaystyle }{D}_{\mathrm{v}}^{*}{\mathbf{grad}}^{*}{\mathit{\rho }}_{\mathrm{v}}^{*}\cdot {\mathit{n}}_{\mathrm{i}}=\left[W_{\mathrm{R}}\right]{\mathit{\alpha }}^{*}{w}_{\mathrm{k}}^{*}\left[{\mathit{\rho }}_{\mathrm{v}}^{*}-\left[R\right]{\mathit{\rho }}_{\mathrm{vs}}^{*}\left(T_{\mathrm{a}}^{*}\right)\right]\mathrm{on}\mathrm{\Gamma }.\end{array}$$

This dimensionless description introduces seven dimensionless numbers that characterize the relative intensity of the physical processes at the pore scale. These dimensionless numbers are defined as

$$\begin{array}{rl}\left[F_{\mathrm{i}}^T\right]& ={\displaystyle \frac{l^{\mathrm{2}}{\mathit{\rho }}_{i_{\mathrm{c}}}{C}_{i_{\mathrm{c}}}}{t_{\mathrm{c}}{k}_{i_{\mathrm{c}}}}},\left[F_{\mathrm{a}}^T\right]={\displaystyle \frac{l^{\mathrm{2}}{\mathit{\rho }}_{a_{\mathrm{c}}}{C}_{a_{\mathrm{c}}}}{t_{\mathrm{c}}{k}_{a_{\mathrm{c}}}}},\\ & \left[F_{\mathrm{a}}^{\mathit{\rho }}\right]={\displaystyle \frac{l^{\mathrm{2}}}{D_{v_{\mathrm{c}}}{t}_{\mathrm{c}}}},\left[K\right]={\displaystyle \frac{k_{i_{\mathrm{c}}}}{k_{a_{\mathrm{c}}}}},\\ & \left[W_{\mathrm{R}}\right]={\displaystyle \frac{l{\mathit{\alpha }}_{\mathrm{c}}{w}_{{\mathrm{k}}_{\mathrm{c}}}}{D_{v_{\mathrm{c}}}}},\left[R\right]={\displaystyle \frac{{\mathit{\rho }}_{v{s}_{\mathrm{c}}}\left(T_{a_{\mathrm{c}}}\right)}{{\mathit{\rho }}_{v_{\mathrm{c}}}}},\end{array}$$

$$\begin{array}{}\text{(19)}& \begin{array}{rl}\left[H\right]& ={\displaystyle \frac{l{L}_{{\mathrm{sg}}_{\mathrm{c}}}{w}_{n_{\mathrm{c}}}}{k_{a_{\mathrm{c}}}{T}_{a_{\mathrm{c}}}}},\mathrm{with}{w}_{n_{\mathrm{c}}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{c}}{w}_{{\mathrm{k}}_{\mathrm{c}}}}{{\mathit{\rho }}_{i_{\mathrm{c}}}}}({\mathit{\rho }}_{v_{\mathrm{c}}}-{\mathit{\rho }}_{v{s}_{\mathrm{c}}}(T_{a_{\mathrm{c}}}\left)\right)\\ & ={\displaystyle \frac{D_{v_{\mathrm{c}}}{\mathit{\rho }}_{v_{\mathrm{c}}}}{l{\mathit{\rho }}_{i_{\mathrm{c}}}}}.\end{array}\end{array}$$

Dimensionless numbers $\left[F_{\mathrm{i}}^T\right]$ and $\left[F_{\mathrm{a}}^T\right]$ correspond to the inverse of the Fourier number in Ω_i and Ω_a, respectively. They characterize the ratio between the thermal energy storage rate and the heat conduction rate. $\left[F_{\mathrm{a}}^{\mathit{\rho }}\right]$ is an analogous inverse Fourier number for the transient water vapor transfer by diffusion in Ω_a. Dimensionless numbers [K], [R], [H], and [W_R] are defined at the ice–air interface. [H] characterizes the ratio between the heat flux induced by deposition and sublimation and the heat flux induced by conduction in the air phase. The above analysis differs slightly from the one presented in Calonne et al. (2014b). Indeed, two new dimensionless parameters are introduced, [W_R] and [R], to better capture the effect of α on the macroscopic models. Finally, let us remark that Eq. (18) defined at the ice–air interface corresponds to a Robin boundary condition, i.e., a weighted combination of a Dirichlet boundary condition and a Neumann boundary condition. Hence, when [W_R] tends towards zero, Eq. (18) is equivalent to a Neumann boundary condition ($D_{\mathrm{v}}^{*}{\mathbf{grad}}^{*}{\mathit{\rho }}_{\mathrm{v}}^{*}\cdot {\mathit{n}}_{\mathrm{i}}=\mathrm{0}$), whereas when [W_R] tends towards infinite (or is very large), Eq. (18) is equivalent to a Dirichlet boundary condition (${\mathit{\rho }}_{\mathrm{v}}^{*}={\mathit{\rho }}_{\mathrm{vs}}^{*}\left(T_{\mathrm{a}}^{*}\right)$).

2.4 Estimation of the dimensionless numbers

The next key step is to estimate the above six dimensionless numbers with respect to the separation of scale parameter $\mathit{\epsilon }=l/L$ in order to weigh the relative importance of the physical phenomena coming from the pore-scale description. In practice, l and L correspond to the order of magnitude of the typical snow grain size and the thickness of a snow layer, respectively. In what follows, we assumed that $l\approx \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ m and L≈0.1 m, leading to $\mathit{\epsilon }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$. The characteristic value of each variable in the dimensionless numbers is summarized in Table 1. These values were evaluated for a temperature of −10 °C and come from the literature (Massman, 1998; Kaempfer and Plapp, 2009). According to these characteristic values, it can be first shown (Calonne et al., 2014b) that the thermal diffusivity in the ice phase $D_{i_{\mathrm{c}}}=k_{i_{\mathrm{c}}}/\left(C_{i_{\mathrm{c}}}{\mathit{\rho }}_{i_{\mathrm{c}}}\right)$ and in the air phase $D_{a_{\mathrm{c}}}=k_{a_{\mathrm{c}}}/(C_{a_{\mathrm{c}}}{\mathit{\rho }}_{a_{\mathrm{c}}}$ ) are of the same order of magnitude as the vapor diffusion coefficient $D_{v_{\mathrm{c}}}$. Thus, the characteristic time t_c associated with these transfers through the snowpack are of the same order of magnitude, where $t_{\mathrm{c}}=\mathcal{O}(L^{\mathrm{2}}/D_{i_{\mathrm{c}}})=\mathcal{O}(L^{\mathrm{2}}/D_{a_{\mathrm{c}}})=\mathcal{O}(L^{\mathrm{2}}/D_{v_{\mathrm{c}}})$. Hence, from Eq. (19), we get $\left[F_{\mathrm{i}}^T\right]=\mathcal{O}\left(\left[F_{\mathrm{a}}^T\right]\right)=\mathcal{O}\left(\left[F_{\mathrm{a}}^{\mathit{\rho }}\right]\right)=\mathcal{O}\left({\mathit{\epsilon }}^{\mathrm{2}}\right)$. At the ice pore interface, from Eq. (19), we have [K]=𝒪(1) and [R]=𝒪(1). The latter estimation implies that the supersaturation $\mathit{\sigma }=({\mathit{\rho }}_{\mathrm{v}}-{\mathit{\rho }}_{\mathrm{vs}})/{\mathit{\rho }}_{\mathrm{vs}}$ varies between −1 and 13, which is consistent with the range of values classically considered (Libbrecht and Rickerby, 2013). The dimensionless number [W_R] can be written as

$$\left[W_{\mathrm{R}}\right]={\displaystyle \frac{l{\mathit{\alpha }}_{\mathrm{c}}{w}_{{\mathrm{k}}_{\mathrm{c}}}}{D_{v_{\mathrm{c}}}}}={\displaystyle \frac{l^{\mathrm{2}}}{D_{v_{\mathrm{c}}}}}{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{c}}{w}_{{\mathrm{k}}_{\mathrm{c}}}}{l}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}}{{\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}}},$$

where ${\mathit{\tau }}_{\mathrm{d}}=l^{\mathrm{2}}/D_{v_{\mathrm{c}}}$ is the characteristic time associated with water vapor diffusion at the pore scale, and ${\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}=l/\left({\mathit{\alpha }}_{\mathrm{c}}{w}_{{\mathrm{k}}_{\mathrm{c}}}\right)$ is the characteristic time associated with the sublimation–deposition process. This result shows that this ratio can take different orders of magnitude, depending on the value of α_c. Using the characteristic values given in Table 1, this ratio is equal to 1 for a particular value of α_c, which is denoted ${\mathit{\alpha }}_{\mathrm{T}}=D_{v_{\mathrm{c}}}/\left(l{w}_{{\mathrm{k}}_{\mathrm{c}}}\right)\approx \mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$. This value decreases when the characteristic length l increases, such that the values range between 10⁻³ for small grains (∼ 0.1 mm) and 10⁻⁵ for very large grains (∼ 5 mm). It also depends on the temperature, but the influence is negligible in the −30 to 0 °C range. The α_T value characterizes the transition between two mechanisms which drive the water vapor transfer at the pore scale. When ${\mathit{\tau }}_{\mathrm{d}}\ll {\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}$, i.e., for α_c≪α_T, the water vapor flux is limited by sublimation–deposition processes. This case is also called the “slow kinetics case” in Fourteau et al. (2021a). When ${\mathit{\tau }}_{\mathrm{d}}\gg {\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}$, i.e., for α_c≫α_T, the water vapor transfer is mainly limited by diffusion, which is called the “fast kinetics case” in Fourteau et al. (2021a). For intermediate cases, both mechanisms may be in competition.

Table 1 Characteristic values of the properties evaluated at −10 °C from the literature (Massman, 1998; Kaempfer and Plapp, 2009).

Figure 2 Estimation of the dimensionless number [W_R] with respect to α, which leads to several cases of macroscopic modeling to be considered (cases A to D2). The α_T value characterizes the transition between two cases presenting different limiting processes for the water vapor transfer at the pore scale, so that ${\mathit{\tau }}_{\mathrm{d}}<{\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}$ or ${\mathit{\tau }}_{\mathrm{d}}>{\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}$, with τ_d as the characteristic time associated with water vapor diffusion and ${\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}$ as the characteristic time associated with the sublimation–deposition process. α_T was estimated based on the characteristic values given in Table 1.

Estimations of the dimensionless number [H] are not as straightforward, as they depend on the intensity of the interface normal growth velocity $w_{n_{\mathrm{c}}}$. When α_c is small (typically smaller than α_T), [W_R] is also small, and Eq. (18) implies that $\mathrm{\Delta }{\mathit{\rho }}_{v_{\mathrm{c}}}$ has a finite value (𝒪(ρ_vc)). Thus, this dimensionless number [H] can be also written as

$$\left[H\right]={\displaystyle \frac{l{L}_{{\mathrm{sg}}_{\mathrm{c}}}{\mathit{\alpha }}_{\mathrm{c}}{w}_{{\mathrm{k}}_{\mathrm{c}}}{\mathit{\rho }}_{v_{\mathrm{c}}}}{k_{a_{\mathrm{c}}}{T}_{a_{\mathrm{c}}}{\mathit{\rho }}_{i_{\mathrm{c}}}}}.$$

In that case, it increases when α_c increases, and according to the characteristic values given in Table 1, it is of the same order as [W_R]. For large values of α_c (typically larger than α_T), Eq. (18) implies that ${\mathit{\rho }}_{\mathrm{v}}^{*}\approx {\mathit{\rho }}_{\mathrm{vs}}^{*}\left(T_{\mathrm{a}}^{*}\right)$. As a consequence, from Eqs. (17) and (19), [H] can be rewritten as

$$\left[H\right]={\displaystyle \frac{l{L}_{{\mathrm{sg}}_{\mathrm{c}}}{D}_{\mathrm{vc}}\mathit{\gamma }\left(T_{a_{\mathrm{c}}}\right)T_{a_{\mathrm{c}}}}{l{\mathit{\rho }}_{i_{\mathrm{c}}}{k}_{a_{\mathrm{c}}}{T}_{a_{\mathrm{c}}}}}={\displaystyle \frac{L_{{\mathrm{sg}}_{\mathrm{c}}}{D}_{v_{\mathrm{c}}}\mathit{\gamma }\left(T_{a_{\mathrm{c}}}\right)}{{\mathit{\rho }}_{i_{\mathrm{c}}}{k}_{a_{\mathrm{c}}}}}={\displaystyle \frac{k_{{\mathrm{dif}}_{\mathrm{c}}}}{k_{a_{\mathrm{c}}}}},$$

where $\mathit{\gamma }\left(T_{a_{\mathrm{c}}}\right)=\mathrm{d}{\mathit{\rho }}_{\mathrm{vs}}\left(T_{a_{\mathrm{c}}}\right)/\mathrm{d}{T}_{a_{\mathrm{c}}}$ is the derivative of the Clausius–Clapeyron law, and $k_{{\mathrm{dif}}_{\mathrm{c}}}=L_{{\mathrm{sg}}_{\mathrm{c}}}{D}_{\mathrm{vc}}\mathit{\gamma }\left(T_{a_{\mathrm{c}}}\right)/{\mathit{\rho }}_{i_{\mathrm{c}}}$ can be seen as an enhancement of the air thermal conductivity. Using the characteristic values given in Table 1 and the Clausius–Clapeyron equation, Eq. (9), for large values of α_c, [H]=𝒪(1). According to the above analysis, several cases must be considered, depending on the value of the condensation coefficient α_c (Fig. 2).

• Case A is ${\mathit{\tau }}_{\mathrm{d}}=\mathcal{O}\left({\mathit{\epsilon }}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}\right)$, i.e., [W_R]=𝒪(ε²) and [H]=𝒪(ε²).

• Case B is ${\mathit{\tau }}_{\mathrm{d}}=\mathcal{O}\left(\mathit{\epsilon }{\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}\right)$, i.e., [W_R]=𝒪(ε) and [H]=𝒪(ε).

• Case C is ${\mathit{\tau }}_{\mathrm{d}}=\mathcal{O}\left({\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}\right)$, i.e., [W_R]=𝒪(1) and [H]=𝒪(1).

• Case D1 is ${\mathit{\tau }}_{\mathrm{d}}=\mathcal{O}\left({\mathit{\epsilon }}^{-\mathrm{1}}{\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}\right)$, i.e., $\left[W_{\mathrm{R}}\right]=\mathcal{O}\left({\mathit{\epsilon }}^{-\mathrm{1}}\right)$ and [H]=𝒪(1).

• Case D2 is ${\mathit{\tau }}_{\mathrm{d}}=\mathcal{O}\left({\mathit{\epsilon }}^{-\mathrm{2}}{\mathit{\tau }}_{\mathrm{sub}/\mathrm{dep}}\right)$, i.e., $\left[W_{\mathrm{R}}\right]=\mathcal{O}\left({\mathit{\epsilon }}^{-\mathrm{2}}\right)$ and [H]=𝒪(1).

Cases A and B correspond to $\mathrm{0}⩽\mathit{\alpha }\ll {\mathit{\alpha }}_{\mathrm{T}}$, whereas cases D1 and D2 correspond to ${\mathit{\alpha }}_{\mathrm{T}}\ll \mathit{\alpha }⩽\mathrm{1}$. Case C ensures the transition between cases B and D1, when α≈α_T.

2.5 Asymptotic analysis

The next step is to introduce multiple-scale coordinates (Bensoussan et al., 1978; Sanchez-Palencia, 1980; Auriault, 1991). The two characteristic lengths L and l introduce two dimensionless space variables, ${\mathit{x}}^{*}=\mathbf{X}/L$ and ${\mathit{y}}^{*}=\mathbf{X}/l$, where X is the physical space variable. The macroscopic (or slow) dimensionless space variable x^* is related to the microscopic (or fast) dimensionless space variable y^* by ${\mathit{x}}^{*}=\mathit{\epsilon }{\mathit{y}}^{*}$. When l is used as the characteristic length, the dimensionless derivative operator grad^* becomes (${\mathbf{grad}}_{y^{*}}+\mathit{\epsilon }{\mathbf{grad}}_{x^{*}})$, where the subscripts ${}_{x^{*}}$ and ${}_{y^{*}}$ denote the derivatives with respect to the variables x^* and y^*, respectively. Following the multiple-scale expansion technique (Bensoussan et al., 1978; Sanchez-Palencia, 1980; Auriault, 1991), the ice temperature $T_{\mathrm{i}}^{*}$, the air temperature $T_{\mathrm{a}}^{*}$, and the water vapor ${\mathit{\rho }}_{\mathrm{v}}^{*}$ are sought in the form of the asymptotic expansion of the powers of ε, as follows:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\mathit{\phi }}^{*}({\mathit{x}}^{*},{\mathit{y}}^{*},t)& ={\mathit{\phi }}^{*\left(\mathrm{0}\right)}({\mathit{x}}^{*},{\mathit{y}}^{*},t)+\mathit{\epsilon }{\mathit{\phi }}^{*\left(\mathrm{1}\right)}({\mathit{x}}^{*},{\mathit{y}}^{*},t)\\ & +{\mathit{\epsilon }}^{\mathrm{2}}{\mathit{\phi }}^{*\left(\mathrm{2}\right)}({\mathit{x}}^{*},{\mathit{y}}^{*},t)+\mathrm{\dots },\end{array}\end{array}$$

where ${\mathit{\phi }}^{*}=T_{\mathrm{i}}^{*},T_{\mathrm{a}}^{*},{\mathit{\rho }}_{\mathrm{v}}^{*}$, and the corresponding φ^*(i) are periodic functions of period Ω with respect to the space variable y^*. Substituting these expansions in the set (13)–(18) gives, through the identification of identical powers of ε, successive boundary value problems to be investigated. All the details concerning this asymptotic analysis are presented in the Supplement. The main results are summarized in the following section.

2.6 Macroscopic-equivalent descriptions

2.6.1 Case A

Case A corresponds to the model presented in Calonne et al. (2014b). According to the order of magnitude of the dimensionless numbers, and notably $\left[H\right]=\mathcal{O}\left({\mathit{\epsilon }}^{\mathrm{2}}\right),\left[W_{\mathrm{R}}\right]=\mathcal{O}\left({\mathit{\epsilon }}^{\mathrm{2}}\right)$, the asymptotic analysis presented in the Supplement (Sect. S1) shows that the heat transfer and the water vapor diffusion at the macroscopic scale are described by Eqs. (S.A.44) and (S.A.47). Returning dimensional variables, the macroscopic model is written as

$$\begin{array}{}\text{(21)}& {\displaystyle }\begin{array}{rl}(\mathit{\rho }C{)}^{\mathrm{eff}}{\displaystyle \frac{\partial T^{\left(\mathrm{0}\right)}}{\partial t}}& -\mathrm{div}\left({\mathit{k}}^{\mathrm{eff}}\mathbf{grad}{T}^{\left(\mathrm{0}\right)}\right)={\mathrm{SSA}}_{\mathrm{V}}{L}_{\mathrm{sg}}{w}_{\mathrm{n}}^{\left(\mathrm{0}\right)}\\ & =-L_{\mathrm{sg}}\dot{\mathit{\varphi }},\end{array}\text{(22)}& {\displaystyle }\mathit{\varphi }{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}}{\partial t}}-\mathrm{div}\left({\mathit{D}}^{\mathrm{eff}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}\right)=-{\mathrm{SSA}}_{\mathrm{V}}{\mathit{\rho }}_{\mathrm{i}}{w}_{\mathrm{n}}^{\left(\mathrm{0}\right)}={\mathit{\rho }}_{\mathrm{i}}\dot{\mathit{\varphi }},\end{array}$$

where $w_{\mathrm{n}}^{\left(\mathrm{0}\right)}$ is given by the Hertz–Knudsen equation (Eq. S.A.43) and the Clausius–Clapeyron law (Eq. S.A.42; see Sect. S1 in the Supplement),

$$\begin{array}{}\text{(23)}& {\displaystyle }{w}_{\mathrm{n}}^{\left(\mathrm{0}\right)}={\displaystyle \frac{\mathit{\alpha }}{{\mathit{\rho }}_{\mathrm{i}}}}{w}_{\mathrm{k}}\left[{\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}-{\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right)\right],\text{(24)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right)={\mathit{\rho }}_{\mathrm{vs}}^{\mathrm{ref}}\exp \left[{\displaystyle \frac{L_{\mathrm{sg}}m}{{\mathit{\rho }}_{\mathrm{i}}k}}\left({\displaystyle \frac{\mathrm{1}}{T^{\mathrm{ref}}}}-{\displaystyle \frac{\mathrm{1}}{T^{\left(\mathrm{0}\right)}}}\right)\right],\end{array}$$

and where ϕ is the porosity and $\dot{\mathit{\varphi }}$ its total time derivative. ${\mathrm{SSA}}_{\mathrm{V}}=\left|\mathrm{\Gamma }\right|/\left|\mathrm{\Omega }\right|$ is the specific surface area per unit volume, defined as the ice surface area over the snow volume (in m⁻¹). The SSA can also be defined per unit mass, with ${\mathrm{SSA}}_{\mathrm{V}}=\mathrm{SSA}\times {\mathit{\rho }}_{\mathrm{i}}$. (ρC)^eff is the effective thermal capacity (Eq. S.A.45), k^eff is the effective thermal conductivity tensor (Eq. S.A.46), and D^eff is the effective diffusion tensor (Eq. S.A.48). These effective properties are defined as

$$\begin{array}{}\text{(25)}& {\displaystyle }(\mathit{\rho }C{)}^{\mathrm{eff}}=(\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}_i{C}_i+\mathit{\varphi }{\mathit{\rho }}_a{C}_{\mathrm{a}},\text{(26)}& {\displaystyle }{\mathit{k}}^{\mathrm{eff}}={\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\left(\underset{{\mathrm{\Omega }}_{\mathrm{a}}}{\int }{k}_{\mathrm{a}}(\mathbf{grad}{\mathit{t}}_{\mathrm{a}}+\mathbf{I})\mathrm{d}\mathrm{\Omega }+\underset{{\mathrm{\Omega }}_{\mathrm{i}}}{\int }{k}_{\mathrm{i}}(\mathbf{grad}{\mathit{t}}_{\mathrm{i}}+\mathbf{I})\mathrm{d}\mathrm{\Omega }\right),\text{(27)}& {\displaystyle }{\mathit{D}}^{\mathrm{eff}}={\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\underset{{\mathrm{\Omega }}_a}{\int }{D}_{\mathrm{v}}(\mathbf{grad}{\mathit{g}}_{\mathrm{v}}+\mathbf{I})\mathrm{d}\mathrm{\Omega },\end{array}$$

where t_a and t_i are two periodic vectors. The solution of the following boundary value problem over the REV Eqs. (S.A.20)–(S.A.24) is as follows:

$$\begin{array}{}\text{(28)}& {\displaystyle }\mathrm{div}\left(k_{\mathrm{i}}\right(\mathbf{grad}{\mathit{t}}_{\mathrm{i}}+\mathbf{I}\left)\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{i}},\text{(29)}& {\displaystyle }\mathrm{div}\left(k_{\mathrm{a}}\right(\mathbf{grad}{\mathit{t}}_{\mathrm{a}}+\mathbf{I}\left)\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\text{(30)}& {\displaystyle }{\mathit{t}}_{\mathrm{i}}={\mathit{t}}_{\mathrm{a}}\mathrm{on}\mathrm{\Gamma },\text{(31)}& {\displaystyle }\left(k_{\mathrm{i}}\right(\mathbf{grad}{\mathit{t}}_{\mathrm{i}}+\mathbf{I})-k_{\mathrm{a}}(\mathbf{grad}{\mathit{t}}_{\mathrm{a}}+\mathbf{I}\left)\right)\cdot {\mathit{n}}_{\mathrm{i}}=\mathrm{0}\mathrm{on}\mathrm{\Gamma },\text{(32)}& {\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\underset{\mathrm{\Omega }}{\int }({\mathit{t}}_{\mathrm{a}}+{\mathit{t}}_{\mathrm{i}})\mathrm{d}\mathrm{\Omega }=\mathbf{0}.\end{array}$$

Moreover, g_v is a periodic vector solution of the following boundary value problem over the REV Eqs. (S.A.35)–(S.A.37):

$$\begin{array}{}\text{(33)}& {\displaystyle }\mathrm{div}\left(D_{\mathrm{v}}\right(\mathbf{grad}{\mathit{g}}_{\mathrm{v}}+\mathbf{I}\left)\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\text{(34)}& {\displaystyle }{D}_{\mathrm{v}}(\mathbf{grad}{\mathit{g}}_{\mathrm{v}}+\mathbf{I})\cdot {\mathit{n}}_{\mathrm{i}}=\mathrm{0}\mathrm{on}\mathrm{\Gamma },\text{(35)}& {\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\underset{{\mathrm{\Omega }}_{\mathrm{a}}}{\int }{\mathit{g}}_{\mathrm{v}}\mathrm{d}\mathrm{\Omega }=\mathbf{0}.\end{array}$$

In that case, the above macroscopic equivalent description shows that, at the first order, the heat and water vapor transfer are described by two equations which are coupled through a source term proportional to the Hertz–Knudsen equation (Eq. 23) and the Clausius–Clapeyron law (Eq. 24) but expressed with respect to the two macroscopic variables T⁽⁰⁾ and ${\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}$. These equations involve two effective parameters, namely the effective thermal conductivity ${\mathit{k}}^{\mathrm{eff}}={\mathit{k}}^{\mathrm{eff}}(k_{\mathrm{i}},k_{\mathrm{a}},\mathrm{microstructure})$ and the effective diffusion ${\mathit{D}}^{\mathrm{eff}}={\mathit{D}}^{\mathrm{eff}}(D_{\mathrm{v}},\mathrm{microstructure})$.

2.6.2 Case B

According to the order of magnitude of the dimensionless numbers in case B, the asymptotic analysis presented in the Supplement (Sect. S2) shows that the heat transfer and the water vapor diffusion at the macroscopic scale are described by Eqs. (S.B.29) and (S.B.44). Returning dimensional variables, the macroscopic model is written as

$$\begin{array}{}\text{(36)}& {\displaystyle }(\mathit{\rho }C{)}^{\mathrm{eff}}{\displaystyle \frac{\partial T^{\left(\mathrm{0}\right)}}{\partial t}}-\mathrm{div}({\mathit{k}}^{\mathrm{eff}}\mathbf{grad}{T}^{\left(\mathrm{0}\right)})=-L_{\mathrm{sg}}\dot{\mathit{\varphi }},\text{(37)}& {\displaystyle }\mathit{\varphi }{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}}{\partial t}}-\mathrm{div}\left({\mathit{D}}^{\mathrm{eff}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\right)={\mathit{\rho }}_{\mathrm{i}}\dot{\mathit{\varphi }},\end{array}$$

with

$$\begin{array}{}\text{(38)}& {\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}={\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right),\end{array}$$

where (ρC)^eff is the effective thermal capacity, k^eff is the effective thermal conductivity tensor, and D^eff is the effective diffusion tensor, as defined in case A. The above macroscopic equivalent description shows that, at the first order, the heat and vapor transfer are only driven by the temperature field, since the water vapor density ${\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}={\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right)$ is directly given by the Clausius–Clapeyron law (Eq. 24). Consequently, from Eq. (37), we have the following:

$$\begin{array}{}\text{(39)}& \begin{array}{rl}\dot{\mathit{\varphi }}& =-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{i}}}}\left(\mathrm{div}\left({\mathit{D}}^{\mathrm{eff}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\right(T^{\left(\mathrm{0}\right)}\left)\right)-\mathit{\varphi }{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}}{\partial t}}\right)\\ & =-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{i}}}}\left(\mathrm{div}\left(\mathit{\gamma }\right(T^{\left(\mathrm{0}\right)}\left){\mathit{D}}^{\mathrm{eff}}\mathbf{grad}{T}^{\left(\mathrm{0}\right)}\right)-\mathit{\varphi }\mathit{\gamma }\left(T^{\left(\mathrm{0}\right)}\right){\displaystyle \frac{\partial T^{\left(\mathrm{0}\right)}}{\partial t}}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(40)}& \begin{array}{rl}\mathit{\gamma }\left(T^{\left(\mathrm{0}\right)}\right)& ={\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right)}{\mathrm{d}{T}^{\left(\mathrm{0}\right)}}}\\ & ={\mathit{\rho }}_{\mathrm{vs}}^{\mathrm{ref}}{\displaystyle \frac{L_{\mathrm{sg}}m}{{\mathit{\rho }}_{\mathrm{i}}k}}{\displaystyle \frac{\mathrm{1}}{(T^{\left(\mathrm{0}\right)}{)}^{\mathrm{2}}}}\exp \left[{\displaystyle \frac{L_{\mathrm{sg}}m}{{\mathit{\rho }}_{\mathrm{i}}k}}\left({\displaystyle \frac{\mathrm{1}}{T^{\mathrm{ref}}}}-{\displaystyle \frac{\mathrm{1}}{T^{\left(\mathrm{0}\right)}}}\right)\right]\\ & ={\displaystyle \frac{L_{\mathrm{sg}}m}{{\mathit{\rho }}_{\mathrm{i}}k}}{\displaystyle \frac{\mathrm{1}}{(T^{\left(\mathrm{0}\right)}{)}^{\mathrm{2}}}}{\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right).\end{array}\end{array}$$

Taking this result into account, the macroscopic heat transfer equation (Eq. 36) is written as

$$\begin{array}{}\text{(41)}& \left((\mathit{\rho }C{)}^{\mathrm{eff}}+\mathit{\varphi }\mathit{\gamma }(T^{\left(\mathrm{0}\right)}){\displaystyle \frac{L_{\mathrm{sg}}}{{\mathit{\rho }}_{\mathrm{i}}}}\right){\displaystyle \frac{\partial T^{\left(\mathrm{0}\right)}}{\partial t}}-\mathrm{div}\left({\stackrel{\mathrm{̃}}{\mathit{k}}}^{\mathrm{B}}\mathbf{grad}{T}^{\left(\mathrm{0}\right)}\right)=\mathrm{0}.\end{array}$$

In this latter equation,

$$\begin{array}{}\text{(42)}& {\stackrel{\mathrm{̃}}{\mathit{k}}}^{\mathrm{B}}={\mathit{k}}^{\mathrm{eff}}+{\displaystyle \frac{\mathit{\gamma }\left(T^{\left(\mathrm{0}\right)}\right)L_{\mathrm{sg}}}{{\mathit{\rho }}_{\mathrm{i}}}}{\mathit{D}}^{\mathrm{eff}},\end{array}$$

which appears as an apparent thermal conductivity of the snow which depends non-linearly on the temperature through γ(T⁽⁰⁾). Our results show that this is valid if [W_R]=𝒪(ε), i.e., typically for α values ranging from around 2 × 10⁻⁷ to 2 × 10⁻⁵. Finally, let us remark that (i) model B can be also seen as a particular case of model A when ${\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}$ tends towards ${\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right)$ by increasing α, and (ii) the apparent thermal conductivity of the snow ${\stackrel{\mathrm{̃}}{\mathit{k}}}^{\mathrm{B}}$ can also be written as

$$\begin{array}{}\text{(43)}& {\stackrel{\mathrm{̃}}{\mathit{k}}}^{\mathrm{B}}={\mathit{k}}^{\mathrm{eff}}+{\displaystyle \frac{\mathit{\gamma }\left(T^{\left(\mathrm{0}\right)}\right)L_{\mathrm{sg}}{D}_{\mathrm{v}}}{{\mathit{\rho }}_{\mathrm{i}}}}{\displaystyle \frac{{\mathit{D}}^{\mathrm{eff}}}{D_{\mathrm{v}}}}={\mathit{k}}^{\mathrm{eff}}+k_{\mathrm{dif}}{\displaystyle \frac{{\mathit{D}}^{\mathrm{eff}}}{D_{\mathrm{v}}}},\end{array}$$

where $k_{\mathrm{dif}}=\mathit{\gamma }\left(T^{\left(\mathrm{0}\right)}\right)L_{\mathrm{sg}}{D}_{\mathrm{v}}/{\mathit{\rho }}_{\mathrm{i}}$ corresponds to an enhancement of the air thermal conductivity, as defined in Sect. 2.4. However, in that case, γ(T⁽⁰⁾) depends on the macroscopic temperature T⁽⁰⁾.

2.6.3 Case C

According to the order of magnitude of the dimensionless numbers in case C, the asymptotic analysis presented in the Supplement (Sect. S3) shows that the heat transfer and the water vapor diffusion at the macroscopic scale are described by Eqs. (S.C.42) and (S.C.47). Returning dimensional variables, the macroscopic model is written as

$$\begin{array}{}\text{(44)}& {\displaystyle }(\mathit{\rho }C{)}^{\mathrm{eff}}{\displaystyle \frac{\partial T^{\left(\mathrm{0}\right)}}{\partial t}}-\mathrm{div}({\mathit{k}}^{\mathrm{C}}\mathbf{grad}{T}^{\left(\mathrm{0}\right)})=-L_{\mathrm{sg}}\dot{\mathit{\varphi }},\text{(45)}& {\displaystyle }\mathit{\varphi }{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}}{\partial t}}-\mathrm{div}\left({\mathit{D}}^{\mathrm{C}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\right(T^{\left(\mathrm{0}\right)}\left)\right)={\mathit{\rho }}_{\mathrm{i}}\dot{\mathit{\varphi }},\end{array}$$

with

$$\begin{array}{}\text{(46)}& {\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}={\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right),\end{array}$$

where (ρC)^eff and k^C are the dimensionless effective thermal capacity and the effective dimensionless thermal conductivity, respectively, and defined as

$$\begin{array}{}\text{(47)}& {\displaystyle }(\mathit{\rho }C{)}^{\mathrm{eff}}=(\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}_{\mathrm{i}}{C}_i+\mathit{\varphi }{\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{a}},\text{(48)}& {\displaystyle }{\mathit{k}}^{\mathrm{C}}={\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\left(\underset{{\mathrm{\Omega }}_{\mathrm{a}}}{\int }{k}_{\mathrm{a}}(\mathbf{grad}{\mathit{s}}_{\mathrm{a}}+\mathbf{I})\mathrm{d}\mathrm{\Omega }+\underset{{\mathrm{\Omega }}_{\mathrm{i}}}{\int }{k}_{\mathrm{i}}(\mathbf{grad}{\mathit{s}}_{\mathrm{i}}+\mathbf{I})\mathrm{d}\mathrm{\Omega }\right),\end{array}$$

where D^C is the effective diffusion tensor defined as

$$\begin{array}{}\text{(49)}& {\mathit{D}}^{\mathrm{C}}={\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\underset{{\mathrm{\Omega }}_a}{\int }{D}_{\mathrm{v}}\left(\mathbf{grad}\right(\mathit{d}+{\mathit{s}}_{\mathrm{a}})+\mathbf{I})\mathrm{d}\mathrm{\Omega }.\end{array}$$

Moreover, s_i, s_a, and d are the periodic vector solutions of the following coupled boundary value problem in a compact form (see Sect. S3 and Eqs. S.C.30–S.C.37 in the Supplement):

$$\begin{array}{}\text{(50)}& {\displaystyle }\mathrm{div}\left(k_{\mathrm{i}}\right(\mathbf{grad}{\mathit{s}}_{\mathrm{i}}+\mathbf{I}\left)\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{i}},\text{(51)}& {\displaystyle }\mathrm{div}\left(k_{\mathrm{a}}\right(\mathbf{grad}{\mathit{s}}_{\mathrm{a}}+\mathbf{I}\left)\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\text{(52)}& {\displaystyle }{\mathit{s}}_{\mathrm{i}}={\mathit{s}}_{\mathrm{a}}\mathrm{on}\mathrm{\Gamma },\text{(53)}& {\displaystyle }\begin{array}{rl}\left(k_{\mathrm{i}}\right(\mathbf{grad}{\mathit{s}}_{\mathrm{i}}+\mathbf{I})& -k_{\mathrm{a}}(\mathbf{grad}{\mathit{s}}_{\mathrm{a}}+\mathbf{I}))\cdot {\mathit{n}}_{\mathrm{i}}\\ & ={\displaystyle \frac{L_{\mathrm{sg}}}{{\mathit{\rho }}_{\mathrm{i}}}}\mathit{\alpha }{w}_{\mathrm{k}}\mathit{\gamma }\left(T^{\left(\mathrm{0}\right)}\right)\mathit{d}\mathrm{on}\mathrm{\Gamma },\end{array}\text{(54)}& {\displaystyle }\mathrm{div}\left(D_{\mathrm{v}}\right(\mathbf{grad}(\mathit{d}+{\mathit{s}}_{\mathrm{a}})+\mathbf{I}\left)\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\text{(55)}& {\displaystyle }\left(\mathbf{grad}\right(\mathit{d}+{\mathit{s}}_{\mathrm{a}})+\mathbf{I})\cdot {\mathit{n}}_{\mathrm{i}}=\mathit{\alpha }{w}_{\mathrm{k}}\mathit{d}\mathrm{on}\mathrm{\Gamma },\end{array}$$

with

$$\begin{array}{}\text{(56)}& {\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\underset{\mathrm{\Omega }}{\int }({\mathit{s}}_{\mathrm{a}}+{\mathit{s}}_{\mathrm{i}})\mathrm{d}\mathrm{\Omega }=\mathbf{0},\text{(57)}& {\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\underset{\mathrm{\Gamma }}{\int }\mathit{d}\mathrm{d}\mathrm{\Gamma }=\mathbf{0}.\end{array}$$

The vectors s_i, s_a, and d depend on the values of α and the temperature (notably through γ(T⁽⁰⁾)). As for model B, the macroscopic heat transfer equation (Eq. 44) can also be written as

$$\begin{array}{}\text{(58)}& \left((\mathit{\rho }C{)}^{\mathrm{eff}}+\mathit{\varphi }\mathit{\gamma }(T^{\left(\mathrm{0}\right)}){\displaystyle \frac{L_{\mathrm{sg}}}{{\mathit{\rho }}_{\mathrm{i}}}}\right){\displaystyle \frac{\partial T^{\left(\mathrm{0}\right)}}{\partial t}}-\mathrm{div}\left({\stackrel{\mathrm{̃}}{\mathbf{k}}}^{\mathrm{C}}\mathbf{grad}{T}^{\left(\mathrm{0}\right)}\right)=\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(59)}& {\stackrel{\mathrm{̃}}{\mathbf{k}}}^{\mathrm{C}}={\mathit{k}}^{\mathrm{C}}+{\displaystyle \frac{\mathit{\gamma }\left(T^{\left(\mathrm{0}\right)}\right)L_{\mathrm{sg}}}{{\mathit{\rho }}_{\mathrm{i}}}}{\mathit{D}}^{\mathrm{C}}={\mathit{k}}^{\mathrm{C}}+k_{\mathrm{dif}}{\displaystyle \frac{{\mathit{D}}^{\mathrm{C}}}{D_{\mathrm{v}}}},\end{array}$$

which appears as an apparent thermal conductivity of the snow. In that case k^C and D^C both depend on k_a, k_i, k_dif, and α, and we have

$$\begin{array}{}\text{(60)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathbf{k}}}^{\mathrm{C}}& ={\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\left(\underset{{\mathrm{\Omega }}_{\mathrm{a}}}{\int }(k_{\mathrm{a}}+k_{\mathrm{dif}})\left(\mathbf{grad}\right({\mathit{s}}_{\mathrm{a}}+\mathit{d})+\mathbf{I})\mathrm{d}\mathrm{\Omega }\right.\\ & \left.+\underset{{\mathrm{\Omega }}_{\mathrm{i}}}{\int }{k}_{\mathrm{i}}(\mathbf{grad}{\mathit{s}}_{\mathrm{i}}+\mathbf{I})\mathrm{d}\mathrm{\Omega }\right).\end{array}\end{array}$$

This model C is valid for [W_R]=𝒪(1), i.e., for ${\mathit{\epsilon }}^{\mathrm{1}/\mathrm{2}}<\left[W_{\mathrm{R}}\right]<{\mathit{\epsilon }}^{-\mathrm{1}/\mathrm{2}}$ and thus for α values in the range $({\mathit{\epsilon }}^{\mathrm{1}/\mathrm{2}}{D}_{\mathrm{vc}}/(l{w}_{\mathrm{kc}}\left)\right)<\mathit{\alpha }<({\mathit{\epsilon }}^{-\mathrm{1}/\mathrm{2}}{D}_{\mathrm{vc}}/(l{w}_{\mathrm{kc}}\left)\right)$ or, typically, $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}<\mathit{\alpha }<\mathrm{4}\times {\mathrm{10}}^{-\mathrm{3}}$.

2.6.4 Cases D1 and D2

According to the order of magnitude of the dimensionless numbers in cases D1 and D2, the asymptotic analysis presented in the Supplement (Sect. S4) shows that these two cases lead to the same macroscopic description. Returning dimensional variables, the macroscopic model (Eqs. S.D1.41–S.D1.45 or Eqs. S.D2.41–S.D2.45) is written as

$$\begin{array}{}\text{(61)}& {\displaystyle }(\mathit{\rho }C{)}^{\mathrm{eff}}{\displaystyle \frac{\partial T^{\left(\mathrm{0}\right)}}{\partial t}}-\mathrm{div}({\mathit{k}}^{\mathrm{D}}\mathbf{grad}{T}^{\left(\mathrm{0}\right)})=-L_{\mathrm{sg}}\dot{\mathit{\varphi }},\text{(62)}& {\displaystyle }\mathit{\varphi }{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}}{\partial t}}-\mathrm{div}\left({\mathit{D}}^{\mathrm{D}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\right)={\mathit{\rho }}_{\mathrm{i}}\dot{\mathit{\varphi }},\text{(63)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}={\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right).\end{array}$$

(ρC)^eff is the classical dimensionless effective thermal capacity. The macroscopic thermal conductivity tensor k^D and the macroscopic diffusion tensor D^D are defined as

$$\begin{array}{}\text{(64)}& {\displaystyle }{\mathit{k}}^{\mathrm{D}}={\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\left(\underset{{\mathrm{\Omega }}_{\mathrm{a}}}{\int }{k}_{\mathrm{a}}(\mathbf{grad}{\mathit{r}}_{\mathrm{a}}+\mathbf{I})\mathrm{d}\mathrm{\Omega }+\underset{{\mathrm{\Omega }}_{\mathrm{i}}}{\int }{k}_{\mathrm{i}}(\mathbf{grad}{\mathit{r}}_{\mathrm{i}}+\mathbf{I})\mathrm{d}\mathrm{\Omega }\right),\text{(65)}& {\displaystyle }{\mathit{D}}^{\mathrm{D}}={\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\underset{{\mathrm{\Omega }}_a}{\int }{D}_{\mathrm{v}}(\mathbf{grad}{\mathit{r}}_{\mathrm{a}}+\mathbf{I})\mathrm{d}\mathrm{\Omega },\end{array}$$

where r_a and r_i are two periodic vectors. The solution of the boundary value problem over the REV Eqs. (S.D1.30)–(S.D1.34) is as follows:

$$\begin{array}{}\text{(66)}& {\displaystyle }\mathrm{div}\left(k_{\mathrm{i}}\right(\mathbf{grad}{\mathit{r}}_{\mathrm{i}}+\mathbf{I}\left)\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{i}},\text{(67)}& {\displaystyle }\mathrm{div}\left(\right(k_{\mathrm{a}}+k_{\mathrm{dif}}\left)\right(\mathbf{grad}{\mathit{r}}_{\mathrm{a}}+\mathbf{I}\left)\right)=\mathrm{0}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\text{(68)}& {\displaystyle }{\mathit{r}}_{\mathrm{i}}={\mathit{r}}_{\mathrm{a}}\mathrm{on}\mathrm{\Gamma },\text{(69)}& {\displaystyle }\left(k_{\mathrm{i}}\right(\mathbf{grad}{\mathit{r}}_{\mathrm{i}}+\mathbf{I})-(k_{\mathrm{a}}+k_{\mathrm{dif}}\left)\right(\mathbf{grad}{\mathit{r}}_{\mathrm{a}}+\mathbf{I}\left)\right)\cdot {\mathit{n}}_{\mathrm{i}}=\mathrm{0}\mathrm{on}\mathrm{\Gamma },\text{(70)}& {\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\underset{\mathrm{\Omega }}{\int }({\mathit{r}}_{\mathrm{a}}+{\mathit{r}}_{\mathrm{i}})\mathrm{d}\mathrm{\Omega }=\mathbf{0}.\end{array}$$

As for models B and C, the macroscopic heat transfer equation (Eq. 61) can be also written as

$$\begin{array}{}\text{(71)}& \left((\mathit{\rho }C{)}^{\mathrm{eff}}+\mathit{\varphi }\mathit{\gamma }(T^{\left(\mathrm{0}\right)}){\displaystyle \frac{L_{\mathrm{sg}}}{{\mathit{\rho }}_{\mathrm{i}}}}\right){\displaystyle \frac{\partial T^{\left(\mathrm{0}\right)}}{\partial t}}-\mathrm{div}\left({\stackrel{\mathrm{̃}}{\mathbf{k}}}^{\mathrm{D}}\mathbf{grad}{T}^{\left(\mathrm{0}\right)}\right)=\mathrm{0}.\end{array}$$

In this latter equation,

$$\begin{array}{}\text{(72)}& {\stackrel{\mathrm{̃}}{\mathbf{k}}}^{\mathrm{D}}={\mathit{k}}^{\mathrm{D}}+{\displaystyle \frac{\mathit{\gamma }\left(T^{\left(\mathrm{0}\right)}\right)L_{\mathrm{sg}}}{{\mathit{\rho }}_{\mathrm{i}}}}{\mathit{D}}^{\mathrm{D}}={\mathit{k}}^{\mathrm{D}}+k_{\mathrm{dif}}{\displaystyle \frac{{\mathit{D}}^{\mathrm{D}}}{D_{\mathrm{v}}}},\end{array}$$

which appears as an apparent thermal conductivity of the snow. In that case, k^D and D^D both depend on k_a, k_i, and k_dif, and we have

$$\begin{array}{}\text{(73)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathbf{k}}}^{\mathrm{D}}& ={\displaystyle \frac{\mathrm{1}}{\left|\mathrm{\Omega }\right|}}\left(\underset{{\mathrm{\Omega }}_{\mathrm{a}}}{\int }(k_{\mathrm{a}}+k_{\mathrm{dif}})(\mathbf{grad}{\mathit{r}}_{\mathrm{a}}+\mathbf{I})\mathrm{d}\mathrm{\Omega }\right.\\ & \left.+\underset{{\mathrm{\Omega }}_{\mathrm{i}}}{\int }{k}_{\mathrm{i}}(\mathbf{grad}{\mathit{r}}_{\mathrm{i}}+\mathbf{I})\mathrm{d}\mathrm{\Omega }\right).\end{array}\end{array}$$

This model D corresponds to the one derived by Moyne et al. (1988), assuming that ρ_v=ρ_vs(T) is on the interface at the microscopic scale and using the volume averaging method. This model is also similar to the one derived by Hansen and Foslien (2015), assuming that $\mathit{\alpha }\approx {\mathrm{10}}^{-\mathrm{2}}$. In that case, we show that this model is valid for $\left[W_{\mathrm{R}}\right]=\mathcal{O}\left({\mathit{\epsilon }}^{-\mathrm{1}}\right)$ or 𝒪(ε⁻²), i.e., typically for α values ranging from around 5 × 10⁻³ to 1.

2.7 Macroscopic equivalent descriptions – synthesis

Figure 3 presents a summary of the four macroscopic models of heat and vapor transport in dry snow derived above, together with their domain of validity according to the value of α. As already mentioned, model A is the one already derived in Calonne et al. (2014b), whereas model D is equivalent to the model derived by Moyne et al. (1988) and Hansen and Foslien (2015).

Figure 3 Definition of the three different macroscopic models and their domain of validity with respect to α. The value of α_T was estimated based on the characteristic values given in Table 1.

A first important outcome is that the hypothesis ρ_v=ρ_vs(T), which is often made, appears to be a good approximation for α values larger than 10⁻⁶, as seen for models B, C, and D. The asymptotic analysis shows that in the range of α [10⁻⁶, α_T], this approximation is of the order of 𝒪(ε) since ${\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}={\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right)$, i.e., $\mathit{\sigma }=({\mathit{\rho }}_{\mathrm{v}}-{\mathit{\rho }}_{\mathrm{vs}})/{\mathit{\rho }}_{\mathrm{vs}}\approx \mathcal{O}\left(\mathit{\epsilon }\right)$. In the range [α_T, 1], this approximation is of the order of 𝒪(ε²), since ${\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}={\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right)$ and ${\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{1}\right)}={\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{1}\right)}\left(T^{\left(\mathrm{1}\right)}\right)$, i.e., σ≈𝒪(ε²). This result implies that models B, C, and D can be written in the same form and be reduced to a simple heat transfer equation. This heat transfer equation involves an apparent thermal conductivity ${\stackrel{\mathrm{̃}}{k}}^{\mathit{\beta }}$, which differs from one model to another. In contrast, model A does not presume any assumption about the vapor saturation.

A second point concerns the relative role of the vapor diffusion and of the sublimation–condensation in the vapor transport, which directly impacts the model formulation. In models A and B, the water vapor transfer is mainly limited by the sublimation–deposition at the ice–air interfaces. Model A consists of two equations of temperature and vapor density coupled through source terms that are proportional to α. Classically, heat and vapor transport are driven by the effective thermal conductivity of snow k^eff and the effective vapor diffusivity of snow D^eff, respectively. Both properties depend on the intrinsic properties of ice and air (k_i,k_a, and D_v) and on the snow microstructure. Model B can be seen as a particular case of model A, assuming that ρ_v tends towards ρ_vs(T) at the macro-scale, leading to a simplification in the form of one heat transfer equation in which the apparent thermal conductivity ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}$ can be easily expressed with respect to k^eff (k_i, k_a, and microstructure) and D^eff (D_v, and microstructure).

In contrast, in model D, the water vapor transfer is mainly limited by the diffusion process at the micro-scale. In that case, the model consists of a single heat transfer equation in which α is not involved but instead driven by an apparent thermal conductivity ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$, which can be expressed with respect to the macroscopic thermal conductivity k^D (k_i, k_a, k_dif, and microstructure) and to the macroscopic diffusion D^D (D_v, k_i, k_a, k_dif, and microstructure), with the latter appearing as a thermodiffusion coefficient. Note that these parameters depend on the air thermal properties k_a+k_dif that are enhanced by the phase change through k_dif.

The transition between a diffusion-limited and sublimation–deposition-limited vapor transport is captured by model C. This transition appears around a transition value α_T estimated here at $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$, yet we recall that it can vary between 10⁻⁵ and 10⁻³, depending on the grain size and, to a lesser extent, on temperature. Model C is of the same form as models B and D, but the apparent thermal conductivity ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{C}}$ can be expressed with respect to the macroscopic thermal conductivity k^C (k_i, k_a, k_dif, α, and microstructure) and the macroscopic diffusion D^C (D_v, k_i, k_a, k_dif, α, and microstructure). Unlike the other models, both macroscopic parameters k^C and D^C depend on α. These parameters tend towards k^eff and D^eff when α tends towards 10⁻⁵, thus recovering model B. They also tend towards k^D and D^D when α tends towards 1, thus recovering model D. Consequently, in practice, models A and C are sufficient to describe the heat and vapor transfer in the whole range of α (see Sect. 3.3).

3 Application to analytical and numerical cases

In this section, two simple snow microstructures, a bilayer snowpack and an assemblage of spherical grains and pores, are first considered to illustrate the influence of the microstructure and of the parameters taken at the pore scale on the macroscopic parameters of models A, B, and D (Sect. 3.2 and 3.1).

Then, a simplified 2D snow microstructure is considered to evaluate the models by comparing simulation results obtained with the pore-scale description and with the macroscopic modeling (Sect. 3.3).

Figure 4 Illustration of the bilayer snowpack problem (a) at the macroscopic scale and (b) at the scale of a representative elementary volume (REV).

3.1 The bilayer snowpack: upper and lower bounds

As a first example, we consider the classical bilayer material problem, and the snowpack is seen as a succession of horizontal layers of pure air and of pure ice, as illustrated in Fig. 4. In this case, the macroscopic parameters arising in models A, B, and D can be analytically determined and constitute the upper and lower bounds of these parameters for any anisotropic snow microstructure. The boundary value problems in (Eqs. 33–35), (28–32), and (66–70) have been solved analytically on the REV (Fig. 4b) in Auriault et al. (2009). Taking into account those results and using Eqs. (27) and (26), we have, for models A and B,

$$\begin{array}{}\text{(74)}& {\displaystyle }{\mathit{D}}^{\mathrm{eff}}=\left(\begin{array}{cc}{D}_{\mathrm{11}}^{\mathrm{eff}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right)D_{\mathrm{11}}^{\mathrm{eff}}=\mathit{\varphi }{D}_{\mathrm{v}},\text{(75)}& {\displaystyle }\begin{array}{rl}{\mathit{k}}^{\mathrm{eff}}& =\left(\begin{array}{cc}{k}_{\mathrm{11}}^{\mathrm{eff}}& \mathrm{0}\\ \mathrm{0}& k_{\mathrm{22}}^{\mathrm{eff}}\end{array}\right)k_{\mathrm{11}}^{\mathrm{eff}}=\mathit{\varphi }{k}_{\mathrm{a}}+(\mathrm{1}-\mathit{\varphi })k_{\mathrm{i}},\\ & k_{\mathrm{22}}^{\mathrm{eff}}={\displaystyle \frac{k_{\mathrm{i}}{k}_{\mathrm{a}}}{(\mathrm{1}-\mathit{\varphi })k_{\mathrm{a}}+\mathit{\varphi }{k}_{\mathrm{i}}}}.\end{array}\end{array}$$

Thus, it follows that

$$\begin{array}{}\text{(76)}& {\stackrel{\mathrm{̃}}{k}}_{\mathrm{11}}^{\mathrm{B}}=\mathit{\varphi }(k_{\mathrm{a}}+k_{\mathrm{dif}})+(\mathrm{1}-\mathit{\varphi })k_{\mathrm{i}},{\stackrel{\mathrm{̃}}{k}}_{\mathrm{22}}^{\mathrm{B}}={\displaystyle \frac{k_{\mathrm{i}}{k}_{\mathrm{a}}}{(\mathrm{1}-\mathit{\varphi })k_{\mathrm{a}}+\mathit{\varphi }{k}_{\mathrm{i}}}}.\end{array}$$

These results imply that the macroscopic properties (D^eff,k^eff, and ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}$) of any anisotropic snow verify the following bounds:

$$\begin{array}{}\text{(77)}& \mathrm{0}⩽D^{\mathrm{eff}}⩽\mathit{\varphi }{D}_{\mathrm{v}},\end{array}$$

and

$$\begin{array}{}\text{(78)}& \begin{array}{rl}{\displaystyle \frac{k_{\mathrm{i}}{k}_{\mathrm{a}}}{(\mathrm{1}-\mathit{\varphi })k_{\mathrm{a}}+\mathit{\varphi }{k}_{\mathrm{i}}}}& ⩽k^{\mathrm{eff}}⩽\mathit{\varphi }{k}_{\mathrm{a}}+(\mathrm{1}-\mathit{\varphi })k_{\mathrm{i}},\\ & {\displaystyle \frac{k_{\mathrm{i}}{k}_{\mathrm{a}}}{(\mathrm{1}-\mathit{\varphi })k_{\mathrm{a}}+\mathit{\varphi }{k}_{\mathrm{i}}}}⩽{\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}⩽\mathit{\varphi }(k_{\mathrm{a}}+k_{\mathrm{dif}})\\ & +(\mathrm{1}-\mathit{\varphi })k_{\mathrm{i}}.\end{array}\end{array}$$

For model D, from Eqs. (65) and (64), we have

$$\begin{array}{}\text{(79)}& {\displaystyle }\begin{array}{rl}{\mathit{D}}^{\mathrm{D}}& =\left(\begin{array}{ll}{D}_{\mathrm{11}}^{\mathrm{D}}& \mathrm{0}\\ \mathrm{0}& D_{\mathrm{22}}^{\mathrm{D}}\end{array}\right)D_{\mathrm{11}}^{\mathrm{D}}=\mathit{\varphi }{D}_{\mathrm{v}},\\ & D_{\mathrm{22}}^{\mathrm{D}}=\mathit{\varphi }{D}_{\mathrm{v}}{\displaystyle \frac{k_{\mathrm{i}}}{(\mathrm{1}-\mathit{\varphi })(k_{\mathrm{a}}+k_{\mathrm{dif}})+\mathit{\varphi }{k}_{\mathrm{i}}}},\end{array}\text{(80)}& {\displaystyle }\begin{array}{rl}{\mathit{k}}^{\mathrm{D}}& =\left(\begin{array}{cc}{k}_{\mathrm{11}}^{\mathrm{D}}& \mathrm{0}\\ \mathrm{0}& k_{\mathrm{22}}^{\mathrm{D}}\end{array}\right)k_{\mathrm{11}}^{\mathrm{D}}=\mathit{\varphi }{k}_{\mathrm{a}}+(\mathrm{1}-\mathit{\varphi })k_{\mathrm{i}},\\ & k_{\mathrm{22}}^{\mathrm{D}}={\displaystyle \frac{k_{\mathrm{i}}(k_{\mathrm{a}}+(\mathrm{1}-\mathit{\varphi }\left)k_{\mathrm{dif}}\right)}{(\mathrm{1}-\mathit{\varphi })(k_{\mathrm{a}}+k_{\mathrm{dif}})+\mathit{\varphi }{k}_{\mathrm{i}}}}.\end{array}\end{array}$$

Thus, it follows that

$$\begin{array}{}\text{(81)}& \begin{array}{rl}{k}_{\mathrm{11}}^{\mathrm{D}}& =\mathit{\varphi }(k_{\mathrm{a}}+k_{\mathrm{dif}})+(\mathrm{1}-\mathit{\varphi })k_{\mathrm{i}},\\ & k_{\mathrm{22}}^{\mathrm{D}}={\displaystyle \frac{k_{\mathrm{i}}(k_{\mathrm{a}}+k_{\mathrm{dif}})}{(\mathrm{1}-\mathit{\varphi })(k_{\mathrm{a}}+k_{\mathrm{dif}})+\mathit{\varphi }{k}_{\mathrm{i}}}}.\end{array}\end{array}$$

In that case, these results imply that the macroscopic properties (D^D,k^D, and ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$) of any anisotropic snow verify the following bounds:

$$\begin{array}{}\text{(82)}& \mathit{\varphi }{D}_{\mathrm{v}}⩽D^{\mathrm{D}}⩽\mathit{\varphi }{D}_{\mathrm{v}}{\displaystyle \frac{k_{\mathrm{i}}}{(\mathrm{1}-\mathit{\varphi })(k_{\mathrm{a}}+k_{\mathrm{dif}})+\mathit{\varphi }{k}_{\mathrm{i}}}},\end{array}$$

and

$$\begin{array}{}\text{(83)}& \begin{array}{rl}{\displaystyle \frac{k_{\mathrm{i}}(k_{\mathrm{a}}+(\mathrm{1}-\mathit{\varphi }\left)k_{\mathrm{dif}}\right)}{(\mathrm{1}-\mathit{\varphi })(k_{\mathrm{a}}+k_{\mathrm{dif}})+\mathit{\varphi }{k}_{\mathrm{i}}}}& ⩽k^{\mathrm{D}}⩽\mathit{\varphi }{k}_{\mathrm{a}}+(\mathrm{1}-\mathit{\varphi })k_{\mathrm{i}},\\ & {\displaystyle \frac{k_{\mathrm{i}}(k_{\mathrm{a}}+k_{\mathrm{dif}})}{(\mathrm{1}-\mathit{\varphi })(k_{\mathrm{a}}+k_{\mathrm{dif}})+\mathit{\varphi }{k}_{\mathrm{i}}}}⩽{\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}\\ & ⩽\mathit{\varphi }(k_{\mathrm{a}}+k_{\mathrm{dif}})+(\mathrm{1}-\mathit{\varphi })k_{\mathrm{i}}.\end{array}\end{array}$$

The above results show that, as already emphasized in Calonne et al. (2014b), Moyne et al. (1988), and Fourteau et al. (2021b), the bounds in Eqs. (77) and (82) of both the effective diffusion coefficients D^eff and D^D are always smaller than D_v and D^D>D^eff, whatsoever the α value. Moreover, according to the definition of D^eff (Eq. 74) and D^D (Eq. 79), if a vertical macroscopic temperature gradient is applied along e₂, model A (or B) will not predict any porosity variation along that direction because of the pore geometry. By contrast, model D, where the sublimation–deposition process is faster than diffusion, can predict mass transport along e₂, since $D_{\mathrm{22}}^{\mathrm{D}}\ne \mathrm{0}$, and thus a variation in the porosity along e₂.

Figure 5 Evolution of the SC estimates of the thermal conductivities $k_{\mathrm{SC}}^{\mathrm{eff}}$, ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{B}}$, and ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ with respect to porosity at four temperatures (a, c, e) and with respect to temperature for four porosities (b, d, f). The vertical dotted gray lines indicate the four temperature and porosity values considered.

3.2 Assemblage of spherical grains and pores: self-consistent estimates

The next analytical model is the self-consistent model (Bruggeman, 1935; Hill, 1965; Budiansky, 1965; Torquato, 2002). Previous works showed that self-consistent (SC) estimates provide good estimations of the macroscopic properties of heat and vapor transport in dry snow (Calonne et al., 2014b, a, 2019). In this model, the snow microstructure is considered a macroscopically isotropic material made of an assemblage of spherical inclusions of air or ice. Each type of inclusion is embedded in a homogeneous equivalent material, allowing us to account for the connectivity of both phases. The equivalent material corresponds to an infinite matrix for which the effective properties are the unknown to be calculated. The solution of the equations for an isolated inclusion then gives an implicit relation which can be solved for this effective property.

Figure 6 Evolution of the normalized SC estimates $D_{\mathrm{SC}}^{\mathrm{eff}}/D_{\mathrm{v}}$ and $D_{\mathrm{SC}}^{\mathrm{D}}/D_{\mathrm{v}}$ with respect to porosity at different temperatures (solid lines). Results of the numerical computations of $D_{\mathrm{SC}}^{\mathrm{D}}/D_{\mathrm{v}}$ at two temperatures from Fourteau et al. (2021a) are also shown (symbols).

Figure 7 Evolution of the thermal conductivity with temperature for k_i (Huang et al., 2013), k_a (Haynes, 2016), k_dif, and k_a+k_dif.

For model A, the SC estimate of the effective thermal conductivity of snow $k_{\mathrm{SC}}^{\mathrm{eff}}$ and of the effective diffusion coefficient $D_{\mathrm{SC}}^{\mathrm{eff}}$ verify the following implicit relation (Torquato, 2002):

$$\begin{array}{}\text{(84)}& {\displaystyle }\begin{array}{rl}{k}_{\mathrm{SC}}^{\mathrm{eff}}& ={\displaystyle \frac{\mathit{\beta }+\sqrt{{\mathit{\beta }}^{\mathrm{2}}+\mathrm{8}{k}_{\mathrm{i}}{k}_{\mathrm{a}}}}{\mathrm{4}}}\mathrm{with}\mathit{\beta }=k_{\mathrm{i}}\left(\mathrm{3}\right(\mathrm{1}-\mathit{\varphi })-\mathrm{1})\\ & +k_{\mathrm{a}}(\mathrm{3}\mathit{\varphi }-\mathrm{1}),\end{array}\text{(85)}& {\displaystyle }{D}_{\mathrm{SC}}^{\mathrm{eff}}=D_{\mathrm{v}}{\displaystyle \frac{(\mathrm{3}\mathit{\varphi }-\mathrm{1})}{\mathrm{2}}}.\end{array}$$

For model B, the SC estimate of the thermal conductivity ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{B}}$ is obtained by replacing the effective properties by their SC estimates in Eq. (41) and reads

$$\begin{array}{}\text{(86)}& {\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{B}}=k_{\mathrm{SC}}^{\mathrm{eff}}+k_{\mathrm{dif}}{\displaystyle \frac{D_{\mathrm{SC}}^{\mathrm{eff}}}{D_{\mathrm{v}}}}.\end{array}$$

Finally, for model D, the SC estimate of thermal conductivity ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ can be obtained by replacing k_a in Eq. (84) with k_a+k_dif as follows:

$$\begin{array}{}\text{(87)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}& ={\displaystyle \frac{\mathit{\beta }+\sqrt{{\mathit{\beta }}^{\mathrm{2}}+\mathrm{8}{k}_{\mathrm{i}}(k_{\mathrm{a}}+k_{\mathrm{dif}})}}{\mathrm{4}}}\\ & \mathrm{with}\mathit{\beta }=k_{\mathrm{i}}\left(\mathrm{3}\right(\mathrm{1}-\mathit{\varphi })-\mathrm{1})+(k_{\mathrm{a}}+k_{\mathrm{dif}})(\mathrm{3}\mathit{\varphi }-\mathrm{1}).\end{array}\end{array}$$

For the diffusion coefficient $D_{\mathrm{SC}}^{\mathrm{D}}$, it can be shown that (Auriault et al., 2009)

$$\begin{array}{}\text{(88)}& D_{\mathrm{SC}}^{\mathrm{D}}=\mathit{\varphi }{D}_{\mathrm{v}}{\displaystyle \frac{\mathrm{3}{\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}}{(k_{\mathrm{a}}+k_{\mathrm{dif}})+\mathrm{2}{\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}}}.\end{array}$$

The above SC estimates of the thermal conductivity and diffusion coefficient are presented in Figs. 5 and 6, and the impact of snow porosity and temperature is shown. To do so, we used the relationships of the thermal conductivity of ice k_i(T) and of air k_a(T), with the temperature from Huang et al. (2013) and Haynes (2016), respectively.

For thermal conductivity, the SC estimates $k_{\mathrm{SC}}^{\mathrm{eff}}$, ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{B}}$, and ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ at a given temperature are similar and follow the classical exponential evolution with snow porosity. Overall, estimates vary between about 0.06 W m⁻¹ K⁻¹ for porosity of 0.5 and 0.01 W m⁻¹ K⁻¹ for porosity of 1 (Fig. 5a, c, and e). More differences between the estimates can be seen for the normalized diffusion coefficient. The effective coefficient $D_{\mathrm{SC}}^{\mathrm{eff}}/D_{\mathrm{v}}$ is overall much smaller than $D_{\mathrm{SC}}^{\mathrm{D}}/D_{\mathrm{v}}$ and evolves linearly from 0.25 to 1 when porosity varies from 0.5 to 1 (Fig. 6). In contrast, $D_{\mathrm{SC}}^{\mathrm{D}}/D_{\mathrm{v}}$ shows a non-linear evolution from 0.7 to 1, with values close to 1 for porosity above 0.8. The non-linearity and the high values of $D_{\mathrm{SC}}^{\mathrm{D}}/D_{\mathrm{v}}$ come from the contribution of the heat conduction, through ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$, and of the latent heat, through k_dif. Finally, those estimates are in good agreement with the computed values from 3D images of snow in Fourteau et al. (2021a).

Next, we look at the impact of temperature on the properties. The impact is weaker than the one of porosity and more complex to understand, as dependencies are multiple. To help understand, we first break down the dependencies and show in Fig. 7 how the variables k_dif and k_a+k_dif and the thermal conductivity of ice k_i and air k_a evolve with temperature. When temperature increases from 210 to 273 K, the thermal conductivity of ice decreases, and the one of air slightly increase, with both evolutions being quasi-linear. Non-linearity is introduced with the parameter k_dif, which increases exponentially with temperature. Values for this parameter are small, even smaller than the air thermal conductivity, and are close to 0 W m⁻¹ K⁻¹ at −60 °C and reach 0.02 W m⁻¹ K⁻¹ at −3 °C. Finally, the term k_a+k_dif evolves in the same way as k_dif (non-linear), but the values are increased by k_a.

Keeping the above considerations in mind, the evolution of $k_{\mathrm{SC}}^{\mathrm{eff}}$, ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{B}}$, and ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ with temperature is presented in Fig. 5b, d, and f. For $k_{\mathrm{SC}}^{\mathrm{eff}}$, the SC estimates basically follow a monotonous decrease in the thermal conductivity with increasing temperature. This decrease is less pronounced for high porosity, and vice versa. These features directly result from the impact of the evolution of the ice and air thermal conductivity with temperature. The evolution of ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{B}}$ and ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ with temperature is more complex as the impact of k_dif superimposes. They show a non-linear evolution with temperature, with an evolution similar to $k_{\mathrm{SC}}^{\mathrm{eff}}$ for the lower temperatures and transitioning to an exponential increase for the higher temperatures, where the latter is driven by k_dif. We see that this non-linearity is even more important for ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ than for ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{B}}$, as k_dif appears in the definition of ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ several times. Finally, estimates of diffusion coefficient ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ show a slight influence of temperature through k_i, k_a, and k_dif and increase with decreasing temperature, in agreement with Fourteau et al. (2021a).

3.3 Numerical evaluation on a simplified 2D geometry

We perform a numerical evaluation of the obtained macroscopic models on a simplified 2D snow microstructure, as in Calonne et al. (2014b). We compare simulations of heat and water vapor transfer in snow obtained with the pore-scale description and with the macroscopic modeling.

3.3.1 Case study definition

Finite element numerical simulations were performed using the code from COMSOL Multiphysics software on a 2D vertical snow layer of 10 cm height and 0.5 cm width (Fig. 8). A constant temperature gradient of 100 or 500 K m⁻¹ was applied across the layer. Temperatures at the top T_top and at the bottom T_bottom are imposed, and T_bottom is kept at 273 K. For the water vapor conditions at the top and bottom, a null vapor flux is applied. Symmetry conditions are imposed on the lateral sides of the snow layer. Simulations were run in a steady state.

Figure 8 Illustration of the 2D geometry for the pore-scale modeling and the macroscopic equivalent modeling.

At the pore scale, the snow layer consists of 200 periodic cells of 0.5 × 0.5 mm²; each periodic cell (REV) is composed of an ice grain with a diameter of 0.3 mm surrounded by air, as shown in Fig. 8. The snow porosity is 0.71, which corresponds to a density of 266 kg m⁻³. The heat and the mass transfer is described by the set of Eqs. (1)–(12), where T_i, T_a, and ρ_v are the unknowns. This set of equations was numerically solved using the material parameter values presented in Table 1 and for different α values in the range of 10⁻¹⁰ to 1. For the sake of simplicity, the thermal conductivities k_i and k_a were kept constant in all the simulations. The chosen values correspond to a temperature of −10 °C (Table 1).

At the macroscopic scale, the snow layer is seen as a continuous equivalent medium. The heat and the mass transfer is described by the homogenized equations Eqs. (21)–(23) for model A, Eqs. (41) and (37) for model B, Eqs. (58) and (45) for model C, and Eqs. (71) and (62) for model D, where T⁽⁰⁾ and ${\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}$ are the macroscopic unknowns. These macroscopic descriptions involve different parameters and properties defined over the REV, which need to be provided. The porosity and the specific surface area SSA_V were equal to 0.71 and 3770 m⁻¹, respectively. The effective properties k^eff and D^eff were computed over the REV composed of a unique cell by solving the boundary value problems in Eqs. (33)–(35) and Eqs. (28)–(32), respectively. Given the symmetry of the REV, all the tensors involved in the macroscopic descriptions are isotropic. We found that k^eff=0.04243 W m⁻¹ K⁻¹ and $D^{\mathrm{eff}}=\mathrm{1.156}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹. The apparent thermal conductivity ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}$ was analytically deduced using Eq. (41). Its value depends on the temperature through the term k_dif(T). ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{C}}$ and D^C were computed over the REV by solving the boundary value problem in Eqs. (50)–(57) at different temperatures by varying the term k_a+k_dif(T) for α values in the range 10⁻⁶ to 1. Finally, ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$ and D^D were computed over the REV by solving the boundary value problem in Eqs. (66)–(70) at different temperatures in a similar way. In the considered temperature range, D^D is almost constant and equal to $\mathrm{2.01}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹.

Figure 9 Evolution of the thermal conductivities k^eff, ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}$, and ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$ with temperature. For ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$, the blue dots represent the numerical estimates over the REV, and the blue line is the fit.

Figure 9 presents the evolution of k^eff, ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}$, and ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$ with temperature. As expected, ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}$ and ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$ evolve non-linearly with T. To perform the simulations, the computed values of ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$ were fitted by the following relation: ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}=\mathrm{46.064}(T/\mathrm{273}{)}^{\mathrm{4}}-\mathrm{156.05}(T/\mathrm{273}{)}^{\mathrm{3}}+\mathrm{198.7}(T/\mathrm{273}{)}^{\mathrm{2}}-\mathrm{112.68}(T/\mathrm{273})+\mathrm{24.045}$ (blue line in Fig. 9). Figure 10 shows the evolution of the dimensionless diffusion coefficients $D^{\mathrm{eff}}/D_{\mathrm{v}}$, $D^{\mathrm{C}}/D_{\mathrm{v}}$, and $D^{\mathrm{D}}/D_{\mathrm{v}}$ and of the macroscopic thermal conductivities k^eff, k^C, and k^D with respect to α and for two temperatures of 270 and 250 K. As expected, only the parameters of model C ($D^{\mathrm{C}}/D_{\mathrm{v}}$ and k^C) vary with α and ensure a continuous transition between the parameters of model A ($D^{\mathrm{eff}}/D_{\mathrm{v}}$ and k^eff) and the ones of model D ($D^{\mathrm{D}}/D_{\mathrm{v}}$ and k^D). When fitting the numerical estimates, such a transition can be described by a simple function as follows:

$$\begin{array}{}\text{(89)}& {\displaystyle \frac{D^{\mathrm{C}}\left(\mathit{\alpha }\right)-D^{\mathrm{eff}}}{D^{\mathrm{D}}-D^{\mathrm{eff}}}}={\displaystyle \frac{k^{\mathrm{C}}\left(\mathit{\alpha }\right)-k^{\mathrm{eff}}}{k^{\mathrm{D}}-k^{\mathrm{eff}}}}={\displaystyle \frac{{\stackrel{\mathrm{̃}}{k}}^{\mathrm{C}}\left(\mathit{\alpha }\right)-{\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}}{{\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}-{\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}}}={\displaystyle \frac{A\mathit{\alpha }}{\mathrm{1}+A\mathit{\alpha }}},\end{array}$$

where A=1200 is a constant. This fit is shown with black lines in Fig. 10. The figure also includes the numerical estimations of the diffusion coefficient on 3D snow microstructures of different densities from Fourteau et al. (2021b), which are in good agreement with the proposed function in Eq. (89).

Figure 10 Evolution of the dimensionless diffusion coefficients $D^{\mathrm{eff}}/D_{\mathrm{v}}$, $D^{\mathrm{C}}/D_{\mathrm{v}}$, and $D^{\mathrm{D}}/D_{\mathrm{v}}$ and of the macroscopic thermal conductivities k^eff, k^C, and k^D with respect to α and for two temperatures (270 and 250 K). The black lines represent the proposed function in Eq. (89) to describe the parameters of model C. Numerical estimates of the diffusion coefficient on 3D snow microstructures of different densities from Fourteau et al. (2021b) are represented by the dots.

3.3.2 Comparison between pore-scale and macro-scale simulations

Results between pore-scale and macro-scale simulations are compared in terms of temperature, vapor density, and mass change rate. At the pore scale, the average values of each variable were taken over the cell and computed as follows:

$$\begin{array}{}\text{(90)}& \begin{array}{rl}\langle T\rangle & ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Omega }}}\left(\underset{{\mathrm{\Omega }}_{\mathrm{i}}}{\int }{T}_{\mathrm{i}}\mathrm{d}\mathrm{\Omega }+\underset{{\mathrm{\Omega }}_{\mathrm{a}}}{\int }{T}_{\mathrm{a}}\mathrm{d}\mathrm{\Omega }\right),\\ & \langle {\mathit{\rho }}_{\mathrm{v}}\rangle ={\displaystyle \frac{\mathrm{1}}{{\mathrm{\Omega }}_{\mathrm{a}}}}\underset{{\mathrm{\Omega }}_{\mathrm{a}}}{\int }{\mathit{\rho }}_{\mathrm{v}}\mathrm{d}\mathrm{\Omega }\langle {\mathit{\rho }}_{\mathrm{vs}}\left(T\right)\rangle ={\displaystyle \frac{\mathrm{1}}{{\mathrm{\Omega }}_{\mathrm{a}}}}\underset{{\mathrm{\Omega }}_{\mathrm{a}}}{\int }{\mathit{\rho }}_{\mathrm{vs}}\left(T_{\mathrm{a}}\right)\mathrm{d}\mathrm{\Omega },\\ & \langle \dot{\mathit{\varphi }}\rangle ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Omega }}}\underset{\mathrm{\Gamma }}{\int }{w}_{\mathrm{n}}\mathrm{d}\mathrm{\Gamma }.\end{array}\end{array}$$

Thus, we compare the vertical profiles of the pore-scale variables 〈T〉, 〈ρ_v〉, 〈ρ_vs〉, and $\langle \dot{\mathit{\varphi }}\rangle$ with the vertical profiles of the macroscopic variables T⁽⁰⁾, ${\mathit{\rho }}_{\mathrm{v}}^{\left(\mathrm{0}\right)}$, ${\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right)$, and $\dot{\mathit{\varphi }}$. As the obtained simulated temperature profiles were close to each other, to ease the comparison, we also use the temperature deviation ΔT, which represents the deviation of the simulated temperature profile from the linear temperature profile imposed by T_top and T_bottom, as illustrated in Fig. 11. In the same vein, we use the water vapor supersaturation, which is the difference between the simulated water vapor density and the saturated water vapor density ρ_v−ρ_vs(T). Figure 12 shows the vertical profiles of ΔT, of ρ_v−ρ_vs(T), and of $\dot{\mathit{\varphi }}$ from the pore-scale simulations (dots) and the macroscopic models (lines), considering a temperature gradient of 100 and 500 K m⁻¹. For the pore-scale simulations, values of α from 10⁻⁹ to 1 were used. For the macroscopic models, results are only shown in their domain of validity with respect to α. To further highlight the impact of α, Fig. 13 presents the evolution of T, ρ, and ρ_vs with α for the specific middle cell of the snow layer, which is here the 100th cell from the bottom ($x=l/\mathrm{2}$ and $y=\mathrm{100}l-l/\mathrm{2}$). Again, pore-scale simulations (dots) and the macroscopic models (line) are compared.

Figure 11 Simplified example of the transition from the x to Δx notation in a concave (red) and a convex (blue) case.

Figure 12 Vertical profiles of ΔT, ρ_v−ρ_vs(T), and $\dot{\mathit{\varphi }}$ from the pore-scale simulations (dots) and from the macroscopic models A (grey lines), B (orange lines), C (magenta lines), and D (blue lines), considering a temperature gradient of 100 and 500 K m⁻¹ and for different values of α. ΔT represents the deviation of the temperature profile from a linear temperature profile. Predictions of ρ_v−ρ_vs(T) from models C and D are not shown in panels (c) and (d), as they are superimposed with the pore-scale simulation results for α = 1 (dotted yellow lines).

Figure 13 Temperature and water vapor density in the middle of the snow layer as a function of α, obtained from the pore-scale simulations (dots) and from the macroscopic models A (grey lines), B (orange lines), C (magenta lines), and D (blue lines) at 100 and 500 K m⁻¹. The models are only shown for the α values within their domain of validity. Values of saturation water vapor density ρ_vs from the pore-scale simulations and from model A are also presented.

We describe first the main features observed in the pore-scale simulations. All the variables show an impact of the α value. The temperature deviation ΔT is overall mainly positive (Fig. 12a and b), which reflects the presence of a heat source by non-conductive processes such as latent heat from deposition. This temperature deviation increases with α and with the temperature gradient. This is also reflected in the temperature of the middle cell that overall increases with increasing α (Fig. 13a and b). This increase is not uniform, and two plateaus are observed – one between ${\mathrm{10}}^{-\mathrm{6}}⩽\mathit{\alpha }⩽{\mathit{\alpha }}_{\mathrm{T}}$ and the other one between ${\mathrm{10}}^{-\mathrm{1}}⩽\mathit{\alpha }⩽\mathrm{1}$. The largest ΔT value is reached in the center of the snow layer and is around 0.4 K at 100 K m⁻¹ and 4 K at 500 K m⁻¹. Looking at the lower part of the layer, negative ΔT values can be found for $\mathit{\alpha }⩽{\mathit{\alpha }}_{\mathrm{T}}\sim \mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ and indicate a heat loss by non-conductive processes such as latent heat from sublimation. This feature vanishes for the large temperature gradient. In terms of water vapor supersaturation ρ_v−ρ_vs(T), we observe positive values (over-saturation) in the upper part of the snow layer, values close to zero in the central part (at saturation), and negative values (under-saturation) in the lower part (Fig. 12c, d). The largest over-saturation and under-saturation values are shown for low α values when phase changes are very limited. With increasing α values, values close to saturation get predominant, and the over-saturation and under-saturation zones become localized near the top and bottom, respectively. This is confirmed in Fig. 13c and d, where ρ_v≈ρ_vs(T) for ${\mathrm{10}}^{-\mathrm{6}}⩽\mathit{\alpha }⩽\mathrm{1}$. Similar to the temperature, two plateaus are shown where ρ_v−ρ_vs(T) evolve little with α. All of these results are consistent with our theoretical analysis presented in Sect. 2. Last, the vertical profiles of $\dot{\mathit{\varphi }}$ are consistent with the ones of supersaturation, showing deposition in the upper part where the porosity decreases and sublimation in the lower part where the porosity increases (Fig. 12e and f). As α increases, those transitions become sharper and sharper, like a front. For ${\mathit{\alpha }}_{\mathrm{T}}⩽\mathit{\alpha }⩽\mathrm{1}$, most values become negative, indicating overall deposition in the snow layer. A sublimation zone is still visible at the bottom of the snow layer, but its thickness is typically of the order of a few REV or smaller. Finally, as the difference ρ_v−ρ_vs is directly related to the interface growth velocity w_n (see Eq. 10), and as it could be useful to compare it with experimental estimates (e.g., Flin and Brzoska, 2008; Brzoska et al., 2008; Pinzer et al., 2012; Libbrecht and Rickerby, 2013), we provide below the mean values of w_n computed for the bottom and middle cell of α = 10⁻⁶. For 100 K m⁻¹, a value of 5.9 × 10⁻¹³ m s⁻¹ and of −2.7 × 10⁻¹¹ m s⁻¹ is found in the middle and bottom cell, respectively. For 500 K m⁻¹, a value of 4.5 × 10⁻¹² m s⁻¹ and of −1.1 × 10⁻¹⁰ m s⁻¹ is found in the middle and bottom cell, respectively.

Next we compare the different macroscopic models to the pore-scale simulations. In both Figs. 12 and Fig. 13, the comparison shows different behaviors, depending on α. For $\mathit{\alpha }⩽{\mathrm{10}}^{-\mathrm{5}}$, model A reproduces precisely all the features shown at the pore scale. Models B and D are independent of α and provide one estimate of the temperature, and thus the water vapor density, for all the α values in their domain of validity. These estimates are only able to reproduce the plateau values observed in the pore-scale simulations, i.e., temperatures for ${\mathrm{10}}^{-\mathrm{7}}⩽\mathit{\alpha }⩽{\mathrm{10}}^{-\mathrm{5}}$ for model B and ${\mathrm{10}}^{-\mathrm{2}}⩽\mathit{\alpha }⩽\mathrm{1}$ for model D (Fig. 13). Model C, which depends on α, allows us to reproduce the main features shown at the pore scale for α values in the range ${\mathrm{10}}^{-\mathrm{5}}⩽\mathit{\alpha }⩽\mathrm{1}$. This model predicts ΔT values higher than the pore-scale simulations for a given α. Moreover, models B, C, and D predict only the deposition in snow, with negative $\dot{\mathit{\varphi }}$ values throughout the layer (Fig. 12e and f). They do not capture the sublimation front at the bottom of the snow layer, in contrast to model A. The mass balance between sublimation–deposition over the whole snow layer is not well satisfied; the Dirichlet boundary condition imposed on the temperature field at the bottom and the top of the snow cannot ensure that the water vapor flux is null at the same time since ρ_v=ρ_vs(T). This limitation can explain the slight differences that we can observe between model C predictions and the pore-scale simulations. In order to overcome this limit, specific boundary conditions should be introduced to allow the description of mass variations near the interfaces.

4 Application to experimental data

This section presents the evaluation of the macroscopic models A, B, and D, based on observations of natural snow evolution from three cold-laboratory TGM experiments. We first introduce the experimental data (Sect. 4.1); then we define the estimates to be taken for the input parameters of the models (Sect. 4.2); and, finally, we present the simulation results with the models and their comparison with the experiments (Sect. 4.3).

Table 2 Overview of the experimental settings used in the simulations.

4.1 Experimental datasets

We used the datasets provided by Bouvet et al. (2023), consisting of two experiments, referred to as Bouvet A and Bouvet B in their paper and hereafter, and the data from Kamata and Sato (2007), referred to as Kamata. These experiments provide the required data to evaluate our models, namely the time series of the vertical profiles of temperature and density of a snow layer evolving under a temperature gradient in a controlled environment. The main characteristics of the three experiments are summarized in Table 2. Bouvet A is a TGM experiment on a 13.5 cm high snow layer for which a temperature gradient (TG) of 93 K m⁻¹ was applied during 20 d. X-ray tomography was done at regular time intervals resulting in 9 large 3D images of the whole vertical dimension of the snow layer at a resolution of 21 µm and 17 small 3D images of the top or bottom part of the layer at a resolution of 8 µm. For the large images, the first few millimeters at the base of the layer are lacking, due to the snow sampling procedure, so no data are available for this area. This experiment also includes monitoring of the temperature profile of the snow layer, measured using seven Pt100 sensors. Bouvet B is a TGM experiment on a 7.7 cm high snow layer for which a TG of 103 K m⁻¹ was applied during 28 d. Four tomography images of the first lower 4.2 cm of the snow layer are provided at a resolution of 10 µm. For both experiments, Bouvet A and Bouvet B, the vertical profiles of the snow density computed from the 3D images are provided, and a vertical mass redistribution can be analyzed. Finally, the Kamata experiment is a TGM experiment on a snow layer of 10 cm height for which an extreme TG of 530 K m⁻¹ was applied during 5.5 d (133 h). The vertical mass redistribution was estimated by measuring snow density for four sections of the snow layer. For that, the snow layer was separated into vertical compartments of about 2.5 cm height each, using horizontal nylon meshes, which enables water vapor to get through. Each compartment was weighed at the initial and final stage of the experiment. In addition, the temperature was recorded at six vertical locations.

4.2 Effective properties and parameters

Next we study the estimates of the effective properties and other input parameters required to run models A, B, and D. For the sake of simplicity, model C is not systematically shown. For each model, these properties are computed from the 3D images of snow of the experiment Bouvet A and Bouvet B. Those values are then compared to different parameterizations from the literature or fitted regressions, and we select the more suited ones to be used later in the models.

4.2.1 Model A

Model A involves three effective parameters that are the effective thermal conductivity k^eff, the effective vapor diffusivity D^eff, and the SSA (Sect. 2.6.1). These parameters were estimated in the case of Bouvet A and Bouvet B by numerical computations on the 3D tomographic images available. The SSA was computed per unit of mass based on the voxel projection approach (Flin et al., 2011; Dumont et al., 2021). k^eff and D^eff were computed with the GeoDict software (Thoemen et al., 2008) by solving the boundary value problems in Eqs. (28)–(32) and Eqs. (33)–(35) on the 3D images and by applying periodic boundary conditions on the external boundaries, as described in Calonne et al. (2011, 2014b). Values of k_i and k_a at −10 °C were used for the computation of k^eff. The obtained 3D tensors of both properties show negligible non-diagonal terms. In the following, we refer to k^eff and D^eff as the average of the diagonal terms of the tensors.

Figure 14 Average values of effective conductivity and normalized effective diffusivity as a function of density and SSA as a function of time, as computed from the tomography images of Bouvet A (symbols; a–c) and Bouvet B (symbols; d–f). The error bars represent the standard deviation of the parameter along the image height. Comparison with the SC model, classical parameterizations, and fits are shown (solid lines).

Figure 14 presents the results of the image-based computations of k^eff, D^eff, and SSA for the experiment of Bouvet A and Bouvet B. To compare them, we show the estimates of k^eff and D^eff by the SC model presented in Sect. 3.2, the density-based parameterizations of k^eff from Calonne et al. (2011) and Riche and Schneebeli (2013), and a fitted regression of the SSA values as a function of time, referred as SSA_Fit(t), based on a logarithmic function as formerly proposed by Legagneux et al. (2004). For thermal conductivity, the parameterization of Calonne et al. (2011) is in good agreement with the image-based computations in both experiments, whereas the parameterization of Riche and Schneebeli (2013), which specifically describes the case of depth hoar, predicts slightly larger values. The SC model largely underestimates the values, which are about 2 to 4 times smaller than the image-based computations. For the vapor diffusion coefficient, the SC model provides overall fair estimates, which are slightly overestimated, especially towards the end of the experiments, as reported for depth hoar and faceted crystals in Calonne et al. (2014b). For SSA, the fit reproduces the SSA evolution well for the experiment of Bouvet A. For Bouvet B, the SSA does not follow the classic exponential decrease, but it increases after 7 d and until the end of the experiment; this increase is specific to hard-depth hoar formation (Bouvet et al., 2023). This feature is not predicted by the applied fit, yet it provides fair estimates of the SSA values.

Given the above considerations, we selected two sets of parameters to simulate Bouvet A, Bouvet B, and Kamata with model A, which are summarized in Table 3. In the set called Calonne, k^eff is estimated with the parameterization of Calonne et al. (2011). In the second set called SC, k^eff is given by the self-consistent estimates. In both sets, D^eff is estimated with the SC model and the SSA with the logarithmic fit, which is specific for Bouvet A and Bouvet B. In the case of the Kamata experiment, we cannot test the proposed estimates of k^eff, D^eff, and SSA against reference data as such data are not available from Kamata and Sato (2007). The SSA evolution was reproduced based on the logarithmic fit from Bouvet A, as both experiments are the closest in terms of initial snow type, grain size, and density. In the model evaluation that follows (Sect. 4.3), the Calonne set is the one used by default to evaluate model A. Results with the SC set are also presented to illustrate an alternative choice of parameters, which, although less accurate, allows for consistent and analytically based estimates for all properties.

Table 3 Summary of the effective parameters used in the simulations.

4.2.2 Models B and D

Models B and D only involve the apparent thermal conductivities of snow ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}$ and ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$, respectively. Estimating ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}$ comes down to estimating k^eff and D^eff, as it is defined as ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}={\mathit{k}}^{\mathrm{eff}}+k_{\mathrm{dif}}{\mathit{D}}^{\mathrm{eff}}/D_{\mathrm{v}}$. For that reason, we use the same estimates of k^eff and D^eff selected for model A as described above. So two sets of input parameters were used for model B, namely the Calonne set, from which the model's performances are evaluated, and the alternative SC set (Table 3). Figure 15 presents the evolution of ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Calonne}}^{\mathrm{B}}$ and ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{B}}$ with temperature for each experiment, taking the mean snow density of the experiments (see Table 2). Both estimates show similar trends but, as in model A, the SC estimate predicts lower values than when using the parameterization of Calonne et al. (2011).

Figure 15 The thermal conductivity estimates for model B (${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{B}}$ and ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Calonne}}^{\mathrm{B}}$) and for model D (${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ and ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Fit}}^{\mathrm{D}}$) are presented as a function of temperatures (solid and dotted lines). The parameters are presented for the Bouvet A, Bouvet B, and Kamata experiments, and the green areas represent the temperature ranges of each experiment. The computed values on 3D images used to derive ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Fit}}^{\mathrm{D}}$ are shown by blue dots. For Bouvet A, ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Fit}}^{\mathrm{D}}\left(T\right)=\mathrm{5.4485}\times {\mathrm{10}}^{-\mathrm{9}}\times T^{\mathrm{4}}-\mathrm{4.8119}\times {\mathrm{10}}^{-\mathrm{6}}\times T^{\mathrm{3}}+\mathrm{1.5965}\times {\mathrm{10}}^{-\mathrm{3}}\times T^{\mathrm{2}}-\mathrm{2.3581}\times {\mathrm{10}}^{-\mathrm{1}}\times T+\mathrm{1.3195}$. For Bouvet B, ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Fit}}^{\mathrm{D}}\left(T\right)=\mathrm{6.0212}\times {\mathrm{10}}^{-\mathrm{9}}\times T^{\mathrm{4}}-\mathrm{5.2974}\times {\mathrm{10}}^{-\mathrm{6}}\times T^{\mathrm{3}}+\mathrm{1.7523}\times {\mathrm{10}}^{-\mathrm{3}}\times T^{\mathrm{2}}-\mathrm{2.5868}\times {\mathrm{10}}^{-\mathrm{1}}\times T+\mathrm{14.6338}$. For Kamata, ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Fit}}^{\mathrm{D}}\left(T\right)=\mathrm{5.1386}\times {\mathrm{10}}^{-\mathrm{9}}\times T^{\mathrm{4}}-\mathrm{4.5612}\times {\mathrm{10}}^{-\mathrm{6}}\times T^{\mathrm{3}}+\mathrm{1.5206}\times {\mathrm{10}}^{-\mathrm{3}}\times T^{\mathrm{2}}-\mathrm{2.2553}\times {\mathrm{10}}^{-\mathrm{1}}\times T+\mathrm{12.6279}$.

For model D, ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$ was computed from the 3D images from Bouvet A and Bouvet B by solving the boundary problem in Eqs. (66)–(70) with the GeoDict software. As above, only diagonal terms of the tensor were considered, and ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$ refers to the average value of the diagonal terms. Here, computations were performed on only one REV from each experiment and for 10 temperatures ranging from 210 to 273 K. We selected the image at 14 d for Bouvet A (cropped between 5.8 and 6.7 cm height) and the image at 7 d for Bouvet B (cropped between 1.7 and 2.5 cm height). To be able to estimate ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$ for the Kamata experiment, we took a 3D image of snow showing similar characteristics, namely a sample of depth hoar with a density of 165 kg m⁻¹ from Fourteau et al. (2021a). Results of the image-based computations are presented in Fig. 15, as well as estimates from the SC model ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ presented in Sect. 3.2. For each experiment, a fit was performed on the computed data and is referred to as ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Fit}}^{\mathrm{D}}$ in the following. Again, the SC estimates for model D capture the trend but largely underestimate the values. In what follows, simulations were performed using the fitted values ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Fit}}^{\mathrm{D}}$, referred to as the Fit set in Table 3, upon which the evaluation of model D is based. As for the other models, simulations with the SC estimates ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$ are also presented for the sake of comparison as they allow for independent and consistent estimates.

4.3 Comparisons between models and experiments

In this section, we compare simulations from models A, B, and D with the measurements from the three experiments Bouvet A, Bouvet B, and Kamata. The simulations were performed with the software COMSOL Multiphysics by resolving the homogenized equations on a 1D geometry that corresponds to the snow layer of the experiments. We used Eqs. (21)–(22) for model A, Eqs. (41) and (37) for model B, and Eqs. (71) and (62) for model D. For model A, the boundary conditions in temperature are the top and bottom imposed temperatures of the experiments. In terms of vapor density, the conditions correspond to zero flux at the top and bottom. Additionally, the source term is forced to zero in the simulation nodes where a density of zero is reached. For models B and D, the boundary conditions are the imposed temperatures. The models were run using the sets of input parameters described in Table 3 and considering the experimental conditions summarized in Table 2. Comparisons between measurements and simulations are performed based on temperature and mass change variables.

4.3.1 Temperature

Figure 16 presents the measured and simulated vertical profiles of ΔT with models A, B, and D for the three experiments, taking different α values from 10⁻⁹ to 10⁻⁴ for model A. The ΔT values analyzed here are the ones computed at the beginning of the experiments when the temperature gradient is well established but before the formation of the air gap. Overall, profiles of ΔT have similar shapes as the ones simulated on the simplified 2D microstructure (Sect. 3.3), describing right-headed curves indicating that processes apart from heat conduction, such as phase change, occur and result in a heat source in the snow layer. In the center part of the layer, a maximum deviation of 1.15 K was measured in the Bouvet A experiment and of 5.9 K in the Kamata experiment. The negative ΔT value in the upper part of the layer in Bouvet A is attributed to the temperature sensor uncertainty (Bouvet et al., 2023).

Figure 16 Vertical steady-state profiles of ΔT simulated with model A, with α ranging from 10⁻⁹ to 10⁻⁴, and with models B and D for the Bouvet A, Bouvet B, and Kamata experiments. Simulations using the Calonne set (solid lines) are shown for model A, and simulations using both Calonne (solid lines) and SC sets (dashed lines) are shown for models B and D. The experimental profiles are shown with dashed black lines for Bouvet A and Kamata.

Looking at the models, the main observation is that they all underestimate ΔT. In more detail, model A predicts negative ΔT in the lower part of the snowpack and is positive otherwise, reflecting a heat sink attributed to more sublimation in the lower part and a heat source attributed to more deposition in the rest of the layer. With increasing α, the positive values of ΔT increase, and the negative ones tends to vanish, so that the shape of the simulated curve becomes closer to the experimental one. The maximum ΔT predicted by model A is reached for the highest α = 10⁻⁴ and is of 0.11 K for Bouvet A and of 0.49 K for Kamata, which corresponds only to 10 % and 8 % of the experimental value, respectively.

Models B and D show a unique ΔT profile valid over their domain of validity, namely ${\mathrm{10}}^{-\mathrm{6}}⩽\mathit{\alpha }⩽{\mathit{\alpha }}_{\mathrm{T}}$ and ${\mathit{\alpha }}_{\mathrm{T}}⩽\mathit{\alpha }⩽\mathrm{1}$, respectively. The profile shape is in agreement with the measurements, showing only positive values throughout the layer. ΔT values of model B correspond to the upper limit of model A. Model D is the closest to the experimental data. Still, values are largely underestimated and reach at most 0.29 K (25 % of the experimental data) for Bouvet A and 1.4 K (23 % of the experimental data) for Kamata. In both models B and D, slightly better results are found when using the SC set of input parameters, even though it corresponds to underestimated estimates, as seen in Sect. 4.2. This better agreement with the SC estimate is somehow artificial and comes from the fact that the lower values of ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{B}}$ and ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{SC}}^{\mathrm{D}}$, compared to ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Calonne}}^{\mathrm{B}}$ and ${\stackrel{\mathrm{̃}}{k}}_{\mathrm{Fit}}^{\mathrm{D}}$, respectively, lead to reduced overall heat conduction through snow and allow for higher ΔT, as well as the fact that the SC estimates allow for a slightly higher sensitivity (steeper slope) of the thermal conductivity to temperature in the temperature range of the considered experiments (see the green areas in Fig. 15). A final interesting point is the strong impact of the density on ΔT, which can be seen by comparing simulations of Bouvet A and Bouvet B for which temperature gradients were very close, but snow density was 210 and 287 kg m⁻³, respectively. For the same temperature gradient, the higher the snow density, the higher the heat conduction through snow and the lower the ΔT. For example, in lighter snow (Bouvet A), a maximum ΔT of 0.29 K is predicted by model D against 0.06 K for the denser snow (Bouvet B).

4.3.2 Mass change

Next, we evaluate the models regarding mass changes across the vertical dimension of the snow layer. We look at the vertical density profile of snow, as well as the height of the air gap formed at the base of the layer at the end of the experiments, caused by an upward mass transfer during the TGM. For model A, the air gap height is defined as the highest height value at which the density is zero. For models B, C, and D, the vertical profile of density cannot be evaluated because they only predict deposition and thus density increase, due to the boundary condition issue already described in Sect. 3.3. Still, to allow for a comparison with measurements, we derived a rough estimate of the air gap by considering that all the mass gained in the snow layer over the whole experiment duration is balanced by a mass loss localized at the very bottom of the snow layer, leading to a sharp air gap described as

$$\begin{array}{}\text{(91)}& \begin{array}{rl}{h}_{\text{air gap}}& =H{\displaystyle \frac{\stackrel{\mathrm{‾}}{\dot{\mathit{\varphi }}}{t}_{\exp }}{(\stackrel{\mathrm{‾}}{\dot{\mathit{\varphi }}}{t}_{\exp }+{\mathit{\varphi }}_{\mathrm{init}}-\mathrm{1})}},\\ & \mathrm{with}\stackrel{\mathrm{‾}}{\dot{\mathit{\varphi }}}={\displaystyle \frac{\mathrm{1}}{H}}\underset{\mathrm{0}}{\overset{H}{\int }}\dot{\mathit{\varphi }}\left(z\right)\mathrm{d}z\\ & =-{\displaystyle \frac{\mathrm{1}}{H}}\underset{\mathrm{0}}{\overset{H}{\int }}{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{i}}}}{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{\partial D{\mathit{\rho }}_{\mathrm{vs}}^{\left(\mathrm{0}\right)}\left(T\right)}{\partial z}}\right)\mathrm{d}z,\end{array}\end{array}$$

with h_{air gap} as the height of the air gap (m), H as the total height of the snow layer (m), ϕ_init as the initial porosity (–), and t_exp as the total duration of the experiment (s) (Table 2). D is the diffusivity coefficient and corresponds to $D_{\mathrm{SC}}^{\mathrm{eff}}$ for model B and $D_{\mathrm{SC}}^{\mathrm{D}}$ for model D. For model C, D corresponds to D^C(α) using Eq. (89). The air gap for model C can also be directly calculated using only the air gap values of models B and D and a fitting function similar to Eq. (89),

$$\begin{array}{}\text{(92)}& {\displaystyle \frac{h_{\text{air gap}}^{\mathrm{C}}-h_{\text{air gap}}^{\mathrm{B}}}{h_{\text{air gap}}^{\mathrm{D}}-h_{\text{air gap}}^{\mathrm{B}}}}={\displaystyle \frac{A\mathit{\alpha }}{\mathrm{1}+A\mathit{\alpha }}},\end{array}$$

with A=1200.

Figure 17 (a, b, c) Density profiles from the macroscopic models and from the experimental data for the final stage of the Bouvet A, Bouvet B, and Kamata experiments. Results of model A are provided for α values from 10⁻⁵ to 10⁻⁹ and for the parameter set Calonne. The height of the air gaps derived for models B and D are shown with horizontal bars in the insets with an enlarged portion of the graph. Results of model B are provided for the SC set (dashed orange lines) and the Calonne set (solid orange lines). Results of model D are provided for the SC set (dashed blue lines) and the Fit set (solid blue lines). (d, e, f) Air gap height at the final stage of the experiments as a function of α, simulated with models A, B, C, and D using the parameter sets of SC, Calonne, and Fit. For model C, the air gap calculated with Eq. (92) is also shown with black lines.

Figure 17 shows the vertical profile of the density and the height of the air gap simulated and measured in Bouvet A, Bouvet B, and Kamata at the end of the experiments. The Bouvet B and Kamata experiments report a mass loss in the lower part of the snow layer. In Bouvet B, it results in the formation of an air gap of 2.7 mm height at the layer base or the point at which the snow density drops from about 290 to 0 kg m⁻³ within a few millimeters (Fig. 17b). In Kamata, the initial uniform density profile around 165 kg m⁻³ evolved and at the final stage showed a density of 152 kg m⁻³ at the bottom of the layer which is lower than elsewhere, where density is around 170 kg m⁻³ (Fig. 17c). So only a decrease in density at the base was observed and not an air gap. However, this might have been prevented by the vertical resolution of the density measurement of 2.5 cm in Kamata, at which point the detection of a millimeter-scale air gap is not possible. To provide an estimation of the height of the potential air gap, we converted the density decrease in the first 2.5 cm at the bottom into a air gap. This would lead to a 2.6 mm height air gap, similar to the one measured for Bouvet B for a much lower temperature gradient. Finally, as already mentioned, the experiment of Bouvet A does not include the first millimeters at the base of the snow layer, so a comparison with the simulations is not possible.

Looking at the experiments Bouvet A and Bouvet B, a first description of model A is given, with α ranging from 10⁻⁹ to 10⁻⁵. The model predicts similar mass transport for both experiments; there is a mass gain in the upper part of the layer and a mass loss in the lower part, with the latter feature being consistent with the measurements. In more detail, and as for temperature, the impact of α is clearly shown. For the lowest α value of 10⁻⁹, the density profile is almost linear, so the mass redistribution is even throughout the layer. As α increases, the area of mass loss and mass gain becomes more localized near the base and top of the layer, and the density transitions become sharper. From $\mathit{\alpha }={\mathrm{10}}^{-\mathrm{6}}$ and above, an air gap is simulated with density values reaching 0 kg m⁻³ at the bottom. The air gap closest to the experiments is obtained with the highest $\mathit{\alpha }={\mathrm{10}}^{-\mathrm{5}}$ and reaches 2.3 mm height for Bouvet A and 1.65 mm for Bouvet B, which corresponds fairly well to the measured air gap, yet is slightly underestimated (75 %) (Fig. 17d and e). Approximations of the air gap for models B and D are close to the ones simulated by model A so that all of the models seem to underestimate the air gap. For Bouvet B, an air gap of 1.7 mm is estimated for model D (63 % of the experimental air gap) and of 1.1 mm for model B (41 % of the experimental air gap). Finally, for all of the models, using the alternative SC set of input parameters (Table 3) has little impact on the air gap and leads mostly to a slight reduction in height (dashed lines in Fig. 17).

Simulations of the Kamata experiment with model A differ from the ones of Bouvet A and Bouvet B. Indeed, a mass gain is not predicted in the upper part of the snow layer but instead in a zone right above the mass loss region. This is particularly visible for $\mathit{\alpha }={\mathrm{10}}^{-\mathrm{5}}$, where the air gap, located in the first 2.1 mm, is directly surmounted by the densest part of the snow layer, located around 4 mm, with a density reaching 175 kg m⁻³. Simulations seem to show that the mass redistribution in the Kamata experiment occurs mostly in the lower part of the snow layer, which might be due to the impact of temperature on the simulated heat and mass transport processes so that they are reduced in the very cold upper part (−65 °C). This effect of temperature can also be seen in the $\dot{\mathit{\varphi }}$ simulations on the simplified microstructure for the temperature gradient of 500 K m⁻¹ for which the imposed temperature conditions were close to the ones of Kamata, as presented in Fig. 12f. Considering that this effect applies in reality, it would imply that the final density measured in Kamata in the first 2.5 cm at the bottom is the result of both the mass loss and mass gain and thus that the above estimation of an air gap of 2.6 mm might be an underestimation. Comparisons of air gaps for the Kamata experiment should be looked at with the above consideration in mind. When averaging the simulated density values of model A over a 2.5 cm step, as done in the measurements, a value of 158 kg m⁻³ is found for the first 2.5 cm, in agreement with the measured one of 152 kg m⁻³. Coming back to the air gap comparisons, model A predicts the air gap estimated for the Kamata experiment well when the highest α value is considered. At $\mathit{\alpha }={\mathrm{10}}^{-\mathrm{5}}$, the simulated air gap is 2 mm; again, this close to the estimated one of 2.6 mm despite being underestimated (77 %). Unlike for Bouvet A and Bouvet B, models B and D stand out from model A, and their approximations of the air gap are significantly larger at between 3 mm and 4 mm. They would thus predict a larger air gap than the one observed in the experiment, which is up to twice the height.

5 Discussion

5.1 Modeling heat and mass transfer with models A, B, C, or D

In the present work, macroscopic models for heat and mass transfer in dry snow have been derived by homogenization from the physics at the pore scale for different values of the condensation coefficient α in the range [10⁻¹⁰, 1]. The latter was assumed to be constant in the whole modeled snow layer. The Robin boundary equation for the water vapor at the ice–air interface allowed us to define a transition value α_T, which equals $\approx \mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ for a typical snow grain size of around 0.5 mm that characterizes the transition between the two main mechanisms driving the water vapor transfer through the snowpack, namely diffusion and sublimation–deposition. The homogenization process allowed us (i) to retrieve three different models already proposed in the literature (Calonne et al., 2014b; Hansen and Foslien, 2015; Moyne et al., 1988); to specify their domains of validity according to the α values; and (ii) to show that the hypothesis ρ_v=ρ_vs(T), which is often made, is a good approximation for α values larger than 10⁻⁶.

At the macroscopic scale, model A (Calonne et al., 2014b), valid for α values in the range [10⁻¹⁰, 10⁻⁵], is described by two coupled equations, with one for the temperature field and one for the water vapor field. They are coupled by a source term that reflects the sublimation–deposition process and depends on α. In this model, the induced porosity variation in the snow layer can be easily computed. In the case of models B, C, and D (Moyne et al., 1988; Hansen and Foslien, 2015), the physics at the macro-scale is driven by the temperature field only as ρ_v=ρ_vs(T). Because the models only solve the temperature field, it is not as straightforward to access the porosity variation. In our case, both models do not satisfy mass conservation and predict only the deposition over the whole snow layer and also the sublimation front occurring at the bottom of the snow layer, as seen in the comparison between pore-scale and macroscopic-scale simulations (Fig. 12). In the future, a more reliable boundary condition, such as a Stefan boundary condition, should be introduced to better describe the evolution of the sublimation front.

According to their definition, models B and D do not depend on α and are able to describe the observed plateaus in a limited range of α. These two models can be seen as particular cases of models A and C. Consequently, in practice, models A and C, which depend on α, are sufficient to describe the macroscopic heat and mass transfer through the snowpack for α values in the range 10⁻¹⁰ to 1. Let us remark that model C requires solving a fully coupled problem at the REV scale in order to determine the effective parameters.

5.2 On the comparisons between simulations and measurements

5.2.1 Summary of the models' evaluation

When comparing experiments and simulations with the three models, it appears that they are able to reproduce the main features of the heat and mass transport during the TGM, including the non-linear temperature profile and, for model A, the upward vapor transport with, eventually, the formation of a millimeter-scale basal air gap. However, a major discrepancy lies in the fact that temperature values are underestimated by all the models. More precisely, the heat source inducing the non-linearity in the temperature profile seems underestimated. The best predictions of the temperature deviation ΔT are obtained by model D and correspond to only about 25 % of the experimental data, translating into temperature differences of around 1 and 5 K for the Bouvet A and Kamata experiments, respectively. To a much lesser extent, upward vapor transport seems slightly underestimated, and the heights of the basal air gaps simulated by model A correspond to about 75 % of the experimental ones, leading to small differences in height of 1 and 0.6 mm for the Bouvet B and Kamata experiments, respectively. Similar conclusions can be drawn for models B and D, based on rough approximations of the air gap. Possible causes of the differences between experiments and simulations are explored in the following.

5.2.2 Uncertainties in the experimental data

Temperature measurements in Bouvet A were performed with Pt100 sensors with an accuracy of ±0.2 °C (Bouvet et al., 2023). Copper Constantan thermocouples were used in the experiment of Kamata and Sato (2007) and are known to be very stable at low temperatures, with an accuracy of ±0.5 °C. In both cases, these uncertainties are smaller than the discrepancies between the measured and modeled ΔT. For density, the experimental setup of Bouvet B ensures precise monitoring of the mass change over time by tomography (Bouvet et al., 2023). Air gaps similar to the one in Bouvet B were reported by Wiese (2017) during temperature gradient experiments. For the Kamata experiment, the reliability of the compartment method is less obvious, and the vertical resolution of 2.5 cm is rather poor to assess the presence of an air gap.

5.2.3 Uncertainties in the numerical simulation input

The macroscopic modeling of heat and mass transport in dry snow relies on the effective parameters of k^eff, D^eff, and the SSA for model A; on ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{B}}$ for model B; and on ${\stackrel{\mathrm{̃}}{k}}^{\mathrm{D}}$ for model D. For thermal conductivity, the different estimates used in the simulations are overall in good agreement with the values computed from the experimental 3D images. A possible way to improve could be to account for the anisotropy of the property so that there is, e.g., an enhanced thermal conductivity in the vertical direction, as observed for snow evolving under a high TGM (e.g., Calonne et al., 2011). However, an increase in the thermal conductivity of models A, B, and D leads to a decrease in both the temperature deviation and air gap height and thus a degradation in the model performance, as illustrated in Fig. 16. Concerning the vapor diffusion coefficient, the SC models were used in the simulations, which provide slightly overestimated estimates. However, improving these estimates, through taking lower values of the diffusion coefficient, leads to a decrease in both temperature deviation and air gap height, which, again, degrades the model results. To illustrate this behavior in the case of model A, by lowering the SC estimate of D^eff by 10 %, a ΔT of 0.096 K and an air gap of 1.17 mm is simulated for $\mathit{\alpha }={\mathrm{10}}^{-\mathrm{5}}$, compared to 0.102 K and 2.3 mm with the initial value. Finally, potential errors in the SSA parameter would affect the source term of model A and would only translate in small variations in α. To conclude, the uncertainties linked to the estimation of the effective parameters cannot be responsible for the reported differences between experiments and simulations.

5.2.4 Models limitations and potential improvements

As the points raised above do not seem sufficient to explain the model errors, a plausible cause remains to be investigated and is the definition of the model itself, i.e., the definition of the physics at the pore scale considered for the homogenization. A first element concerns the source terms in model A, which are derived from the Hertz–Knudsen equation and rely on a condensation coefficient α (Eq. 7). Here, this coefficient was taken to be constant and uniform over the snow layer and considered equal for both condensation and sublimation. In their review, Persad and Ward (2016) explore the expressions of the evaporation coefficient and of the condensation coefficient in the Hertz–Knudsen equation for the water–air interface. They conclude that most errors come from assuming the evaporation and condensation coefficients to be equal and assuming thermal equilibrium across the liquid–vapor interface (Eq. 4 in this study). Moreover, as mentioned in the Introduction, the condensation parameter α depends on many parameters, such as the vapor supersaturation, which can lead to a non-linear expression of the Hertz–Knudsen equation. Hence, refining the Hertz–Knudsen equation could add non-linearity in the source terms of model A, which could enhance the contribution of latent heat and thus increase the temperature deviation ΔT, which would improve the model's prediction.

Another point is that the natural convection was not taken into account at the pore scale. This process was, however, hypothesized to be key for heat and mass transport of snow under strong temperature gradients, such as Arctic and sub-Arctic ones (e.g., Sturm and Johnson, 1991; Domine et al., 2018). To include natural convection, fluxes of temperature (J_T) and water vapor ($J_{{\mathit{\rho }}_{\mathrm{v}}}$) should be expressed at the pore scale as follows:

$$\begin{array}{}\text{(93)}& \begin{array}{rl}{J}_{\mathrm{T}}& =-k_{\mathrm{a}}\mathbf{grad}{T}_{\mathrm{a}}+{\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{a}}{v}_{\mathrm{a}}{T}_{\mathrm{a}}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\\ & \mathrm{and}{J}_{{\mathit{\rho }}_{\mathrm{v}}}=-D_{\mathrm{v}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{v}}+v_{\mathrm{a}}{\mathit{\rho }}_{\mathrm{v}}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\end{array}\end{array}$$

with v_a as the air velocity. A numerical study was recently presented by Jafari et al. (2022) using a macroscopic model similar to model A. They show that the occurrence and intensity of natural convection in snow depends on the Rayleigh number, which is defined as

$$\begin{array}{}\text{(94)}& \mathit{Ra}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{\mathit{\beta }}_{\mathrm{T}}g(T_{\mathrm{bottom}}-T_{\mathrm{top}})HK}{\left(\right({\mathit{\mu }}_{\mathrm{a}}{k}^{\mathrm{eff}})/({\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{a}}\left)\right)}},\end{array}$$

where H is the height of the snow layer, K is the snow permeability, g=9.81 m s⁻² is the gravity, ${\mathit{\mu }}_{\mathrm{a}}=\mathrm{17.29}\times {\mathrm{10}}^{-\mathrm{6}}$ Pa s is the air viscosity, and β_T=0.0036 K⁻¹ is the thermal expansion coefficient. Their simulations indicate that, for Ra>50 and H>25 cm, the natural convection could generate an upward air flux from the warmer region to the colder one, and vice versa. We estimated the Rayleigh number for the three experiments used in this study. Using the values in Tables 1 and 2, and using the parameterization of Calonne et al. (2012) for the snow permeability, the Rayleigh number is typically 0.15, 0.02, and 0.85 for the experiments Bouvet A, Bouvet B and Kamata, respectively. These values are much smaller than the threshold value presented by Jafari et al. (2022), which would indicate that natural convection is negligible in our cases. Moreover, Kamata et al. (1999) present a symmetric TGM experiment with warmer conditions at the base and top of the snow layer and colder conditions imposed in the middle using a cold plate. The snow layer was thus under a positive TG in one part and under a negative TG in the other part, with both having the same intensity. The temperature profiles were recorded in both parts of the snow layer and show similar non-linear curves in both cases, although natural convection could only occur in the bottom area, where the temperature conditions are unstable. The authors conclude that natural convection does not seem to impact their temperature fields. This would be consistent with the small Rayleigh number that we estimated to be 0.16 for this experiment.

Finally, the cross-coupling effects between the temperature and water vapor density, such as the Soret and Dufour effects, were not considered in the physics at the pore scale. The effect of the vapor density gradient on the heat flux, called the Dufour effect, is characterized by the thermodiffusion coefficient D^Tv, and the effect of temperature gradient on the vapor density flux, called the Soret effect, is characterized by the thermodiffusion coefficient D^vT. Taking these effects into account, the temperature and water vapor flux can be expressed as follows:

$$\begin{array}{}\text{(95)}& \begin{array}{rl}{J}_{\mathrm{T}}& =-k_{\mathrm{a}}\mathbf{grad}{T}_{\mathrm{a}}-D^{\mathrm{Tv}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{v}}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}},\\ & \mathrm{and}{J}_{{\mathit{\rho }}_{\mathrm{v}}}=-D_{\mathrm{v}}\mathbf{grad}{\mathit{\rho }}_{\mathrm{v}}-D^{\mathrm{vT}}\mathbf{grad}{T}_{\mathrm{a}}\mathrm{in}{\mathrm{\Omega }}_{\mathrm{a}}.\end{array}\end{array}$$

For porous media, the Dufour effect is neglected in most cases, whereas the Soret effect is often taken into account (e.g., Davarzani et al., 2010; Häussling Löwgren et al., 2020) and can be measured using the Soret coefficient defined as $S_{\mathrm{T}}=D^{\mathrm{vT}}/D_{\mathrm{v}}$. This coefficient is positive when the heaviest species in the pore space move toward the colder regions and is negative when they move toward the warmer regions. However, this coefficient could change sign when the temperature is lowered (Caldwell, 1973; Chapman and Cowling, 1990). When the temperature is positive, the Soret coefficient for water vapor is supposed to be positive. To the best of our knowledge, there are no data concerning this coefficient when temperature is negative. The Soret effect can be easily introduced in pore-scale simulations in the case of the 2D simplified microstructure as presented in Sect. 3.3. By doing so, we found that negative S_T coefficients lead to an increase the simulated temperature deviation ΔT, and vice versa. For example, for α=0.1 and a temperature gradient of 500 K m⁻¹, a maximum value of ΔT of 6.5, 4.2, and 2.8 K is simulated for S_T equal to −2 × 10⁻⁴, 0, and 2 × 10⁻⁴, respectively. The Soret effect can also be introduced in model D, by replacing k_a+k_diff with $k_{\mathrm{a}}+k_{\mathrm{diff}}+S_{\mathrm{T}}{D}_{\mathrm{v}}{L}_{\mathrm{sg}}/{\mathit{\rho }}_{\mathrm{i}}$. Using the self-consistent estimate of thermal conductivity (Eq. 87) and for values of S_T at −2 × 10⁻⁴, 0, and 2 × 10⁻⁴, the maximum simulated values of ΔT for the Kamata experiment are 6.3, 3.1, and 1.95 K, respectively, whereas the experimental value is around 6 K. A negative Soret coefficient seems thus suitable to improve the temperature simulations and better describe the experimental data. However, the influence of the Soret effect on the air gap is not straightforward, as it seems to induce a downward movement of vapor molecules and is thus opposed to the formation of a basal air gap. These preliminary results show that the introduction of such coupling effects (Soret and/or Dufour) between the temperature and the water vapor density in the modeling of heat and water vapor transfer in snowpacks is interesting and that future work would be needed to investigate such a hypothesis.

6 Conclusion

This paper presents the definition and evaluation of the equivalent macroscopic modeling of heat and mass transport during the TGM in dry snow. First, we applied the homogenization process to retrieve the macroscopic models valid for condensation coefficients α ranging from 10⁻¹⁰ to 1. We showed that, at a transition value ${\mathit{\alpha }}_{\mathrm{T}}\approx \mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$, the modeling changes from vapor transport limited by sublimation–deposition (models A and B) to vapor transport limited by diffusion (model D). The homogenization process allowed us to retrieve different models proposed in the literature (Calonne et al., 2014b; Hansen and Foslien, 2015; Moyne et al., 1988) and to clarify their domains of validity according to the value of α. Models A and C are sufficient to describe the heat and mass transfer in the whole range of α values between 10⁻¹⁰ and 1. For α values between 10⁻¹⁰ and 10⁻⁵, model A consists of two equations of temperature and water vapor density coupled through the source terms, which are proportional to the Hertz–Knudsen equation and therefore to α. This model does not presume any assumption on the saturation of the vapor density. For α values between 10⁻⁵ and 1, model C consists of one temperature equation which involves α, since the hypothesis ρ_v=ρ_vs(T), which is often made, is valid in that range. Concerning the two other models, model B can be seen as a particular case of model A for α values in the range from 10⁻⁷ to 10⁻⁵, whereas model D can be seen as a particular case of model C for α values in the range from 10⁻³ to 1.

In the second part of the paper, we evaluated the homogenized models A, B, C, and D by comparison with three laboratory experiments of the TGM of snow (Kamata and Sato, 2007; Bouvet et al., 2023), as well as by a numerical evaluation for a 2D simplified microstructure. Evaluations were performed based on the temperature and density profiles of snow and, more precisely, on the ability to reproduce two main features reported in the TGM experiments, namely the non-linear concave-shaped temperature profile, characterized by the temperature deviation from a linear gradient ΔT, and the upward vapor transport leading to a mass loss or an air gap at the base of the snow layer. We showed that (i) the four models allow the reproduction of the shape of the temperature profile, but the values were largely underestimated, with the best prediction being obtained with model D and corresponding only to 25 % of the temperature difference observed in the experimental data. This major discrepancy highlights that a process that contributes to heat up the layer is not well captured, if at all; (ii) model A allows us to reproduce the upward vapor transport and the formation of a millimeter-scale basal air gap, with the best result being obtained for the highest α value of 10⁻⁵; and (iii) models B and D (and also C) do not allow the reproduction of mass transport as they predict only the mass gain in the snow layer and as they do not satisfy mass conservation in the present case. Potential improvements were suggested and include the refining or enrichment of the physics at the pore scale considered to derive the models, such as questioning the expression of the Hertz–Knudsen equation or the role of the Soret and/or Dufour effects, as well as improving the boundary conditions to allow for realistic mass transport for models B, C, and D.