# Multiscale Modeling Of Heat And Mass Transfer In Dry Snow: Influence Of The Condensation Coefficient And Comparison With Experiments

Lisa Bouvet, Neige Calonne, Frédéric Flin, Christian Geindreau

► To cite this version:

Lisa Bouvet, Neige Calonne, Frédéric Flin, Christian Geindreau. Multiscale Modeling Of Heat And Mass Transfer In Dry Snow: Influence Of The Condensation Coefficient And Comparison With Experiments. The Cryosphere, 2024, 18, pp.4285-4313. 10.5194/tc-18-4285-2024. insu-04752890

HAL Id: insu-04752890

https://insu.hal.science/insu-04752890v1

Submitted on 25 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multiscale modeling of heat and mass transfer in dry snow: influence of the condensation coefficient and comparison with experiments

Lisa Bouvet<fnref n="1" />,<fnref n="2" />, Neige Calonne<fnref n="1" />, Frédéric Flin<fnref n="1" />, and Christian Geindreau<fnref n="2" />

**Correspondence:** Lisa Bouvet (lisa.bouvet@meteo.fr)

Received: 22 September 2023 – Discussion started: 9 November 2023

Revised: 4 June 2024 – Accepted: 2 July 2024 – Published: 19 September 2024

**Abstract.** 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 con- densation 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 mi- crostructures (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 experi- mental 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. How- ever, both features are underestimated overall by the mod- els 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 phe- nomena 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 meta- morphism (TGM), which transforms snow into faceted crys- tals (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 in- duces significant changes in the snow physical and mechan- ical properties, such as thermal conductivity, vapor diffusiv- ity, or elastic properties (e.g., Srivastava et al., 2010; Calonne et al., 2014a; Wautier et al., 2015), affecting the snowpack

<footnote n="1">¹Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Études de la Neige, Grenoble, France</footnote>

<footnote n="2">²Université Grenoble Alpes, CNRS, Grenoble INP, 3SR, Grenoble, France</footnote>
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 $\alpha$, 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^{-4}$ to $10^{-1}$ (see, e.g., Libbrecht and Rickerby, 2013). Indirect estimates from snow modeling at the pore scale range from $10^{-4}$ to $10^{-3}$ (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^{\text{eff}}$, which depends on the ice and air conductivity and on the snow microstructure, and (ii) the effective diffusion $D^{\text{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 $\alpha$. Calonne et al. (2014b) have shown that this model is valid for interface growth velocities below $3 \times 10^{-11} \text{ m s}^{-1}$, 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 $\alpha$ (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 $\alpha = 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 $\alpha$ 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^{\text{eff}}$ for $\alpha$ values smaller than $\approx 10^{-4}$ and then increases with increasing $\alpha$ until it reaches a plateau for $\alpha$ values larger than $10^{-2}$, 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 $\alpha$ values. This should be clarified, especially by estimating the $\alpha$ values for which the assumption of saturated water vapor is theoret-
ically 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^{-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 $\varepsilon = l/L \ll 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 $\varepsilon = l/L \ll 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 $\varepsilon$; (iv) looking for the unknown fields in the form

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Physical phenomena under consideration at the representative elementary volume (REV) scale.</figcaption></figure>

of asymptotic expansions in powers of $\varepsilon$; 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 $\varepsilon = l/L \ll 1$.

In what follows, $\Omega$ is the REV domain, $\Omega_i$ is the ice domain, and $\Omega_a$ is the air domain (Fig. 1). The ice grain interface is denoted $\Gamma$, and $n_i$ is the unit's outward vector of $\Omega_i$. The subscripts i or a are related to quantities defined in $\Omega_i$ and $\Omega_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:

$$ \rho_i C_i \frac{\partial T_i}{\partial t} - \text{div}(k_i \mathbf{grad} T_i) = 0 \quad \text{in } \Omega_i, \qquad (1) $$

$$ \rho_a C_a \frac{\partial T_a}{\partial t} - \text{div}(k_a \mathbf{grad} T_a) = 0 \quad \text{in } \Omega_a, \qquad (2) $$

$$ \frac{\partial \rho_v}{\partial t} - \text{div}(D_v \mathbf{grad} \rho_v) = 0 \quad \text{in } \Omega_a, \qquad (3) $$

$$ T_i = T_a \quad \text{on } \Gamma, \qquad (4) $$

$$ k_i \mathbf{grad} T_i \cdot \mathbf{n}_i - k_a \mathbf{grad} T_a \cdot \mathbf{n}_i = L_{sg} \mathbf{w} \cdot \mathbf{n}_i \quad \text{on } \Gamma, \qquad (5) $$

$$ D_v \mathbf{grad} \rho_v \cdot \mathbf{n}_i = (\rho_i - \rho_v) \mathbf{w} \cdot \mathbf{n}_i \approx \rho_i \mathbf{w} \cdot \mathbf{n}_i \quad \text{on } \Gamma, \qquad (6) $$

where $t$ is the time (s), $T$ is the temperature (K), $k$ is the thermal conductivity (W m⁻¹ K⁻¹), $\rho$ 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⁻³), $\mathbf{w}$ is the interface growth velocity (m s⁻¹), $\rho_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 $\mathbf{X}$, respectively. At the interface, the heat and mass transfer are coupled through the normal interface growth velocity, $w_n = \mathbf{w} \cdot \mathbf{n}_i$, which is given by the Hertz–Knudsen equation,

$$ w_n = \mathbf{w} \cdot \mathbf{n}_i = \frac{1}{\beta} \left[ \frac{\rho_v - \rho_{vs}(T_a)}{\rho_{vs}(T_a)} - d_0 K \right] \quad \text{on } \Gamma, \qquad (7) $$

such that $w_n$ is positive when the ice grain grows and negative when it sublimates. $\beta$ is the interface of the kinetic coefficient (s m⁻¹), $\rho_{vs}$ is the saturation of the water vapor density in air (kg m⁻³), $d_0$ is the capillary length (m), and $K$ is the interface of the mean curvature (m⁻¹). The interface kinetic coefficient $\beta$ is linked to the condensation coefficient $\alpha$ by

$$ \frac{1}{\beta} = \alpha \frac{\rho_{vs}(T_a)}{\rho_i} \sqrt{\frac{k_B T_a}{2\pi m}}, \qquad (8) $$

where $m$ is the mass of a water molecule (kg), and $k_B$ is the Boltzmann constant equal to $1.38 \times 10^{-23} \text{ J K}^{-1}$. As already mentioned, the condensation coefficient $\alpha$ 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 $\rho_{vs}$ at a given air temperature $T_a$ is given by the Clausius–Clapeyron law as follows:

$$ \rho_{vs}(T_a) = \rho_{vs}^{\text{ref}}(T^{\text{ref}}) \exp \left[ \frac{L_{sg} m}{\rho_i k_B} \left( \frac{1}{T^{\text{ref}}} - \frac{1}{T_a} \right) \right]. \qquad (9) $$

In the current work, we chose the reference values $T^{\text{ref}} = 263 \text{ K}$, leading to a $\rho_{vs}^{\text{ref}}(T^{\text{ref}})$ value of $2.173 \times 10^{-3} \text{ kg m}^{-3}$. For simplicity, we assume that none of the material properties ($\rho$, $C$, $k_B$, $D_v$, $\beta$, 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{aligned} w_n &= \mathbf{w} \cdot \mathbf{n}_i = \frac{1}{\beta \rho_{vs}(T_a)} [\rho_v - \rho_{vs}(T_a)] \\ &= \frac{\alpha}{\rho_i} w_k(T_a) [\rho_v - \rho_{vs}(T_a)] \quad \text{on } \Gamma, \end{aligned} \qquad (10) $$

where $w_k = \sqrt{k_B T_a / 2\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{aligned} k_i \mathbf{grad} T_i \cdot \mathbf{n}_i - k_a \mathbf{grad} T_a \cdot \mathbf{n}_i &= L_{sg} \frac{\alpha}{\rho_i} w_k(T_a) [\rho_v - \rho_{vs}(T_a)] \\ &= L_{sg} \frac{D_v}{\rho_i} \mathbf{grad} \rho_v \cdot \mathbf{n}_i \quad \text{on } \Gamma, \end{aligned} \qquad (11) $$

$$ D_v \mathbf{grad} \rho_v \cdot \mathbf{n}_i = \alpha w_k(T_a) [\rho_v - \rho_{vs}(T_a)] \quad \text{on } \Gamma. \qquad (12) $$

## 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 $\varphi$ reads $\varphi = \varphi_c \varphi^*$, 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

$$ [F_i^T] \rho_i^* C_i^* \frac{\partial T_i^*}{\partial t^*} - \text{div}^*(k_i^* \mathbf{grad}^* T_i^*) = 0 \quad \text{in } \Omega_i, \qquad (13) $$

$$ [F_a^T] \rho_a^* C_a^* \frac{\partial T_a^*}{\partial t^*} - \text{div}^*(k_a^* \mathbf{grad}^* T_a^*) = 0 \quad \text{in } \Omega_a, \qquad (14) $$

$$ [F_a^\rho] \frac{\partial \rho_v^*}{\partial t^*} - \text{div}^*(D_v^* \mathbf{grad}^* \rho_v^*) = 0 \quad \text{in } \Omega_a, \qquad (15) $$

$$ T_i^* = T_a^* \quad \text{on } \Gamma, \qquad (16) $$

$$ [K] k_i^* \mathbf{grad}^* T_i^* \cdot \mathbf{n}_i - k_a^* \mathbf{grad}^* T_a^* \cdot \mathbf{n}_i = [H] L_{sg}^* \frac{D_v^*}{\rho_i^*} \mathbf{grad}^* \rho_v^* \cdot \mathbf{n}_i \quad \text{on } \Gamma, \qquad (17) $$

$$ D_v^* \mathbf{grad}^* \rho_v^* \cdot \mathbf{n}_i = [W_R] \alpha^* w_k^* [\rho_v^* - [R] \rho_{vs}^*(T_a^*)] \quad \text{on } \Gamma. \qquad (18) $$

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{aligned}
[F_i^T] &= \frac{l^2 \rho_{ic} C_{ic}}{t_c k_{ic}}, & [F_a^T] &= \frac{l^2 \rho_{ac} C_{ac}}{t_c k_{ac}}, \\
[F_a^\rho] &= \frac{l^2}{D_{vc} t_c}, & [K] &= \frac{k_{ic}}{k_{ac}}, \\
[W_R] &= \frac{l \alpha_c w_{kc}}{D_{vc}}, & [R] &= \frac{\rho_{vs_c}(T_{ac})}{\rho_{vc}},
\end{aligned}
$$

$$ 
\begin{aligned}
[H] &= \frac{l L_{sg_c} w_{nc}}{k_{ac} T_{ac}}, & \text{with } w_{nc} &= \frac{\alpha_c w_{kc}}{\rho_{ic}} (\rho_{vc} - \rho_{vs_c}(T_{ac})) \\
    &= \frac{D_{vc} \rho_{vc}}{l \rho_{ic}}.
\end{aligned}
\quad (19) $$

Dimensionless numbers $[F_i^T]$ and $[F_a^T]$ correspond to the inverse of the Fourier number in $\Omega_i$ and $\Omega_a$, respectively. They characterize the ratio between the thermal energy storage rate and the heat conduction rate. $[F_a^\rho]$ is an analogous inverse Fourier number for the transient water vapor transfer by diffusion in $\Omega_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 $\alpha$ 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_v^* \text{grad}^* \rho_v^* \cdot n_i = 0)$, whereas when $[W_R]$ tends towards infinite (or is very large), Eq. (18) is equivalent to a Dirichlet boundary condition $(\rho_v^* = \rho_{vs}(T_a^*)$).

## 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 $\varepsilon = 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 5 \times 10^{-4}$ m and $L \approx 0.1$ m, leading to $\varepsilon = 5 \times 10^{-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_{ic} = k_{ic}/(C_{ic} \rho_{ic})$ and in the air

phase $D_{ac} = k_{ac}/(C_{ac} \rho_{ac})$ are of the same order of magnitude as the vapor diffusion coefficient $D_{vc}$. Thus, the characteristic time $t_c$ associated with these transfers through the snowpack are of the same order of magnitude, where $t_c = O(L^2/D_{ic}) = O(L^2/D_{ac}) = O(L^2/D_{vc})$. Hence, from Eq. (19), we get $[F_i^T] = O([F_a^T]) = O([F_a^\rho]) = O(\varepsilon^2)$. At the ice pore interface, from Eq. (19), we have $[K] = O(1)$ and $[R] = O(1)$. The latter estimation implies that the supersaturation $\sigma = (\rho_v - \rho_{vs})/\rho_{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

$$ [W_R] = \frac{l \alpha_c w_{kc}}{D_{vc}} = \frac{l^2}{D_{vc}} \frac{\alpha_c w_{kc}}{l} = \frac{\tau_d}{\tau_{\text{sub/dep}}}, $$

where $\tau_d = l^2/D_{vc}$ is the characteristic time associated with water vapor diffusion at the pore scale, and $\tau_{\text{sub/dep}} = l/(\alpha_c w_{kc})$ 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 $\alpha_c$. Using the characteristic values given in Table 1, this ratio is equal to 1 for a particular value of $\alpha_c$, which is denoted $\alpha_T = D_{vc}/(l w_{kc}) \approx 3 \times 10^{-4}$. This value decreases when the characteristic length $l$ increases, such that the values range between $10^{-3}$ for small grains (~0.1 mm) and $10^{-5}$ 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 $\alpha_T$ value characterizes the transition between two mechanisms which drive the water vapor transfer at the pore scale. When $\tau_d \ll \tau_{\text{sub/dep}}$, i.e., for $\alpha_c \ll \alpha_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 $\tau_d \gg \tau_{\text{sub/dep}}$, i.e., for $\alpha_c \gg \alpha_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.

Estimations of the dimensionless number $[H]$ are not as straightforward, as they depend on the intensity of the interface normal growth velocity $w_{nc}$. When $\alpha_c$ is small (typically smaller than $\alpha_T$), $[W_R]$ is also small, and Eq. (18) implies that $\Delta\rho_{vc}$ has a finite value ($O(\rho_{vc})$). Thus, this dimensionless number $[H]$ can be also written as

$$ [H] = \frac{l L_{sg_c} \alpha_c w_{kc} \rho_{vc}}{k_{ac} T_{ac} \rho_{ic}}. $$

In that case, it increases when $\alpha_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 $\alpha_c$ (typically larger than $\alpha_T$), Eq. (18) implies that $\rho_v^* \approx \rho_{vs}(T_a^*)$. As a consequence, from Eqs. (17) and (19), $[H]$ can be rewritten as

$$ [H] = \frac{l L_{sg_c} D_{vc} \gamma(T_{ac}) T_{ac}}{l \rho_{ic} k_{ac} T_{ac}} = \frac{L_{sg_c} D_{vc} \gamma(T_{ac})}{\rho_{ic} k_{ac}} = \frac{k_{\text{difc}}}{k_{ac}}, $$

where $\gamma(T_{ac}) = d\rho_{vs}(T_{ac})/dT_{ac}$ is the derivative of the Clausius–Capeyron law, and $k_{\text{difc}} = L_{sg_c} D_{vc} \gamma(T_{ac})/\rho_{ic}$ can
<table><caption>**Table 1.** Characteristic values of the properties evaluated at -10 °C from the literature (Massman, 1998; Kaempfer and Plapp, 2009).</caption><thead><tr><th>Symbol</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>T<sub>ic</sub>, T<sub>ac</sub></td><td>Temperature of ice, air</td><td>263 K</td></tr><tr><td>k<sub>ic</sub></td><td>Heat conductivity of ice</td><td>2.3 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>k<sub>ac</sub></td><td>Heat conductivity of air</td><td>0.024 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>C<sub>ic</sub></td><td>Specific heat capacity of ice</td><td>2000 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>C<sub>ac</sub></td><td>Specific heat capacity of air</td><td>1005 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>L<sub>sgc</sub></td><td>Latent heat of sublimation of ice</td><td>2.60 &times; 10<sup>9</sup> J m<sup>-3</sup></td></tr><tr><td>D<sub>vc</sub></td><td>Water vapor diffusion coefficient in air</td><td>2.036 &times; 10<sup>-5</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>&rho;<sub>vc</sub></td><td>Water vapor density in air</td><td>0.002 kg m<sup>-3</sup></td></tr><tr><td>&rho;<sub>ic</sub></td><td>Ice density</td><td>917 kg m<sup>-3</sup></td></tr><tr><td>&rho;<sub>ac</sub></td><td>Air density</td><td>1.335 kg m<sup>-3</sup></td></tr><tr><td>l</td><td>Microscopic length</td><td>5 &times; 10<sup>-4</sup> m</td></tr><tr><td>L</td><td>Macroscopic length</td><td>0.1 m</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Estimation of the dimensionless number [$W_R$] with respect to $\alpha$, which leads to several cases of macroscopic modeling to be considered (cases A to D2). The $\alpha_T$ value characterizes the transition between two cases presenting different limiting processes for the water vapor transfer at the pore scale, so that $\tau_d < \tau_{sub/dep}$ or $\tau_d > \tau_{sub/dep}$, with $\tau_d$ as the characteristic time associated with water vapor diffusion and $\tau_{sub/dep}$ as the characteristic time associated with the sublimation-deposition process. $\alpha_T$ was estimated based on the characteristic values given in Table 1.</figcaption></figure>

be seen as an enhancement of the air thermal conductivity. Using the characteristic values given in Table 1 and the Clausius–Capeyron equation, Eq. (9), for large values of $\alpha_c$, [$H$] = $O(1)$. According to the above analysis, several cases must be considered, depending on the value of the condensation coefficient $\alpha_c$ (Fig. 2).

- Case A is $\tau_d = O(\varepsilon^2 \tau_{sub/dep})$, i.e., [$W_R$] = $O(\varepsilon^2)$ and [$H$] = $O(\varepsilon^2)$.

- Case B is $\tau_d = O(\varepsilon \tau_{sub/dep})$, i.e., [$W_R$] = $O(\varepsilon)$ and [$H$] = $O(\varepsilon)$.

- Case C is $\tau_d = O(\tau_{sub/dep})$, i.e., [$W_R$] = $O(1)$ and [$H$] = $O(1)$.

- Case D1 is $\tau_d = O(\varepsilon^{-1} \tau_{sub/dep})$, i.e., [$W_R$] = $O(\varepsilon^{-1})$ and [$H$] = $O(1)$.

- Case D2 is $\tau_d = O(\varepsilon^{-2} \tau_{sub/dep})$, i.e., [$W_R$] = $O(\varepsilon^{-2})$ and [$H$] = $O(1)$.

Cases A and B correspond to $0 \le \alpha \ll \alpha_T$, whereas cases D1 and D2 correspond to $\alpha_T \ll \alpha \ll 1$. Case C ensures the transition between cases B and D1, when $\alpha \approx \alpha_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, $x^* = X/L$ and $y^* = 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 $x^* = \varepsilon y^*$. When l is used as the characteristic length, the dimensionless derivative operator $\mathbf{grad}^*$ becomes $(\mathbf{grad}_{y^*} + \varepsilon \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_i^*$, the air temperature $T_a^*$, and the water vapor $\rho_v^*$ are sought in the form of the asymptotic expansion of the powers of $\varepsilon$, as follows:

$$ \varphi^*(\mathbf{x}^*, \mathbf{y}^*, t) = \varphi^{*(0)}(\mathbf{x}^*, \mathbf{y}^*, t) + \varepsilon\varphi^{*(1)}(\mathbf{x}^*, \mathbf{y}^*, t) + \varepsilon^2\varphi^{*(2)}(\mathbf{x}^*, \mathbf{y}^*, t) + \dots, \quad (20) $$

where $\varphi^* = T_i^*, T_a^*, \rho_v^*$, and the corresponding $\varphi^{*(i)}$ are periodic functions of period $\Omega$ with respect to the space variable $\mathbf{y}^*$. Substituting these expansions in the set (13)–(18) gives, through the identification of identical powers of $\varepsilon$, 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 $[H] = O(\varepsilon^2)$, $[W_R] = O(\varepsilon^2)$, 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

$$ (\rho C)^{\text{eff}} \frac{\partial T^{(0)}}{\partial t} - \text{div}(\mathbf{k}^{\text{eff}} \mathbf{grad} \, T^{(0)}) = \text{SSA}_V L_{sg} w_n^{(0)} \\ = -L_{sg} \dot{\phi}, \quad (21) $$

$$ \phi \frac{\partial \rho_v^{(0)}}{\partial t} - \text{div}(\mathbf{D}^{\text{eff}} \mathbf{grad} \, \rho_v^{(0)}) = -\text{SSA}_V \rho_i w_n^{(0)} = \rho_i \dot{\phi}, \quad (22) $$

where $w_n^{(0)}$ 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),

$$ w_n^{(0)} = \frac{\alpha}{\rho_i} w_k \left[ \rho_v^{(0)} - \rho_{vs}^{(0)} (T^{(0)}) \right], \quad (23) $$

$$ \rho_{vs}^{(0)}(T^{(0)}) = \rho_{vs}^{\text{ref}} \exp \left[ \frac{L_{sg}m}{\rho_i k} \left( \frac{1}{T^{\text{ref}}} - \frac{1}{T^{(0)}} \right) \right], \quad (24) $$

and where $\phi$ is the porosity and $\dot{\phi}$ its total time derivative. $SSA_V = |\Gamma|/|\Omega|$ 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 $SSA_V = SSA \times \rho_i$. $(\rho C)^{\text{eff}}$ is the effective thermal capacity (Eq. S.A.45), $\mathbf{k}^{\text{eff}}$ is the effective thermal conductivity tensor (Eq. S.A.46), and $\mathbf{D}^{\text{eff}}$ is the effective diffusion tensor (Eq. S.A.48). These effective properties are defined as

$$ (\rho C)^{\text{eff}} = (1 - \phi) \rho_i C_i + \phi \rho_a C_a, \quad (25) $$

$$ \mathbf{k}^{\text{eff}} = \frac{1}{|\Omega|} \left( \int_{\Omega_a} k_a (\mathbf{grad} \, \mathbf{t}_a + \mathbf{I}) d\Omega + \int_{\Omega_i} k_i (\mathbf{grad} \, \mathbf{t}_i + \mathbf{I}) d\Omega \right), \quad (26) $$

$$ \mathbf{D}^{\text{eff}} = \frac{1}{|\Omega|} \int_{\Omega_a} D_v (\mathbf{grad} \, \mathbf{g}_v + \mathbf{I}) d\Omega, \quad (27) $$

where $\mathbf{t}_a$ and $\mathbf{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:

$$ \text{div}(k_i(\mathbf{grad} \, \mathbf{t}_i + \mathbf{I})) = 0 \quad \text{in } \Omega_i, \quad (28) $$

$$ \text{div}(k_a(\mathbf{grad} \, \mathbf{t}_a + \mathbf{I})) = 0 \quad \text{in } \Omega_a, \quad (29) $$

$$ \mathbf{t}_i = \mathbf{t}_a \quad \text{on } \Gamma, \quad (30) $$

$$ (k_i(\mathbf{grad} \, \mathbf{t}_i + \mathbf{I}) - k_a(\mathbf{grad} \, \mathbf{t}_a + \mathbf{I})) \cdot \mathbf{n}_i = 0 \quad \text{on } \Gamma, \quad (31) $$

$$ \frac{1}{|\Omega|} \int_{\Omega} (\mathbf{t}_a + \mathbf{t}_i) d\Omega = \mathbf{0}. \quad (32) $$

Moreover, $\mathbf{g}_v$ is a periodic vector solution of the following boundary value problem over the REV Eqs. (S.A.35)–(S.A.37):

$$ \text{div}(D_v(\mathbf{grad} \, \mathbf{g}_v + \mathbf{I})) = 0 \quad \text{in } \Omega_a, \quad (33) $$

$$ D_v(\mathbf{grad} \, \mathbf{g}_v + \mathbf{I}) \cdot \mathbf{n}_i = 0 \quad \text{on } \Gamma, \quad (34) $$

$$ \frac{1}{|\Omega|} \int_{\Omega_a} \mathbf{g}_v d\Omega = \mathbf{0}. \quad (35) $$

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^{(0)}$ and $\rho_v^{(0)}$. These equations involve two effective parameters, namely the effective thermal conductivity $\mathbf{k}^{\text{eff}} = \mathbf{k}^{\text{eff}}(k_i, k_a, \text{microstructure})$ and the effective diffusion $\mathbf{D}^{\text{eff}} = \mathbf{D}^{\text{eff}}(D_v, \text{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

$$ (\rho C)^{\text{eff}} \frac{\partial T^{(0)}}{\partial t} - \text{div}(\mathbf{k}^{\text{eff}} \mathbf{grad} T^{(0)}) = -L_{\text{sg}} \dot{\phi}, \quad (36) $$

$$ \phi \frac{\partial \rho_{\text{vs}}^{(0)}}{\partial t} - \text{div}(\mathbf{D}^{\text{eff}} \mathbf{grad} \rho_{\text{vs}}^{(0)}) = \rho_i \dot{\phi}, \quad (37) $$

with

$$ \rho_{\text{v}}^{(0)} = \rho_{\text{vs}}^{(0)}(T^{(0)}), \quad (38) $$

where $(\rho C)^{\text{eff}}$ is the effective thermal capacity, $\mathbf{k}^{\text{eff}}$ is the effective thermal conductivity tensor, and $\mathbf{D}^{\text{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 $\rho_v^{(0)} = \rho_{vs}^{(0)}(T^{(0)})$ is directly given by the Clausius–Clapeyron law (Eq. 24). Consequently, from Eq. (37), we have the following:

$$ \begin{aligned} \dot{\phi} &= -\frac{1}{\rho_i} \left( \text{div}(\mathbf{D}^{\text{eff}} \mathbf{grad} \rho_{\text{vs}}^{(0)}(T^{(0)})) - \phi \frac{\partial \rho_{\text{vs}}^{(0)}}{\partial t} \right) \\ &= -\frac{1}{\rho_i} \left( \text{div}(\gamma(T^{(0)}) \mathbf{D}^{\text{eff}} \mathbf{grad} T^{(0)}) - \phi \gamma(T^{(0)}) \frac{\partial T^{(0)}}{\partial t} \right), \end{aligned} \quad (39) $$

where

$$ \begin{aligned} \gamma(T^{(0)}) &= \frac{d\rho_{\text{vs}}^{(0)}(T^{(0)})}{dT^{(0)}} \\ &= \rho_{\text{vs}}^{\text{ref}} \frac{L_{\text{sg}}m}{\rho_i k} \frac{1}{(T^{(0)})^2} \exp\left[\frac{L_{\text{sg}}m}{\rho_i k} \left(\frac{1}{T^{\text{ref}}} - \frac{1}{T^{(0)}}\right)\right] \\ &= \frac{L_{\text{sg}}m}{\rho_i k} \frac{1}{(T^{(0)})^2} \rho_{\text{vs}}^{(0)}(T^{(0)}). \end{aligned} \quad (40) $$

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

$$ \left( (\rho C)^{\text{eff}} + \phi \gamma(T^{(0)}) \frac{L_{\text{sg}}}{\rho_i} \right) \frac{\partial T^{(0)}}{\partial t} - \text{div}(\tilde{\mathbf{k}}^{\text{B}} \mathbf{grad} T^{(0)}) = 0. \quad (41) $$

In this latter equation,

$$ \tilde{\mathbf{k}}^{\text{B}} = \mathbf{k}^{\text{eff}} + \frac{\gamma(T^{(0)}) L_{\text{sg}}}{\rho_i} \mathbf{D}^{\text{eff}}, \quad (42) $$

which appears as an apparent thermal conductivity of the snow which depends non-linearly on the temperature through $\gamma(T^{(0)})$. Our results show that this is valid if $[W_R] = O(\epsilon)$, i.e., typically for $\alpha$ values ranging from around $2 \times 10^{-7}$ to $2 \times 10^{-5}$. Finally, let us remark that (i) model B can be also seen as a particular case of model A when $\rho_v^{(0)}$ tends towards $\rho_{vs}^{(0)}(T^{(0)})$ by increasing $\alpha$, and (ii) the apparent thermal conductivity of the snow $\tilde{\mathbf{k}}^{\text{B}}$ can also be written as

$$ \tilde{\mathbf{k}}^{\text{B}} = \mathbf{k}^{\text{eff}} + \frac{\gamma(T^{(0)}) L_{\text{sg}} D_{\text{v}}}{\rho_i} \frac{\mathbf{D}^{\text{eff}}}{D_{\text{v}}} = \mathbf{k}^{\text{eff}} + k_{\text{dif}} \frac{\mathbf{D}^{\text{eff}}}{D_{\text{v}}}, \quad (43) $$

where $k_{\text{dif}} = \gamma(T^{(0)})L_{\text{sg}}D_{\text{v}}/\rho_i$ corresponds to an enhancement of the air thermal conductivity, as defined in Sect. 2.4. However, in that case, $\gamma(T^{(0)})$ depends on the macroscopic temperature $T^{(0)}$.

### 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

$$ (\rho C)^{\text{eff}} \frac{\partial T^{(0)}}{\partial t} - \text{div}(\mathbf{k}^{\text{C}} \mathbf{grad} T^{(0)}) = -L_{\text{sg}} \dot{\phi}, \quad (44) $$

$$ \phi \frac{\partial \rho_{\text{vs}}^{(0)}}{\partial t} - \text{div}(\mathbf{D}^{\text{C}} \mathbf{grad} \rho_{\text{vs}}^{(0)}(T^{(0)})) = \rho_i \dot{\phi}, \quad (45) $$

with

$$ \rho_{\text{v}}^{(0)} = \rho_{\text{vs}}^{(0)}(T^{(0)}), \quad (46) $$

where $(\rho C)^{\text{eff}}$ and $\mathbf{k}^{\text{C}}$ are the dimensionless effective thermal capacity and the effective dimensionless thermal conductivity, respectively, and defined as

$$ (\rho C)^{\text{eff}} = (1 - \phi)\rho_i C_i + \phi \rho_a C_a, \quad (47) $$

$$ \mathbf{k}^{\text{C}} = \frac{1}{|\Omega|} \left( \int_{\Omega_a} k_a(\mathbf{grad}s_a + \mathbf{I})d\Omega + \int_{\Omega_i} k_i(\mathbf{grad}s_i + \mathbf{I})d\Omega \right), \quad (48) $$

where $\mathbf{D}^{\text{C}}$ is the effective diffusion tensor defined as

$$ \mathbf{D}^{\text{C}} = \frac{1}{|\Omega|} \int_{\Omega_a} D_{\text{v}}(\mathbf{grad}(\mathbf{d} + s_a) + \mathbf{I}) d\Omega. \quad (49) $$

Moreover, $s_i$, $s_a$, and $\mathbf{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):

$$ \text{div}(k_i(\mathbf{grad}s_i + \mathbf{I})) = 0 \quad \text{in } \Omega_i, \quad (50) $$

$$ \text{div}(k_a(\mathbf{grad}s_a + \mathbf{I})) = 0 \quad \text{in } \Omega_a, \quad (51) $$

$$ s_i = s_a \quad \text{on } \Gamma, \quad (52) $$

$$ \begin{aligned} & (k_i(\mathbf{grad}s_i + \mathbf{I}) - k_a(\mathbf{grad}s_a + \mathbf{I})) \cdot \mathbf{n}_i \\ & = \frac{L_{\text{sg}}}{\rho_i} \alpha w_k \gamma(T^{(0)}) \mathbf{d} \quad \text{on } \Gamma, \end{aligned} \quad (53) $$

$$ \text{div}(D_{\text{v}}(\mathbf{grad}(\mathbf{d} + s_a) + \mathbf{I})) = 0 \quad \text{in } \Omega_a, \quad (54) $$

$$ (\mathbf{grad}(\mathbf{d} + s_a) + \mathbf{I}) \cdot \mathbf{n}_i = \alpha w_k \mathbf{d} \quad \text{on } \Gamma, \quad (55) $$

with

$$ \frac{1}{|\Omega|} \int_{\Omega} (s_a + s_i) d\Omega = 0, \quad (56) $$
$$ \frac{1}{|\Omega|} \int_{\Gamma} d\, d\Gamma = 0. \qquad (57) $$

The vectors $s_i$, $s_a$, and $d$ depend on the values of $\alpha$ and the temperature (notably through $\gamma(T^{(0)})$). As for model B, the macroscopic heat transfer equation (Eq. 44) can also be written as

$$ \left( (\rho C)^{\text{eff}} + \phi \gamma(T^{(0)}) \frac{L_{\text{sg}}}{\rho_i} \right) \frac{\partial T^{(0)}}{\partial t} - \text{div}(\tilde{\mathbf{k}}^{\text{C}} \text{grad} T^{(0)}) = 0, \quad (58) $$

where

$$ \tilde{\mathbf{k}}^{\text{C}} = \mathbf{k}^{\text{C}} + \frac{\gamma(T^{(0)})L_{\text{sg}}}{\rho_i} \mathbf{D}^{\text{C}} = \mathbf{k}^{\text{C}} + k_{\text{dif}} \frac{\mathbf{D}^{\text{C}}}{D_v}, \quad (59) $$

which appears as an apparent thermal conductivity of the snow. In that case $\tilde{\mathbf{k}}^{\text{C}}$ and $\mathbf{D}^{\text{C}}$ both depend on $k_a$, $k_i$, $k_{\text{dif}}$, and $\alpha$, and we have

$$ \tilde{\mathbf{k}}^{\text{C}} = \frac{1}{|\Omega|} \left( \int_{\Omega_a} (k_a + k_{\text{dif}})(\text{grad}(s_a + d) + \mathbf{I})d\Omega + \int_{\Omega_i} k_i(\text{grad} s_i + \mathbf{I})d\Omega \right). \quad (60) $$

This model C is valid for $[W_R] = O(1)$, i.e., for $\epsilon^{1/2} < [W_R] < \epsilon^{-1/2}$ and thus for $\alpha$ values in the range $(\epsilon^{1/2}D_{\text{vc}}/(lw_{\text{kc}})) < \alpha < (\epsilon^{-1/2}D_{\text{vc}}/(lw_{\text{kc}}))$ or, typically, $3 \times 10^{-5} < \alpha < 4 \times 10^{-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

$$ (\rho C)^{\text{eff}} \frac{\partial T^{(0)}}{\partial t} - \text{div}(\mathbf{k}^{\text{D}} \text{grad} T^{(0)}) = -L_{\text{sg}} \dot{\phi}, \quad (61) $$

$$ \phi \frac{\partial \rho_{\text{vs}}^{(0)}}{\partial t} - \text{div}(\mathbf{D}^{\text{D}} \text{grad} \rho_{\text{vs}}^{(0)}) = \rho_i \dot{\phi}, \quad (62) $$

$$ \rho_v^{(0)} = \rho_{vs}^{(0)}(T^{(0)}). \quad (63) $$

$(\rho C)^{\text{eff}}$ is the classical dimensionless effective thermal capacity. The macroscopic thermal conductivity tensor $\tilde{\mathbf{k}}^{\text{D}}$ and the macroscopic diffusion tensor $\mathbf{D}^{\text{D}}$ are defined as

$$ \tilde{\mathbf{k}}^{\text{D}} = \frac{1}{|\Omega|} \left( \int_{\Omega_a} k_a(\text{grad } \mathbf{r}_a + \mathbf{I})d\Omega + \int_{\Omega_i} k_i(\text{grad } \mathbf{r}_i + \mathbf{I})d\Omega \right), \quad (64) $$

$$ \mathbf{D}^{\text{D}} = \frac{1}{|\Omega|} \int_{\Omega_a} D_v(\text{grad } \mathbf{r}_a + \mathbf{I})d\Omega, \quad (65) $$

where $\mathbf{r}_a$ and $\mathbf{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:

$$ \text{div}(k_i(\text{grad } \mathbf{r}_i + \mathbf{I})) = 0 \quad \text{in } \Omega_i, \quad (66) $$

$$ \text{div}((k_a + k_{\text{dif}})(\text{grad } \mathbf{r}_a + \mathbf{I})) = 0 \quad \text{in } \Omega_a, \quad (67) $$

$$ \mathbf{r}_i = \mathbf{r}_a \quad \text{on } \Gamma, \quad (68) $$

$$ (k_i(\text{grad } \mathbf{r}_i + \mathbf{I}) - (k_a + k_{\text{dif}})(\text{grad } \mathbf{r}_a + \mathbf{I})) \cdot \mathbf{n}_i = 0 \quad \text{on } \Gamma, \quad (69) $$

$$ \frac{1}{|\Omega|} \int_{\Omega} (\mathbf{r}_a + \mathbf{r}_i)d\Omega = 0. \quad (70) $$

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

$$ \left( (\rho C)^{\text{eff}} + \phi \gamma(T^{(0)}) \frac{L_{\text{sg}}}{\rho_i} \right) \frac{\partial T^{(0)}}{\partial t} - \text{div}(\tilde{\mathbf{k}}^{\text{D}} \text{grad} T^{(0)}) = 0. \quad (71) $$

In this latter equation,

$$ \tilde{\mathbf{k}}^{\text{D}} = \mathbf{k}^{\text{D}} + \frac{\gamma(T^{(0)})L_{\text{sg}}}{\rho_i} \mathbf{D}^{\text{D}} = \mathbf{k}^{\text{D}} + k_{\text{dif}} \frac{\mathbf{D}^{\text{D}}}{D_v}, \quad (72) $$

which appears as an apparent thermal conductivity of the snow. In that case, $\tilde{\mathbf{k}}^{\text{D}}$ and $\mathbf{D}^{\text{D}}$ both depend on $k_a$, $k_i$, and $k_{\text{dif}}$, and we have

$$ \tilde{\mathbf{k}}^{\text{D}} = \frac{1}{|\Omega|} \left( \int_{\Omega_a} (k_a + k_{\text{dif}})(\text{grad } \mathbf{r}_a + \mathbf{I})d\Omega + \int_{\Omega_i} k_i(\text{grad } \mathbf{r}_i + \mathbf{I})d\Omega \right). \quad (73) $$

This model D corresponds to the one derived by Moyne et al. (1988), assuming that $\rho_v = \rho_{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 $\alpha \approx 10^{-2}$. In that case, we show that this model is valid for $[W_R] = O(\epsilon^{-1})$ or $O(\epsilon^{-2})$, i.e., typically for $\alpha$ values ranging from around $5 \times 10^{-3}$ 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 $\alpha$. 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><img src="image_4.png" /><figcaption>**Figure 3.** Definition of the three different macroscopic models and their domain of validity with respect to $\alpha$. The value of $\alpha_T$ was estimated based on the characteristic values given in Table 1.</figcaption></figure>

A first important outcome is that the hypothesis $\rho_v = \rho_{vs}(T)$, which is often made, appears to be a good approximation for $\alpha$ values larger than $10^{-6}$, as seen for models B, C, and D. The asymptotic analysis shows that in the range of $\alpha$ [$10^{-6}$, $\alpha_T$], this approximation is of the order of $O(\varepsilon)$ since $\rho_v^{(0)} = \rho_{vs}^{(0)}(T^{(0)})$, i.e., $\sigma = (\rho_v - \rho_{vs})/\rho_{vs} \approx O(\varepsilon)$. In the range [$\alpha_T$, 1], this approximation is of the order of $O(\varepsilon^2)$, since $\rho_v^{(0)} = \rho_{vs}^{(0)}(T^{(0)})$ and $\rho_v^{(1)} = \rho_{vs}^{(1)}(T^{(1)})$, i.e., $\sigma \approx O(\varepsilon^2)$. 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 $\tilde{k}^\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 $\alpha$. Classically, heat and vapor transport are driven by the effective thermal conductivity of snow $k^\text{eff}$ and the effective vapor diffusivity of snow $D^\text{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 $\rho_v$ tends towards $\rho_{vs}(T)$ at the macro-scale, leading to a simplification in the form of one heat transfer equation in which the apparent thermal conductivity $\tilde{k}^B$ can be easily expressed with respect to $k^\text{eff}$ ($k_i$, $k_a$, and microstructure) and $D^\text{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 $\alpha$ is not involved but instead driven by an apparent thermal conductivity $\tilde{k}^D$, which can be expressed with respect to the macroscopic thermal conductivity $k^D$ ($k_i$, $k_a$, $k_{\text{dif}}$, and microstructure) and to the macroscopic diffusion $D^D$ ($D_v$, $k_i$, $k_a$, $k_{\text{dif}}$, and microstructure), with the latter appearing as a thermodiffusion coefficient. Note that these parameters depend on the air thermal properties $k_a + k_{\text{dif}}$ that are enhanced by the phase change through $k_{\text{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 $\alpha_T$ estimated here at $3 \times 10^{-4}$, yet we recall that it can vary between $10^{-5}$ and $10^{-3}$, 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 $\tilde{k}^C$ can be expressed with respect to
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Illustration of the bilayer snowpack problem (a) at the macroscopic scale and (b) at the scale of a representative elementary volume (REV).</figcaption></figure>

the macroscopic thermal conductivity $k^C(k_i, k_a, k_{\text{dif}}, \alpha$, and microstructure) and the macroscopic diffusion $D^C(D_v, k_i, k_a, k_{\text{dif}}, \alpha$, and microstructure). Unlike the other models, both macroscopic parameters $k^C$ and $D^C$ depend on $\alpha$. These parameters tend towards $k^\text{eff}$ and $D^\text{eff}$ when $\alpha$ tends towards $10^{-5}$, thus recovering model B. They also tend towards $k^D$ and $D^D$ when $\alpha$ 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 $\alpha$ (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).

#### 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,

$$D^{\text{eff}} = \begin{pmatrix} D_{11}^{\text{eff}} & 0 \\ 0 & 0 \end{pmatrix} \qquad D_{11}^{\text{eff}} = \phi D_v, \qquad (74)$$

$$k^{\text{eff}} = \begin{pmatrix} k_{11}^{\text{eff}} & 0 \\ 0 & k_{22}^{\text{eff}} \end{pmatrix} \qquad k_{11}^{\text{eff}} = \phi k_a + (1-\phi)k_i, \\ k_{22}^{\text{eff}} = \frac{k_i k_a}{(1-\phi)k_a + \phi k_i}. \qquad (75)$$

Thus, it follows that

$$\tilde{k}_{11}^{\text{B}} = \phi(k_a + k_{\text{dif}}) + (1-\phi)k_i, \qquad \tilde{k}_{22}^{\text{B}} = \frac{k_i k_a}{(1-\phi)k_a + \phi k_i}. \qquad (76)$$

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

$$0 \le D^{\text{eff}} \le \phi D_v, \qquad (77)$$

and

$$\frac{k_i k_a}{(1-\phi)k_a + \phi k_i} \le k^{\text{eff}} \le \phi k_a + (1-\phi)k_i, \\ \frac{k_i k_a}{(1-\phi)k_a + \phi k_i} \le \tilde{k}^{\text{B}} \le \phi(k_a + k_{\text{dif}}) \\ + (1-\phi)k_i. \qquad (78)$$

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

$$D^{\text{D}} = \begin{pmatrix} D_{11}^{\text{D}} & 0 \\ 0 & D_{22}^{\text{D}} \end{pmatrix} \qquad D_{11}^{\text{D}} = \phi D_v, \\ D_{22}^{\text{D}} = \phi D_v \frac{k_i}{(1-\phi)(k_a + k_{\text{dif}}) + \phi k_i}, \qquad (79)$$

$$k^{\text{D}} = \begin{pmatrix} k_{11}^{\text{D}} & 0 \\ 0 & k_{22}^{\text{D}} \end{pmatrix} \qquad k_{11}^{\text{D}} = \phi k_a + (1-\phi)k_i, \\ k_{22}^{\text{D}} = \frac{k_i(k_a + (1-\phi)k_{\text{dif}})}{(1-\phi)(k_a + k_{\text{dif}}) + \phi k_i}. \qquad (80)$$

Thus, it follows that

$$k_{11}^{\text{D}} = \phi(k_a + k_{\text{dif}}) + (1-\phi)k_i, \\ k_{22}^{\text{D}} = \frac{k_i(k_a + k_{\text{dif}})}{(1-\phi)(k_a + k_{\text{dif}}) + \phi k_i}. \qquad (81)$$

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

$$\phi D_v \le D^D \le \phi D_v \frac{k_i}{(1-\phi)(k_a + k_{\text{dif}}) + \phi k_i}, \qquad (82)$$

and

$$\frac{k_i(k_a + (1-\phi)k_{\text{dif}})}{(1-\phi)(k_a + k_{\text{dif}}) + \phi k_i} \le k^D \le \phi k_a + (1-\phi)k_i, \\ \frac{k_i(k_a + k_{\text{dif}})}{(1-\phi)(k_a + k_{\text{dif}}) + \phi k_i} \le \tilde{k}^D \\ \le \phi(k_a + k_{\text{dif}}) + (1-\phi)k_i. \qquad (83)$$

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^{\text{eff}}$ and $D^{\text{D}}$ are always smaller than $D_v$ and $D^{\text{D}} > D^{\text{eff}}$, whatsoever the $\alpha$ value. Moreover, according to the definition of $D^{\text{eff}}$ (Eq. 74) and $D^{\text{D}}$ (Eq. 79), if a vertical macroscopic temperature gradient is applied along $e_2$, 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_2$, since $D_{22}^{\text{D}} \neq 0$, and thus a variation in the porosity along $e_2$.

## 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.

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

$$k_{\text{SC}}^{\text{eff}} = \frac{\beta + \sqrt{\beta^2 + 8k_i k_a}}{4} \quad \text{with } \beta = k_i(3(1-\phi) - 1) + k_a(3\phi - 1), \quad (84)$$

$$D_{\text{SC}}^{\text{eff}} = D_v \frac{(3\phi - 1)}{2}. \quad (85)$$

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

$$\tilde{k}_{\text{SC}}^{\text{B}} = k_{\text{SC}}^{\text{eff}} + k_{\text{dif}} \frac{D_{\text{SC}}^{\text{eff}}}{D_v}. \quad (86)$$

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

$$\tilde{k}_{\text{SC}}^{\text{D}} = \frac{\beta + \sqrt{\beta^2 + 8k_i(k_a + k_{\text{dif}})}}{4} \quad \text{with } \beta = k_i(3(1-\phi) - 1) + (k_a + k_{\text{dif}})(3\phi - 1). \quad (87)$$

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

$$D_{\text{SC}}^{\text{D}} = \phi D_v \frac{3\tilde{k}_{\text{SC}}^{\text{D}}}{(k_a + k_{\text{dif}}) + 2\tilde{k}_{\text{SC}}^{\text{D}}}. \quad (88)$$

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_{\text{SC}}^{\text{eff}}$, $\tilde{k}_{\text{SC}}^{\text{B}}$, and $\tilde{k}_{\text{SC}}^{\text{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_{\text{SC}}^{\text{eff}}/D_v$ is overall much smaller than $D_{\text{SC}}^{\text{D}}/D_v$ and evolves linearly from 0.25 to 1 when porosity varies from 0.5 to 1 (Fig. 6). In contrast, $D_{\text{SC}}^{\text{D}}/D_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_{\text{SC}}^{\text{D}}/D_v$ come from the contribution of the heat conduction, through $\tilde{k}_{\text{SC}}^{\text{D}}$, and of the latent heat, through $k_{\text{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_{\text{dif}}$ and $k_a + k_{\text{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_{\text{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_{\text{dif}}$ evolves in the same way as $k_{\text{dif}}$ (non-linear), but the values are increased by $k_a$.

Keeping the above considerations in mind, the evolution of $k_{\text{SC}}^{\text{eff}}$, $\tilde{k}_{\text{SC}}^{\text{B}}$, and $\tilde{k}_{\text{SC}}^{\text{D}}$ with temperature is presented in Fig. 5b, d, and f. For $k_{\text{SC}}^{\text{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 $\tilde{k}_{\text{SC}}^{\text{B}}$ and $\tilde{k}_{\text{SC}}^{\text{D}}$ with temperature is more complex as the impact of $k_{\text{dif}}$ superimposes. They show a non-linear evolution with temperature, with an evolution similar to $k_{\text{SC}}^{\text{eff}}$ for the lower temperatures and transitioning to an exponential increase for the higher
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Evolution of the SC estimates of the thermal conductivities $k_{\text{SC}}^{\text{eff}}$, $\tilde{k}_{\text{SC}}^{\text{B}}$, and $\tilde{k}_{\text{SC}}^{\text{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.</figcaption></figure>

temperatures, where the latter is driven by $k_{\text{dif}}$. We see that this non-linearity is even more important for $\tilde{k}_{\text{SC}}^{\text{D}}$ than for $\tilde{k}_{\text{SC}}^{\text{B}}$, as $k_{\text{dif}}$ appears in the definition of $\tilde{k}_{\text{SC}}^{\text{D}}$ several times. Finally, estimates of diffusion coefficient $\tilde{k}_{\text{SC}}^{\text{D}}$ show a slight influence of temperature through $k_i$, $k_a$, and $k_{\text{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_{\text{top}}$ and at the bottom $T_{\text{bottom}}$ are imposed, and $T_{\text{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
<figure><img src="image_7.png" /></figure>

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

<figure><img src="image_8.png" /></figure>

**Figure 7.** Evolution of the thermal conductivity with temperature for $k_i$ (Huang et al., 2013), $k_a$ (Haynes, 2016), $k_{\text{dif}}$, and $k_a + k_{\text{dif}}$.

the lateral sides of the snow layer. Simulations were run in a steady state.

At the pore scale, the snow layer consists of 200 periodic cells of $0.5 \times 0.5 \text{ mm}^2$; 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 \text{ kg m}^{-3}$. The heat and the mass transfer is described by the set of Eqs. (1)–(12), where $T_i$, $T_a$, and $\rho_v$ are the unknowns. This set of equations was numerically solved using the material parameter values presented in Table 1 and for different $\alpha$ values in the range of $10^{-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 \text{ °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^{(0)}$ and $\rho_v^{(0)}$ are the macroscopic unknowns. These macroscopic

<figure><img src="image_9.png" /></figure>

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

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 \text{ m}^{-1}$, respectively. The effective properties $k^{\text{eff}}$ and $D^{\text{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^{\text{eff}} = 0.04243 \text{ W m}^{-1} \text{ K}^{-1}$ and $D^{\text{eff}} = 1.156 \times 10^{-5} \text{ m}^2 \text{ s}^{-1}$. The apparent thermal conductivity $\tilde{k}^B$ was analytically deduced using Eq. (41). Its value depends on the temperature through the term $k_{\text{dif}}(T)$. $\tilde{k}^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_{\text{dif}}(T)$ for $\alpha$ values in the range $10^{-6}$ to 1. Finally, $\tilde{k}^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 $2.01 \times 10^{-5} \text{ m}^2 \text{ s}^{-1}$.

Figure 9 presents the evolution of $k^{\text{eff}}$, $\tilde{k}^B$, and $\tilde{k}^D$ with temperature. As expected, $\tilde{k}^B$ and $\tilde{k}^D$ evolve nonlinearly with $T$. To perform the simulations, the computed values of $\tilde{k}^D$ were fitted by the following relation: $\tilde{k}^D = 46.064(T/273)^4 - 156.05(T/273)^3 + 198.7(T/273)^2 - 112.68(T/273) + 24.045$ (blue line in Fig. 9). Figure 10 shows the evolution of the dimensionless diffusion coefficients $D^{\text{eff}}/D_v$, $D^C/D_v$, and $D^D/D_v$ and of the macroscopic
<figure><img src="image_10.png" /><figcaption>Figure 9. Evolution of the thermal conductivities $k^{\text{eff}}$, $\tilde{k}^{\text{B}}$, and $\tilde{k}^{\text{D}}$ with temperature. For $\tilde{k}^{\text{D}}$, the blue dots represent the numerical estimates over the REV, and the blue line is the fit.</figcaption></figure>

thermal conductivities $k^{\text{eff}}$, $k^{\text{C}}$, and $k^{\text{D}}$ with respect to $\alpha$ and for two temperatures of 270 and 250 K. As expected, only the parameters of model C ($D^{\text{C}}/D_{\text{v}}$ and $k^{\text{C}}$) vary with $\alpha$ and ensure a continuous transition between the parameters of model A ($D^{\text{eff}}/D_{\text{v}}$ and $k^{\text{eff}}$) and the ones of model D ($D^{\text{D}}/D_{\text{v}}$ and $k^{\text{D}}$). When fitting the numerical estimates, such a transition can be described by a simple function as follows:

$$ \frac{D^{\text{C}}(\alpha) - D^{\text{eff}}}{D^{\text{D}} - D^{\text{eff}}} = \frac{k^{\text{C}}(\alpha) - k^{\text{eff}}}{k^{\text{D}} - k^{\text{eff}}} = \frac{\tilde{k}^{\text{C}}(\alpha) - \tilde{k}^{\text{B}}}{\tilde{k}^{\text{D}} - \tilde{k}^{\text{B}}} = \frac{A\alpha}{1 + A\alpha}, \quad (89) $$

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).

### 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{align} \langle T \rangle &= \frac{1}{\Omega} \left( \int_{\Omega_i} T_i d\Omega + \int_{\Omega_a} T_a d\Omega \right), \nonumber \\ \langle \rho_v \rangle &= \frac{1}{\Omega_a} \int_{\Omega_a} \rho_v d\Omega, \quad \langle \rho_{vs}(T) \rangle = \frac{1}{\Omega_a} \int_{\Omega_a} \rho_{vs}(T_a) d\Omega, \nonumber \\ \langle \dot{\phi} \rangle &= \frac{1}{\Omega} \int_{\Gamma} w_n d\Gamma. \tag{90} \end{align} $$

Thus, we compare the vertical profiles of the pore-scale variables $\langle T \rangle$, $\langle \rho_v \rangle$, $\langle \rho_{vs} \rangle$, and $\langle \dot{\phi} \rangle$ with the vertical profiles of the macroscopic variables $T^{(0)}$, $\rho_v^{(0)}$, $\rho_{vs}^{(0)}(T^{(0)})$, and $\dot{\phi}$. As the obtained simulated temperature profiles were close to each other, to ease the comparison, we also use the temperature

deviation $\Delta T$, which represents the deviation of the simulated temperature profile from the linear temperature profile imposed by $T_{\text{top}}$ and $T_{\text{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 $\rho_v - \rho_{vs}(T)$. Figure 12 shows the vertical profiles of $\Delta T$, of $\rho_v - \rho_{vs}(T)$, and of $\dot{\phi}$ from the pore-scale simulations (dots) and the macroscopic models (lines), considering a temperature gradient of 100 and 500 K m$^{-1}$. For the pore-scale simulations, values of $\alpha$ from $10^{-9}$ to 1 were used. For the macroscopic models, results are only shown in their domain of validity with respect to $\alpha$. To further highlight the impact of $\alpha$, Fig. 13 presents the evolution of $T$, $\rho$, and $\rho_{vs}$ with $\alpha$ for the specific middle cell of the snow layer, which is here the 100th cell from the bottom ($x = l/2$ and $y = 100l - l/2$). Again, pore-scale simulations (dots) and the macroscopic models (line) are compared.

We describe first the main features observed in the pore-scale simulations. All the variables show an impact of the $\alpha$ value. The temperature deviation $\Delta 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 $\alpha$ and with the temperature gradient. This is also reflected in the temperature of the middle cell that overall increases with increasing $\alpha$ (Fig. 13a and b). This increase is not uniform, and two plateaus are observed – one between $10^{-6} \le \alpha \le \alpha_{\text{T}}$ and the other one between $10^{-1} \le \alpha \le 1$. The largest $\Delta T$ value is reached in the center of the snow layer and is around 0.4 K at 100 K m$^{-1}$ and 4 K at 500 K m$^{-1}$. Looking at the lower part of the layer, negative $\Delta T$ values can be found for $\alpha \le \alpha_{\text{T}} \sim 3 \times 10^{-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 $\rho_v - \rho_{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 $\alpha$ values when phase changes are very limited. With increasing $\alpha$ 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 $\rho_v \approx \rho_{vs}(T)$ for $10^{-6} \le \alpha \le 1$. Similar to the temperature, two plateaus are shown where $\rho_v - \rho_{vs}(T)$ evolve little with $\alpha$. All of these results are consistent with our theoretical analysis presented in Sect. 2. Last, the vertical profiles of $\dot{\phi}$ 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 $\alpha$ increases, those transitions become sharper and sharper, like a front. For $\alpha_{\text{T}} \le \alpha \le 1$, most values become negative, indicating overall deposition in the snow layer. A sublimation zone is still visible at the bottom of the
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Evolution of the dimensionless diffusion coefficients $D^{\text{eff}}/D_v$, $D^C/D_v$, and $D^D/D_v$ and of the macroscopic thermal conductivities $k^{\text{eff}}$, $k^C$, and $k^D$ with respect to $\alpha$ 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.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Simplified example of the transition from the $x$ to $\Delta x$ notation in a concave (red) and a convex (blue) case.</figcaption></figure>

snow layer, but its thickness is typically of the order of a few REV or smaller. Finally, as the difference $\rho_v - \rho_{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 $\alpha = 10^{-6}$. For $100 \text{ K m}^{-1}$, a value of $5.9 \times 10^{-13} \text{ m s}^{-1}$ and of $-2.7 \times 10^{-11} \text{ m s}^{-1}$ is found in the middle and bottom cell, respectively. For $500 \text{ K m}^{-1}$, a value of $4.5 \times 10^{-12} \text{ m s}^{-1}$ and of $-1.1 \times 10^{-10} \text{ m s}^{-1}$ 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 $\alpha$. For $\alpha \le 10^{-5}$, model A reproduces precisely all the features shown at the pore scale. Models B and D are independent of $\alpha$ and provide one estimate of the temperature, and thus the water vapor density, for all the $\alpha$ 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 $10^{-7} \le \alpha \le 10^{-5}$ for model B and $10^{-2} \le \alpha \le 1$ for model D (Fig. 13). Model C, which depends on $\alpha$, allows us to reproduce the main features shown at the pore scale for

$\alpha$ values in the range $10^{-5} \le \alpha \le 1$. This model predicts $\Delta T$ values higher than the pore-scale simulations for a given $\alpha$. Moreover, models B, C, and D predict only the deposition in snow, with negative $\dot{\phi}$ 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 $\rho_v = \rho_{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).

## 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
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Vertical profiles of $\Delta T$, $\rho_v - \rho_{vs}(T)$, and $\dot{\phi}$ 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$^{-1}$ and for different values of $\alpha$. $\Delta T$ represents the deviation of the temperature profile from a linear temperature profile. Predictions of $\rho_v - \rho_{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 $\alpha = 1$ (dotted yellow lines).</figcaption></figure>

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$^{-1}$ 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$^{-1}$ 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 den-
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Temperature and water vapor density in the middle of the snow layer as a function of $\alpha$, 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^{-1}$. The models are only shown for the $\alpha$ values within their domain of validity. Values of saturation water vapor density $\rho_{vs}$ from the pore-scale simulations and from model A are also presented.</figcaption></figure>

<table><caption>**Table 2.** Overview of the experimental settings used in the simulations.</caption><thead><tr><th>Experiment</th><th>Bouvet A</th><th>Bouvet B</th><th>Kamata</th></tr></thead><tbody><tr><td>Initial density (kg m<sup>-3</sup>)</td><td>210</td><td>287</td><td>165</td></tr><tr><td>Snow layer height (cm)</td><td>13.5</td><td>7.7</td><td>10</td></tr><tr><td>Temperature gradient (K m<sup>-1</sup>)</td><td>93</td><td>103</td><td>530</td></tr><tr><td>Snow base temperature (°C)</td><td>-3.1</td><td>-6.5</td><td>-12</td></tr><tr><td>Snow surface temperature (°C)</td><td>-15.6</td><td>-14.5</td><td>-65</td></tr><tr><td>Duration (days)</td><td>20</td><td>28</td><td>5.5</td></tr></tbody></table>

sity 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^{-1}$ 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^{\text{eff}}$, the effective vapor diffusivity $D^{\text{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^{\text{eff}}$ and $D^{\text{eff}}$ were computed with the GeoDict software (Thoe men 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^{\text{eff}}$. The obtained 3D tensors of both properties show negligible non-diagonal terms. In the following, we refer to $k^{\text{eff}}$ and $D^{\text{eff}}$ as the average of the diagonal terms of the tensors.

Figure 14 presents the results of the image-based computations of $k^{\text{eff}}$, $D^{\text{eff}}$, and SSA for the experiment of Bouvet A and Bouvet B. To compare them, we show the estimates of $k^{\text{eff}}$ and $D^{\text{eff}}$ by the SC model presented in Sect. 3.2, the density-based parameterizations of $k^{\text{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 $\text{SSA}_{\text{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^{\text{eff}}$ is estimated with the parameterization of Calonne et al. (2011). In the second set called SC, $k^{\text{eff}}$ is given by the self-consistent estimates. In both sets, $D^{\text{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^{\text{eff}}$, $D^{\text{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.

### 4.2.2 Models B and D

Models B and D only involve the apparent thermal conductivities of snow $\tilde{k}^B$ and $\tilde{k}^D$, respectively. Estimating $\tilde{k}^B$ comes down to estimating $k^{\text{eff}}$ and $D^{\text{eff}}$, as it is defined as $\tilde{k}^B = k^{\text{eff}} + k_{\text{dif}}D^{\text{eff}}/D_v$. For that reason, we use the same estimates of $k^{\text{eff}}$ and $D^{\text{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 $\tilde{k}_{\text{Calonne}}^B$ and $\tilde{k}_{\text{SC}}^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).

For model D, $\tilde{k}^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 $\tilde{k}^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 $\tilde{k}^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 $\tilde{k}_{\text{SC}}^D$ presented in Sect. 3.2. For each experiment, a fit was performed on the computed data and is referred to as $\tilde{k}_{\text{Fit}}^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 $\tilde{k}_{\text{Fit}}^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 $\tilde{k}_{\text{SC}}^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
<figure><img src="image_15.png" /><figcaption>**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).</figcaption></figure>

<table><caption>**Table 3.** Summary of the effective parameters used in the simulations.</caption><tbody><tr><td rowspan="2">Model A</td><td>SC set</td><td><i>k</i><sub>SC</sub><sup>eff</sup>, <i>D</i><sub>SC</sub><sup>eff</sup>, SSA<sub>Fit</sub>(t)</td></tr><tr><td>Calonne set</td><td><i>k</i><sub>Calonne</sub><sup>eff</sup>, <i>D</i><sub>SC</sub><sup>eff</sup>, SSA<sub>Fit</sub>(t)</td></tr><tr><td rowspan="2">Model B</td><td>SC set</td><td><i>k&#x0303;</i><sub>SC</sub><sup>B</sup> = <i>k</i><sub>SC</sub><sup>eff</sup> + <i>k</i><sub>dif</sub><i>D</i><sub>SC</sub><sup>eff</sup>/<i>D</i><sub>v</sub>, <i>D</i><sub>SC</sub><sup>eff</sup></td></tr><tr><td>Calonne set</td><td><i>k&#x0303;</i><sub>Calonne</sub><sup>B</sup> = <i>k</i><sub>Calonne</sub><sup>eff</sup> + <i>k</i><sub>dif</sub><i>D</i><sub>SC</sub><sup>eff</sup>/<i>D</i><sub>v</sub>, <i>D</i><sub>SC</sub><sup>eff</sup></td></tr><tr><td rowspan="2">Model C</td><td>SC set</td><td><i>k&#x0303;</i><sub>SC</sub><sup>C</sup>(<i>&alpha;</i>) and <i>D</i><sub>SC</sub><sup>C</sup>(<i>&alpha;</i>) from Eq. (89)</td></tr><tr><td>Calonne set</td><td><i>k&#x0303;</i><sub>Calonne</sub><sup>C</sup>(<i>&alpha;</i>) and <i>D</i><sub>Calonne</sub><sup>C</sup>(<i>&alpha;</i>) from Eq. (89)</td></tr><tr><td rowspan="2">Model D</td><td>SC set</td><td><i>k&#x0303;</i><sub>SC</sub><sup>D</sup>, <i>D</i><sub>SC</sub><sup>D</sup></td></tr><tr><td>Fit set</td><td><i>k&#x0303;</i><sub>Fit</sub><sup>D</sup>, <i>D</i><sub>SC</sub><sup>D</sup></td></tr></tbody></table>

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 mea-

surements 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 $\Delta T$ with models A, B, and D for the three experiments, taking different $\alpha$ values from $10^{-9}$ to $10^{-4}$ for model A. The $\Delta 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 $\Delta 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
<figure><img src="image_16.png" /><figcaption>**Figure 15.** The thermal conductivity estimates for model B ($\tilde{k}_{\text{SC}}^{\text{B}}$ and $\tilde{k}_{\text{Calonne}}^{\text{B}}$) and for model D ($\tilde{k}_{\text{SC}}^{\text{D}}$ and $\tilde{k}_{\text{Fit}}^{\text{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 $\tilde{k}_{\text{Fit}}^{\text{D}}$ are shown by blue dots. For Bouvet A, $\tilde{k}_{\text{Fit}}^{\text{D}}(T) = 5.4485 \times 10^{-9} \times T^4 - 4.8119 \times 10^{-6} \times T^3 + 1.5965 \times 10^{-3} \times T^2 - 2.3581 \times 10^{-1} \times T + 1.3195$. For Bouvet B, $\tilde{k}_{\text{Fit}}^{\text{D}}(T) = 6.0212 \times 10^{-9} \times T^4 - 5.2974 \times 10^{-6} \times T^3 + 1.7523 \times 10^{-3} \times T^2 - 2.5868 \times 10^{-1} \times T + 14.6338$. For Kamata, $\tilde{k}_{\text{Fit}}^{\text{D}}(T) = 5.1386 \times 10^{-9} \times T^4 - 4.5612 \times 10^{-6} \times T^3 + 1.5206 \times 10^{-3} \times T^2 - 2.2553 \times 10^{-1} \times T + 12.6279$.</figcaption></figure>

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 $\Delta T$ value in the upper part of the layer in Bouvet A is attributed to the temperature sensor uncertainty (Bouvet et al., 2023).

Looking at the models, the main observation is that they all underestimate $\Delta T$. In more detail, model A predicts negative $\Delta 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 $\alpha$, the positive values of $\Delta 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 $\Delta T$ predicted by model A is reached for the highest $\alpha = 10^{-4}$ 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 $\Delta T$ profile valid over their domain of validity, namely $10^{-6} \le \alpha \le \alpha_T$ and $\alpha_T \le \alpha \le 1$, respectively. The profile shape is in agreement with the measurements, showing only positive values throughout the layer. $\Delta 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 $\tilde{k}_{\text{SC}}^{\text{B}}$ and $\tilde{k}_{\text{SC}}^{\text{D}}$, compared to $\tilde{k}_{\text{Calonne}}^{\text{B}}$ and $\tilde{k}_{\text{Fit}}^{\text{D}}$, respectively, lead to reduced overall heat conduction through snow and allow for higher $\Delta 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 $\Delta 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$^{-3}$, respectively. For the same temperature gradient, the higher the snow density, the higher the heat conduction through snow and the lower the $\Delta T$. For example, in lighter snow (Bouvet A), a maximum $\Delta 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
<figure><img src="image_17.png" /><figcaption>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.</figcaption></figure>

$$h_{\text{air gap}} = H \frac{\bar{\phi} t_{\text{exp}}}{(\bar{\phi} t_{\text{exp}} + \phi_{\text{init}} - 1)},$$
$$\text{with } \bar{\phi} = \frac{1}{H} \int_{0}^{H} \dot{\phi}(z) dz$$
$$= -\frac{1}{H} \int_{0}^{H} \frac{1}{\rho_{i}} \frac{\partial}{\partial z} \left( \frac{\partial D \rho_{\text{vs}}^{(0)}(T)}{\partial z} \right) dz, \quad (91)$$

with $h_{\text{air gap}}$ as the height of the air gap (m), $H$ as the total height of the snow layer (m), $\phi_{\text{init}}$ as the initial porosity (–), and $t_{\text{exp}}$ as the total duration of the experiment (s) (Table 2). $D$ is the diffusivity coefficient and corresponds to $D_{\text{SC}}^{\text{eff}}$ for model B and $D_{\text{SC}}^{\text{D}}$ for model D. For model C, $D$ corresponds to $D^{\text{C}}(\alpha)$ 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),

$$\frac{h_{\text{air gap}}^{\text{C}} - h_{\text{air gap}}^{\text{B}}}{h_{\text{air gap}}^{\text{D}} - h_{\text{air gap}}^{\text{B}}} = \frac{A\alpha}{1 + A\alpha}, \quad (92)$$

with $A = 1200$.

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 an 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 α = 10⁻⁶ 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 α = 10⁻⁵ 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
<figure><img src="image_18.png" /><figcaption>**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 $\alpha$ values from $10^{-5}$ to $10^{-9}$ 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 $\alpha$, 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.</figcaption></figure>

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 $\alpha = 10^{-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 \text{ kg m}^{-3}$. 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 im-

pact 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 $\phi$ simulations on the simplified microstructure for the temperature gradient of $500 \text{ K m}^{-1}$ 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 \text{ kg m}^{-3}$ is found for the first 2.5 cm, in agreement with the measured one of $152 \text{ kg m}^{-3}$. Coming back to the air gap comparisons, model
A predicts the air gap estimated for the Kamata experiment well when the highest $\alpha$ value is considered. At $\alpha = 10^{-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 $\alpha$ in the range [$10^{-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 $\alpha_T$, which equals $\approx 3 \times 10^{-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 $\alpha$ values; and (ii) to show that the hypothesis $\rho_v = \rho_{vs}(T)$, which is often made, is a good approximation for $\alpha$ values larger than $10^{-6}$.

At the macroscopic scale, model A (Calonne et al., 2014b), valid for $\alpha$ values in the range [$10^{-10}$, $10^{-5}$], 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 $\alpha$. 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 $\rho_v = \rho_{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 $\alpha$ and are able to describe the observed plateaus in

a limited range of $\alpha$. These two models can be seen as particular cases of models A and C. Consequently, in practice, models A and C, which depend on $\alpha$, are sufficient to describe the macroscopic heat and mass transfer through the snowpack for $\alpha$ values in the range $10^{-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 $\Delta 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 $\pm 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 $\pm 0.5$ °C. In both cases, these uncertainties are smaller than the discrepancies between the measured and modeled $\Delta 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^{\text{eff}}$, $D^{\text{eff}}$, and the SSA for model A; on $\tilde{k}^B$ for model B; and on $\tilde{k}^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^{\text{eff}}$ by 10%, a $\Delta T$ of 0.096 K and an air gap of 1.17 mm is simulated for $\alpha = 10^{-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 $\alpha$. 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 $\alpha$ (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 $\alpha$ 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 $\Delta 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_{\rho v}$) should be expressed at the pore scale as follows:

$$ J_T = -k_a \mathbf{grad} T_a + \rho_a C_a v_a T_a \quad \text{in } \Omega_a, \\ \text{and} \quad J_{\rho v} = -D_v \mathbf{grad} \rho_v + v_a \rho_v \quad \text{in } \Omega_a, \qquad (93) $$

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

$$ Ra = \frac{\rho_a \beta_T g (T_{\text{bottom}} - T_{\text{top}}) H K}{((\mu_a k^{\text{eff}}) / (\rho_a C_a))}, \qquad (94) $$

where $H$ is the height of the snow layer, $K$ is the snow permeability, $g = 9.81 \text{ m s}^{-2}$ is the gravity, $\mu_a = 17.29 \times 10^{-6} \text{ Pa s}$ is the air viscosity, and $\beta_T = 0.0036 \text{ K}^{-1}$ 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^{\text{Tv}}$, and the effect of temperature gradient
on the vapor density flux, called the Soret effect, is characterized by the thermodiffusion coefficient $D^{\text{T}}$. Taking these effects into account, the temperature and water vapor flux can be expressed as follows:

$$ J_{\text{T}} = -k_{\text{a}}\text{grad}T_{\text{a}} - D^{\text{T}}\text{grad}\rho_{\text{v}} \quad \text{in } \Omega_{\text{a}}, \\ \text{and} \quad J_{\rho_{\text{v}}} = -D_{\text{v}}\text{grad}\rho_{\text{v}} - D^{\text{T}}\text{grad}T_{\text{a}} \quad \text{in } \Omega_{\text{a}}. \quad (95) $$

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_{\text{T}} = D^{\text{T}}/D_{\text{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_{\text{T}}$ coefficients lead to an increase the simulated temperature deviation $\Delta T$, and vice versa. For example, for $\alpha = 0.1$ and a temperature gradient of 500 K m$^{-1}$, a maximum value of $\Delta T$ of 6.5, 4.2, and 2.8 K is simulated for $S_{\text{T}}$ equal to $-2 \times 10^{-4}$, 0, and $2 \times 10^{-4}$, respectively. The Soret effect can also be introduced in model D, by replacing $k_{\text{a}}+k_{\text{diff}}$ with $k_{\text{a}}+k_{\text{diff}}+S_{\text{T}}D_{\text{v}}L_{\text{sg}}/\rho_{\text{i}}$. Using the self-consistent estimate of thermal conductivity (Eq. 87) and for values of $S_{\text{T}}$ at $-2 \times 10^{-4}$, 0, and $2 \times 10^{-4}$, the maximum simulated values of $\Delta 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 $\alpha$ ranging from $10^{-10}$ to 1. We showed that, at a transition value $\alpha_{\text{T}} \approx 3 \times 10^{-4}$, the mod-

eling 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 $\alpha$. Models A and C are sufficient to describe the heat and mass transfer in the whole range of $\alpha$ values between $10^{-10}$ and 1. For $\alpha$ values between $10^{-10}$ and $10^{-5}$, 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 $\alpha$. This model does not presume any assumption on the saturation of the vapor density. For $\alpha$ values between $10^{-5}$ and 1, model C consists of one temperature equation which involves $\alpha$, since the hypothesis $\rho_{\text{v}} = \rho_{\text{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 $\alpha$ values in the range from $10^{-7}$ to $10^{-5}$, whereas model D can be seen as a particular case of model C for $\alpha$ values in the range from $10^{-3}$ 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 $\Delta 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 $\alpha$ value of $10^{-5}$; 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.
*Data availability.* The theoretical development of the macroscopic equivalent descriptions are available in the Supplement.

*Supplement.* The supplement related to this article is available online at: https://doi.org/10.5194/tc-18-4285-2024-supplement.

*Author contributions.* CG, NC, and LB conducted the derivation of the macroscopic models. LB applied the models to the experimental data. The analyses and interpretations were carried out by LB, FF, NC, and CG. LB and CG prepared the paper, with contributions from all co-authors.

*Competing interests.* The contact author has declared that none of the authors has any competing interests.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

*Acknowledgements.* We warmly thank the editor, Jürg Schweizer, and the reviewers, Tien Dung Le and the anonymous reviewer, for their insightful comments which helped to significantly improve the quality of the paper.

*Financial support.* The 3SR lab is part of the Labex Tec 21 (Investissements d'Avenir; grant no. ANR-11-LABX-0030). CNRM/CEN is part of Labex OSUG@2020 (Investissements d'Avenir; grant no. ANR-10-LABX-0056). This research has been supported by the Agence Nationale de la Recherche through the MiMESis-3D ANR project (grant no. ANR-19-CE01-0009).

*Review statement.* This paper was edited by Jürg Schweizer and reviewed by Tien Dung Le and one anonymous referee.

## References

Albert, M. R. and McGilvary, W. R.: Thermal effects due to air flow and vapor transport in dry snow, J. Glaciol., 38, 273–281, https://doi.org/10.3189/S0022143000003683, 1992.

Anderson, E. A.: A point energy and mass balance model of a snow cover, Tech. Rep., Office of Hydrology – National Weather Service, https://repository.library.noaa.gov/view/noaa/6392 (last access: 16 July 2024), 1976.

Auriault, J.-L.: Heterogeneous medium. Is an equivalent description possible?, Int. J. Eng. Sci., 29, 785–795, https://doi.org/10.1016/0020-7225(91)90001-J, 1991.

Auriault, J.-L., Boutin, C., and Geindreau, C.: Homogenization of coupled phenomena in heterogenous media, Wiley-ISTE, London, https://doi.org/10.1002/9780470612033, 2009.

Barrett, J. W., Garcke, H., and Nürnberg, R.: Numerical computations of faceted pattern formation in snow crystal growth, Phys. Rev. E, 86, 011604, https://doi.org/10.1103/PhysRevE.86.011604, 2012.

Bensoussan, A., Lions, J.-L., and Papanicolaou, G.: Asymptotic Analysis for periodic structures, North Holland, https://www.ams.org/books/chel/374/chel374-endmatter.pdf (last access: 16 July 2024), 1978.

Bouvet, L., Calonne, N., Flin, F., and Geindreau, C.: Snow equilibrium metamorphism described by a phase-field model applicable on micro-tomographic images: Prediction of microstructural and transport properties, J. Adv. Model. Earth Sy., 14, e2022MS002998, https://doi.org/10.1029/2022MS002998, 2022.

Bouvet, L., Calonne, N., Flin, F., and Geindreau, C.: Heterogeneous grain growth and vertical mass transfer within a snow layer under a temperature gradient, The Cryosphere, 17, 3553–3573, https://doi.org/10.5194/tc-17-3553-2023, 2023.

Bruggeman, V. D.: Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen, Annalen der Physik, 416, 636–664, https://doi.org/10.1002/andp.19354160705, 1935.

Brzoska, J.-B., Flin, F., and Barckicke, J.: Explicit iterative computation of diffusive vapour field in the 3-D snow matrix: preliminary results for low flux metamorphism, Ann. Glaciol., 48, 13–18, https://doi.org/10.3189/172756408784700798, 2008.

Budiansky, B.: On the elastic moduli of some heterogeneous materials, J. Mech. Phys. Solids, 13, 223–227, https://doi.org/10.1016/0022-5096(65)90011-6, 1965.

Caldwell, D. R.: Measurement of negative thermal diffusion coefficients by observing the onset of thermaline convection, J. Phys. Chem., 77, 2004–2008, https://doi.org/10.1021/j100635a018, 1973.

Calonne, N., Flin, F., Morin, S., Lesaffre, B., Rolland du Roscoat, S., and Geindreau, C.: Numerical and experimental investigations of the effective thermal conductivity of snow, Geophys. Res. Lett., 38, L23501, https://doi.org/10.1029/2011GL049234, 2011.

Calonne, N., Geindreau, C., Flin, F., Morin, S., Lesaffre, B., Rolland du Roscoat, S., and Charrier, P.: 3-D image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy, The Cryosphere, 6, 939–951, https://doi.org/10.5194/tc-6-939-2012, 2012.

Calonne, N., Flin, F., Geindreau, C., Lesaffre, B., and Rolland du Roscoat, S.: Study of a temperature gradient metamorphism of snow from 3-D images: time evolution of microstructures, physical properties and their associated anisotropy, The Cryosphere, 8, 2255–2274, https://doi.org/10.5194/tc-8-2255-2014, 2014a.

Calonne, N., Geindreau, C., and Flin, F.: Macroscopic modeling for heat and water vapor transfer in dry snow by homogenization, J. Phys. Chem. B, 118, 13393–13403, https://doi.org/10.1021/jp5052535, 2014b.

Calonne, N., Geindreau, C., and Flin, F.: Macroscopic modeling of heat and water vapor transfer with phase change in dry snow based on an upscaling method: Influence of
air convection, J. Geophys. Res.-Earth, 120, 2476-2497,
https://doi.org/10.1002/2015JF003605, 2015.

Calonne, N., Milliancourt, L., Burr, A., Philip, A., Martin, C. L., Flin, F., and Geindreau, C.: Thermal Conductivity of Snow, Firn, and Porous Ice From 3-D Image-Based Computations, Geophys. Res. Lett., 46, 13079-13089,
https://doi.org/10.1029/2019GL085228, 2019.

Chapman, S. and Cowling, T. G.: The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases, Cambridge university press, ISBN 9780521408448, 1990.

Davarzani, H., Marcoux, M., and Quintard, M.: Theoretical predictions of the effective thermodiffusion coefficients in porous media, Int. J. Heat Mass Tran., 53, 1514-1528,
https://doi.org/10.1016/j.ijheatmasstransfer.2009.11.044, 2010.

de Quervain, M.: On the metamorphism of snow, in: Ice and Snow, edited by: Kingery, W., MIT Press, 377-390, https://www.umr-cnrm.fr/IMG/pdf/dequervain_1963_on_the_metamorphism_of_snow.pdf (last access: 16 July 2024), 1963.

Delaney, L., Houston, R., and Eagleton, L.: The rate of vaporization of water and ice, Chem. Eng. Sci., 19, 105-114,
https://doi.org/10.1016/0009-2509(64)85115-0, 1964.

Domine, F., Belke-Brea, M., Sarrazin, D., Arnaud, L., Barrere, M., and Poirier, M.: Soil moisture, wind speed and depth hoar formation in the Arctic snowpack, J. Glaciol., 64, 990-1002,
https://doi.org/10.1017/jog.2018.89, 2018.

Domine, F., Picard, G., Morin, S., Barrere, M., Madore, J.-B., and Langlois, A.: Major issues in simulating some arc-tic snowpack properties using current detailed snow physics models: consequences for the thermal regime and water budget of permafrost, J. Adv. Model. Earth Sy., 11, 34-44,
https://doi.org/10.1029/2018MS001445, 2019.

Dumont, M., Flin, F., Malinka, A., Brissaud, O., Hagenmuller, P., Lapalus, P., Lesaffre, B., Dufour, A., Calonne, N., Roland du Roscoat, S., and Ando, E.: Experimental and model-based investigation of the links between snow bidirectional reflectance and snow microstructure, The Cryosphere, 15, 3921-3948, https://doi.org/10.5194/tc-15-3921-2021, 2021.

Fierz, C., Armstrong, R. L., Durand, Y., Etchevers, P., Greene, E., McClung, D. M., Nishimura, K., Satyawali, P. K., and Sokratov, S. A.: The international classification for seasonal snow on the ground, IHP-VII Technical Documents in Hydrology, IACS Contribution no. 1, 83, https://unesdoc.unesco.org/ark:/48223/pf0000186462 (last access: 16 July 2024), 2009.

Flin, F.: Description physique des métamorphoses de la neige à partir d'images de microstructures 3D naturelles obtenues par microtomographie X, PhD thesis, Université Joseph Fourier, Grenoble, http://www.theses.fr/2004GRE10006 (last access: 16 July 2024), 2004.

Flin, F. and Brzoska, J.-B.: The temperature gradient metamorphism of snow: vapour diffusion model and application to tomographic images, Ann. Glaciol., 49, 17-21,
https://doi.org/10.3189/17275640878714834, 2008.

Flin, F., Brzoska, J.-B., Lesaffre, B., Coléou, C., and Pieritz, R. A.: Full three-dimensional modelling of curvature-dependent snow metamorphism: first results and comparison with experimental tomographic data, J. Phys. D Appl. Phys., 36, A49-A54,
https://doi.org/10.1088/0022-3727/36/10A/310, 2003.

Flin, F., Lesaffre, B., Dufour, A., Gillibert, L., Hasan, A., Rol-
land du Roscoat, S., Cabanes, S., and Puglièse, P.: On the com-
putations of specific surface area and specific grain contact area
from snow 3D images, in: Physics and Chemistry of Ice, edited
by: Furukawa, Y., Hokkaido University Press, Sapporo, Japan,
321–328, http://www.umr-cnrm.fr/IMG/pdf/flin_etal_2011_ssa_
sgca_publi shed_color.pdf (last access: 16 July 2024), 2011.

Fourteau, K., Domine, F., and Hagenmuller, P.: Impact of wa-
ter vapor diffusion and latent heat on the effective ther-
mal conductivity of snow, The Cryosphere, 15, 2739–2755,
https://doi.org/10.5194/tc-15-2739-2021, 2021a.

Fourteau, K., Domine, F., and Hagenmuller, P.: Macroscopic water vapor diffusion is not enhanced in snow, The Cryosphere, 15, 389-406, https://doi.org/10.5194/tc-15-389-2021, 2021b.

Furukawa, Y.: 25 – Snow and Ice Crystal Growth, in: Hand-
book of Crystal Growth, 2nd edn, edited by: Nishinaga,
T., Elsevier, Boston, 1061–1112, ISBN 978-0-444-56369-9,
https://doi.org/10.1016/B978-0-444-56369-9.00025-3, 2015.

Granger, R., Flin, F., Ludwig, W., Hammad, I., and Geindreau, C.: Orientation selective grain sublimation–deposition in snow under temperature gradient metamorphism observed with diffraction contrast tomography, The Cryosphere, 15, 4381–4398, https://doi.org/10.5194/tc-15-4381-2021, 2021.

Hansen, A. C. and Foslien, W. E.: A macroscale mixture the-
ory analysis of deposition and sublimation rates during heat
and mass transfer in dry snow, The Cryosphere, 9, 1857–1878,
https://doi.org/10.5194/tc-9-1857-2015, 2015.

Häussling Löwgren, B., Bergmann, J., and Alves-Filho,
O.: A Numerical Implementation of the Soret Ef-
fect in Drying Processes, ChemEngineering, 4, 13,
https://doi.org/10.3390/chemengineering4010013, 2020.

Haynes, W. M.: CRC handbook of chemistry and physics, CRC
press, https://doi.org/10.1201/9781315380476, 2016.

Hill, R.: A self-consistent mechanics of composite materials, J. Mech. Phys. Solids, 13, 213-222, https://doi.org/10.1016/0022-5096(65)90010-4, 1965.

Huang, W., Li, Z., Liu, X., Zhao, H., Guo, S., and Jia, Q.: Ef-
fective thermal conductivity of reservoir freshwater ice with
attention to high temperature, Ann. Glaciol., 54, 189-195,
https://doi.org/10.3189/2013AoG62A075, 2013.

Jafari, M., Sharma, V., and Lehning, M.: Convection of water vapour in snowpacks, J. Fluid Mech., 934, A38, https://doi.org/10.1017/jfm.2021.1146, 2022.

Jordan, R.: A one-dimensional temperature model for a snow cover: technical documentation for SNTHERM.89., Tech. Rep. 91-16, U.S. Army Cold Regions Research and Engineering Laboratory, http://hdl.handle.net/11681/11677 (last access: 16 July 2024), 1991.

Kaempfer, T. U. and Plapp, M.: Phase-field modeling of dry snow metamorphism, Phys. Rev. E, 79, 031502,
https://doi.org/10.1103/PhysRevE.79.031502, 2009.

Kamata, Y. and Sato, A.: Water-vapor transport in snow with high temperature gradient, in: Physics and Chemistry of Ice, edited by: Kuhs, W. F., 281-288, http://www.umr-cnrm.fr/IMG/pdf/kamata_and_sato_2007_high_tg.pdf (last access: 16 July 2024), 2007.

Kamata, Y., Sokratov, S. A., and Sato, A.: Temperature and temperature gradient dependence of snow recrystallization in depth hoar snow, in: Advances in Cold-Region Ther-
mal Engineering and Sciences, edited by: Hutter, K., Wang, Y., and Beer, H., Springer Berlin Heidelberg, 395–402, https://doi.org/10.1007/BFb0104197, 1999.

Krol, Q. and Löwe, H.: Analysis of local ice crystal growth in snow, J. Glaciol., 62, 378–390, https://doi.org/10.1017/jog.2016.32, 2016.

Legagneux, L., Taillandier, A.-S., and Domine, F.: Grain growth theories and the isothermal evolution of the specific surface area of snow, J. Appl. Phys., 95, 6175–6184, https://doi.org/10.1063/1.1710718, 2004.

Lehning, M., Bartelt, P., Brown, B., and Fierz, C.: A physical SNOWPACK model for the Swiss avalanche warning Part III: meteorological forcing, thin layer formation and evaluation, Cold Reg. Sci. Technol., 35, 169–184, https://doi.org/10.1016/S0165-232X(02)00072-1, 2002.

Libbrecht, K. G.: The physics of snow crystals, Rep. Prog. Phys., 68, 855–895, https://doi.org/10.1088/0034-4885/68/4/R03, 2005.

Libbrecht, K. G.: A taxonomy of snow crystal growth behaviors: 1. Using c-axis ice needles as seed crystals, arXiv [preprint], https://doi.org/10.48550/ARXIV.2109.00098, 2021.

Libbrecht, K. G. and Rickerby, M. E.: Measurements of surface attachment kinetics for faceted ice crystal growth, J. Cryst. Growth, 377, 1–8, https://doi.org/10.1016/j.jcysgro.2013.04.037, 2013.

Massman, W.: A review of the molecular diffusivities of H₂O, CO₂, CH₄, CO, O₃, SO₂, NH₃, N₂O, NO, and NO₂ in air, O₂ and N₂ near STP, Atmos. Environ., 32, 1111–1127, https://doi.org/10.1016/S1352-2310(97)00391-9, 1998.

Moyne, C., Batsale, J.-C., and Degiovanni, A.: Approche expérimentale et théorique de la conductivité thermique des milieux poreux humides – II. Théorie, Int. J. Heat Mass Tran., 31, 2319–2330, https://doi.org/10.1016/0017-9310(88)90163-9, 1988.

Persad, A. H. and Ward, C. A.: Expressions for the Evaporation and Condensation Coefficients in the Hertz-Knudsen Relation, Chem. Rev., 116, 7727–7767, https://doi.org/10.1021/acs.chemrev.5b00511, 2016.

Pinzer, B. R., Schneebeli, M., and Kaempfer, T. U.: Vapor flux and recrystallization during dry snow metamorphism under a steady temperature gradient as observed by time-lapse micro-tomography, The Cryosphere, 6, 1141–1155, https://doi.org/10.5194/tc-6-1141-2012, 2012.

Powers, D., O'Neill, K., and Colbeck, S. C.: Theory of Natural Convection in Snow, J. Geophys. Res., 90, 10641–10649, https://doi.org/10.1029/JD090iD06p10641, 1985.

Riche, F. and Schneebeli, M.: Thermal conductivity of snow measured by three independent methods and anisotropy considerations, The Cryosphere, 7, 217–227, https://doi.org/10.5194/tc-7-217-2013, 2013.

Sanchez-Palencia, E.: Non-homogeneous media and vibration theory, in: Lectures Notes in Physics, vol. 127, Springer-Verlag, Berlin, https://doi.org/10.1007/3-540-10000-8, 1980.

Srivastava, P., Mahajan, P., Satyawali, P., and Kumar, V.: Observation of temperature gradient metamorphism in snow by X-ray computed microtomography: measurement of microstructure parameters and simulation of linear elastic properties, Ann. Glaciol., 51, 73–82, https://doi.org/10.3189/172756410791386571, 2010.

Sturm, M. and Benson, C. S.: Vapor transport, grain growth and depth-hoar development in the subarctic snow, J. Glaciol., 43, 42–59, https://doi.org/10.3189/S0022143000002793, 1997.

Sturm, M. and Johnson, J.: Natural convection in the subarctic snow cover, J. Geophys. Res.-Sol. Ea., 96, 11657–11671, https://doi.org/10.1029/91JB00895, 1991.

Thoemen, H., Walther, T., and Wiegmann, A.: 3D simulation of macroscopic heat and mass transfer properties from the microstructure of wood fibre networks, Comp. Sci. Techn., 68, 608–616, https://doi.org/10.1016/j.compscitech.2007.10.014, 2008.

Torquato, S.: Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Springer, https://doi.org/10.1115/1.1483342, 2002.

Vetter, R., Sigg, S., Singer, H., Kadau, D., Herrmann, H., and Schneebeli, M.: Simulating isothermal aging of snow, EPL (Europhysics Letters), 89, 26001, https://doi.org/10.1209/0295-5075/89/26001, 2010.

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snow-pack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791, https://doi.org/10.5194/gmd-5-773-2012, 2012.

Wautier, A., Geindreau, C., and Flin, F.: Linking snow microstructure to its macroscopic elastic stiffness tensor: A numerical homogenization method and its application to 3-D images from X-ray tomography, Geophys. Res. Lett., 42, 8031–8041, https://doi.org/10.1002/2015GL065227, 2015.

Wiese, M.: Time-lapse tomography of mass fluxes and microstructural changes in snow, PhD thesis, ETH Zurich, https://doi.org/10.3929/ethz-b-000213853, 2017.

Yosida, Z., Oura, H., Kuroiwa, D., Huzioka, T., Kojima, K., Aoki, S., and Kinosita, S.: Physical Studies on Deposited Snow: I Thermal Properties, Tech. Rep. 7, Institute of Low Temperature Science, Hokkaido University, Sapporo, Japan, http://hdl.handle.net/2115/20216 (last access: 16 July 2024), 1955.