the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Macroscopic water vapor diffusion is not enhanced in snow

### Florent Domine

### Pascal Hagenmuller

Water vapor transport in dry snowpacks plays a significant role for snow metamorphism and the mass and energy balance of snowpacks. The molecular diffusion of water vapor in the interstitial pores is usually considered to be the main or only transport mechanism, and current detailed snow physics models therefore rely on the knowledge of the effective diffusion coefficient of water vapor in snow. Numerous previous studies have concluded that water vapor diffusion in snow is enhanced relative to that in air. Various field observations also indicate that for vapor transport in snow to be explained by diffusion alone, the effective diffusion coefficient should be larger than that in air. Here we show using theory and numerical simulations of idealized and measured snow microstructures that, although sublimation and deposition of water vapor onto snow crystal surfaces do enhance microscopic diffusion in the pore space, this effect is more than countered by the restriction of diffusion space due to ice. The interaction of water vapor with the ice results in water vapor diffusing more than inert molecules in snow but still less than in free air, regardless of the value of the sticking coefficient of water molecules on ice. Our results imply that processes other than diffusion play a predominant role in water vapor transport in dry snowpacks.

When a snowpack is submitted to a temperature gradient, macroscopic water vapor transfer occurs from the warmer to the colder parts of the snowpack in a process sometimes referred to as layer-to-layer vapor flux. This redistribution of mass plays a significant role in the evolution of the snowpack and its physical properties. In the absence of air convection in the snowpack, this macroscopic vapor flux results from the microscopic vapor diffusion occurring in the interstitial pores of snow and is impacted by water sublimation and deposition processes acting as sources and sinks of vapor at the ice–pore interface (Yosida et al., 1955; Colbeck, 1983). The physics at play in the pores is generally agreed upon, even though questions about the precise kinetics of the sublimation and deposition of water molecules onto ice surfaces in snow remain open (Legagneux and Domine, 2005; Pinzer et al., 2012; Calonne et al., 2014; Krol and Löwe, 2016). However, even for investigators assuming the same physics at the microscopic scale, the transition from the microscopic to the macroscopic scale remains a point of contention in the snow community (Giddings and LaChapelle, 1962; Colbeck, 1993; Pinzer et al., 2012; Hansen and Foslien, 2015; Shertzer and Adams, 2018; Hansen, 2019). Yet, a proper understanding of vapor transport in snow at the macroscopic scale is a prerequisite for accurate snowpack physical modeling.

There has notably been a long-standing controversy concerning the magnitude of the macroscopic diffusive fluxes transporting mass from one layer to another and in particular to determine whether they are larger than what would be observed in free air under similar macroscopic vapor gradients. The seminal study of Yosida et al. (1955) set out to measure in the laboratory the macroscopic vapor flux in a pile of snow subjected to a thermal gradient. Their results indicated that contrary to first expectations, the vapor flux was about 3 to 4 times larger than in free air. To explain this enhanced diffusion, Yosida et al. (1955) introduced the “hand-to-hand” delivery mechanism, which notably considers that the deposition of water molecules on one side of an ice grain and the sublimation on another side act as a shortcut in the vapor trajectory. Several subsequent experimental studies have either confirmed (e.g., Sommerfeld et al., 1987) or contradicted (Sokratov and Maeno, 2000) the findings of Yosida et al. (1955) that macroscopic vapor diffusion is significantly larger in snowpacks than in free air. Similarly, several analytical and numerical modeling works have either accepted (Colbeck, 1993; Christon et al., 1994; Gavriliev, 2008; Hansen and Foslien, 2015) or contradicted (Giddings and LaChapelle, 1962; Calonne et al., 2014; Shertzer and Adams, 2018) the results of Yosida et al. (1955) and the hand-to-hand mechanism. As mentioned by Sokratov and Maeno (2000) and Pinzer et al. (2012), the experimental discrepancies can be explained by the difficulty to accurately measure macroscopic vapor fluxes and vapor concentration gradients in snow, either in the field or in the laboratory. Yet, the large disagreement between the various analytical and modeling works, which sometimes differ more than 10-fold (e.g., Colbeck, 1993; Calonne et al., 2014), cannot be explained by experimental errors.

The aim of this paper is to clarify the origin of these discrepancies and to quantify the macroscopic vapor flux based on theoretical and numerical modeling. As the kinetics of sublimation and deposition of water molecules on the ice surfaces in snow is not well constrained, we decided to explore a broad range of possible kinetics in our study. We start by considering in Sect. 2 whether the hand-to-hand mechanism, as originally proposed by Yosida et al. (1955), can indeed explain the large macroscopic vapor fluxes observed in snow. Then in Sect. 3, we recall how the macroscopic vapor flux can be obtained from the microscopic vapor flux occurring at the pore scale. In Sect. 4 we present theoretical work to bound the macroscopic vapor flux in snow by treating two limiting cases of surface kinetics. Finally, numerical simulations are presented in Sect. 5 in order to illustrate the points raised throughout the article and to provide some numerical values of the effective diffusion coefficient.

As previously mentioned, the experiment of Yosida et al. (1955) marks the introduction of the idea of enhanced vapor diffusion due to the hand-to-hand delivery mechanism. Their experimental setup consisted of four stacked cans (3.5 cm in height and 5.5 cm in diameter each) filled with snow and separated with wire meshes that held the snow in place in each can without preventing vapor diffusion between them. A temperature difference was imposed between the top and bottom of the stack in order to create a vertical thermal gradient of about 45 K m^{−1} and thus induce a macroscopic vapor flux. The experiments were carried out with average temperatures of about −4 ^{∘}C and lasted about 5 h. The cans filled with snow were weighed before and after the experiment in order to determine their mass gain or loss, which can be used to estimate the magnitude of the macroscopic vapor flux transporting mass from one can to another. Based on these measurements and assuming that vapor was at saturation concentration, Yosida et al. (1955) concluded that the macroscopic vapor flux was about 3 to 4 times greater than what would be expected in free air for a similar concentration gradient. Noting that this result appears to contradict the idea that the presence of ice would impede the diffusion of vapor in snow, Yosida et al. (1955) proposed the hand-to-hand delivery mechanism as an explanation for this contradiction. This mechanism first states that because of its low thermal conductivity, the pore space of snow tends to concentrate the thermal gradient, leading to a concentrated vapor gradient in the pores. Moreover, Yosida et al. (1955) proposed that “Water vapor needs not force its way through the interspaces between the ice grains composing snow. It needs only condense on one side of an ice grain and evaporate from the other side to condense again on the side facing to it of the next grain. In this way the distance which the water vapor actually traverses by diffusion turns out to be a fraction of the distance of its displacement. Such a situation makes the diffusion of water vapor through snow easier than through open air, which causes D [the effective diffusion coefficient in snow] to appear greater than D_{0} [the diffusion coefficient in free air]”. One should note that this explanation entails more than the simple continuous sublimation of vapor from some interfaces and subsequent deposition on others. Yosida et al. (1955) argued that this is equivalent to a situation in which a molecule depositing on one side of an ice grain reappears as a sublimating molecule on another side.

Our understanding however is that the second part of the mechanism proposed by Yosida et al. (1955) is not physically sound and that the continuous deposition and sublimation of molecules cannot be used to explain their experimental results. A schematic illustration of the experiment is given in Fig. 1, with only two cans for simplicity. The hand-to-hand delivery of water molecules is represented by the orange and red dots, depositing on the lower side and sublimating on the upper side of the ice grain at the interface between the two cans. For this mechanism to explain the experimental observations, the continuous deposition and sublimation should produce a real mass flux from one can to the other, as if the depositing molecule reappeared as the sublimating one. However, what actually happens is that the depositing molecule (represented as an orange dot in Fig. 1) remains incorporated at the bottom of the ice grain, thus remaining in the first can. Similarly, the sublimating molecule (represented as a red dot in Fig. 1) was already present in the second can. The synchronous sublimation and deposition therefore do not lead to a mass transfer between the two cans. This is different from the molecules traversing the boundary in the air space (represented as green dots in Fig. 1) that actually lead to a mass transfer by depleting the first can in favor of the second one. We therefore argue that the hand-to-hand mechanism, as proposed by Yosida et al. (1955), is not physically sound.

If one adopts the hand-to-hand mechanism, such as Hansen (2019) for instance, the idea of water vapor shortcutting the ice may appear supported by the indistinguishability of water molecules. For an observer focused on the pore space, the argument says it really appears as if the water vapor is transported almost instantaneously through the ice as a disappearing water molecule depositing on one side of an ice grain is almost instantaneously replaced by an appearing molecule sublimating on the other side. However, this point of view neglects the fact that the mass leaving a control volume also depends on the gain or loss of the ice during the deposition and sublimation process. As exemplified in the right panel of Fig. 1, for an observer focused on the ice everything appears as if the ice disappearing on the sublimation side reappeared on the depositing side (see for instance the videos in the supplements of Pinzer et al., 2012; Hagenmuller et al., 2019). Because of mass conservation during the sublimation and deposition process, the apparent flux of vapor skipping the ice is compensated by an equal counter-flux of water molecules in the ice space. Therefore, the mass transfer from one control volume to another is solely governed by the diffusion of water molecules in the air (green dots in Fig. 1).

We stress that we do not disagree with the insightful propositions of Yosida et al. (1955) (i) that the vapor flux tends to travel from one ice grain to another and not to go around them and (ii) that the thermal gradient is enhanced in the pore space compared to the macroscopic gradient. The point of contention is that the continuous sublimation and deposition of water molecules do not count as a contribution to the mass flux. This problem with the hand-to-hand mechanism has been previously addressed by Giddings and LaChapelle (1962), when they noted that “The hand-to-hand transfer does not contribute to the flux because this transfer does not shift water molecules across a plane fixed in the solid network”.

The problem at hand is now to quantify the impact of the enhanced thermal gradient in the air space on the macroscopic diffusion of vapor and to determine whether it can account for the large macroscopic vapor fluxes reported in the literature (e.g., Yosida et al., 1955; Sommerfeld et al., 1987) and in particular if the macroscopic diffusion fluxes in snow are larger than the fluxes in free air.

Let us consider a volume of snow (Fig. 2a) subjected to vertical macroscopic temperature and vapor gradients at its boundaries. For this study we consider that the macroscopic water vapor gradient equals the macroscopic gradient of saturated vapor and is therefore driven by the macroscopic temperature gradient (as in Yosida et al., 1955; Colbeck, 1993; Sokratov and Maeno, 2000; Pinzer et al., 2012). A necessary condition to be able to treat this snow sample as an equivalent macroscopic medium is the condition of separation of scales (Auriault, 1991; Auriault et al., 2010). This separation of scale can be expressed as

where *L*_{micro} is the length scale characterizing the size of the representative elementary volume (REV) (Auriault et al., 2010; Calonne et al., 2014) of the microstructure, and *L*_{macro} is the length scale characterizing variations in the snowpack or of the external forcing applied at the macroscopic scale, for instance the change between different snow layers or changes in thermal and vapor gradients (Fig. 2b). In this study we consider snow samples with a size of at least *L*_{micro} but less than *L*_{macro}. In this case, the snow sample is large enough to be treated as an equivalent macroscopic body but not so large that it spans several snow layers and can thus be considered to be macroscopically homogeneous. The relation between the various length scales is exemplified in Fig. 2.

At the microscopic scale, vapor diffuses in response to vapor concentration gradients in the pore space. The resulting microscopic vapor fluxes ** f** are governed by Fick's law: $\mathit{f}=-{D}_{\mathrm{0}}\mathbf{\nabla}c$, with

*D*

_{0}being the diffusion coefficient of vapor in air and ∇

*c*the gradient of vapor concentration in the pore. These microscopic fluxes may result in a net transport of mass at the macroscopic scale, i.e., a macroscopic flux. The magnitude of this macroscopic flux

**corresponds to the mass transported through an orthogonal plane per unit time and per unit surface of snow. This macroscopic flux is the quantity that Yosida et al. (1955) set out to measure.**

*F*This paper, as previous works in the scientific literature, determines the macroscopic flux from the first principles of physics at the pore scale. It is therefore necessary to determine how the macroscopic flux ** F** at the macroscopic scale can be obtained from the microscopic fluxes

**in the pores. One might attempt to compute**

*f***as the quantity of matter transported through an arbitrary plane of the microstructure. In this case,**

*F***would be given as the surface average of the pore-scale flux**

*F***, with the averaging performed over the entire plane, ice included (the vapor flux being 0 in the ice). Yet, this method of computing the macroscopic vapor flux can be problematic. As pointed out by Pinzer et al. (2012) the water vapor fluxes through different horizontal planes of a microstructure are not necessarily all equal. Thus, depending on the plane chosen, the same snow sample could be assigned different macroscopic fluxes contrary to the notion that the snow sample is homogeneous from the macroscopic point of view. To avoid this issue, the macroscopic flux should therefore be computed as the volume-averaged microscopic vapor flux over the entire representative volume of the microstructure (Shertzer and Adams, 2018), which is equivalent to averaging the fluxes through various horizontal planes (Pinzer et al., 2012). Again, the averaging needs to be performed over the total volume, including the ice space, and the macroscopic vapor flux**

*f***is thus given by**

*F*
where *V* and *V*_{a}, respectively, represent the total volume of the snow sample and the pore volume.

We now phenomenologically define the effective diffusion coefficient for vapor *D*_{eff} such that $\mathit{F}=-{D}_{\mathrm{eff}}\mathrm{\nabla}C$, where ∇*C* is the macroscopic vapor concentration gradient (Colbeck, 1993; Shertzer and Adams, 2018). Here, the vapor concentration is expressed in mass per volume of pore space, and the averaging is thus performed in the pore only. The macroscopic vapor gradient is thus given by the difference in average vapor concentration between two opposing sides of the snow sample divided by the size of the sample. This corresponds to the definition implicitly adopted by Yosida et al. (1955). In the snow science community the effective diffusion coefficient *D*_{eff} is usually expected to be independent of the applied thermal and vapor gradients (e.g., Yosida et al., 1955; Colbeck, 1993). In this case, it is possible to treat the problem of macroscopic vapor transport in snow with a generalized Fick's law, where *D*_{eff} is independent of the applied boundary conditions and only depends on the snow microstructure. Such an effective diffusion coefficient does not depend on the external conditions and is then said to be intrinsic (Auriault et al., 2010). However, one should keep in mind that the effective diffusion coefficients computed in this work might depend on the applied vapor and thermal gradients and are therefore not necessarily intrinsic. Moreover the proposed numerical values may also not apply in the case where the macroscopic concentration gradient is decoupled from the macroscopic thermal gradient. Finally, we define the normalized effective diffusion coefficient as ${D}_{\mathrm{eff}}^{\mathrm{norm}}={D}_{\mathrm{eff}}/{D}_{\mathrm{0}}$. Expressing macroscopic water vapor fluxes in snow under the form of normalized diffusion coefficients allows us to easily compare them to free air.

Note that the goal of this work is only to quantify the macroscopic water vapor flux in snow and its associated phenomenological effective diffusion coefficient. Contrary to Calonne et al. (2014) we do not attempt to derive the macroscopic equations governing water vapor at the layer scale.

Let us consider a snow sample of volume *V* subjected to vertical thermal and vapor concentration gradients. For simplicity, we assume the problem to be steady-state. The diffusion of water vapor at the microscopic scale is governed by the following system of equations (Calonne et al., 2014):

where Ω_{a}, Γ, and ** n** represent the pore space, the ice–pore interface, and the normal vector to Γ pointing toward the ice.

*D*

_{0}is the vapor diffusion coefficient in free air,

*c*is the vapor concentration in the pores,

*c*

_{sat}is the vapor saturation concentration at the ice interface, ${v}_{\mathrm{kin}}=\sqrt{\left(kT\right)/\left(\mathrm{2}\mathit{\pi}m\right)}$ is related to the velocity of water molecules in the gas and is referred to as the kinetic velocity (

*k*being Boltzmann's constant and

*m*the mass of a water molecule), and

*α*is the sticking coefficient of water molecules on the ice surface (sometimes referred to as the accommodation coefficient) and is less than or equal to unity. The second equation of the system is the Hertz–Knudsen equation and governs the mass fluxes that are incorporated or released from the ice. In the presence of a large-enough thermal gradient, the dependence of the saturation concentration on the local curvature of the ice surface can be neglected compared to its dependence on temperature (Colbeck, 1983). Under this condition, we can expect

*c*

_{sat}to become a function of temperature only. Moreover, even if curvature effects were not negligible at the microscopic level, it appears unlikely for them to result in a net macroscopic vapor flux as in a homogeneous snow layer curvature differences are distributed isotropically within the microstructure and thus do not result in a net movement of water vapor.

The actual value of the *α* coefficient is not well known and in general will depend on the local saturation of water vapor and on the crystallographic properties of the ice surface (Saito, 1996; Libbrecht and Rickerby, 2013). Yet, two limiting cases, corresponding to the case of infinitely fast surface kinetics and inert ice surfaces, can easily be analyzed. As is empirically verified later, these two cases appear to correspond to the upper and lower bounds of macroscopic vapor fluxes in snow. Solving Eq. (3) we obtain the microscopic vapor fluxes inside the whole microstructure. Using Eq. (2) yields the water vapor flux at the macroscopic scale ** F**.

## 4.1 The infinitely fast surface kinetics case

In the case where the product *α**v*_{kin} is very large, small oversaturations (undersaturations) lead to an abrupt adsorption (desorption) of water molecules, rapidly restoring the saturation value. In the limiting and hypothetical case of infinitely fast surface kinetics (i.e., *α**v*_{kin}→∞), the vapor concentration is constantly at saturation at the ice–pore interface, and the Hertz–Knudsen equation can be replaced by the simpler equality of the vapor concentration with its saturation value at the ice surface. While the infinitely fast kinetics case is strictly theoretical as *α**v*_{kin} is less than or equal to *v*_{kin}, it helps apprehending the macroscopic vapor flux when surface kinetics processes are much faster than diffusion in the air space. Note also that the saturation of water vapor at the interface does not mean that the deposition and sublimation fluxes are 0 at the interface.

As explained by Pinzer et al. (2012), the infinitely fast surface kinetics situation is the case where the microscopic vapor gradients across the pores are maximal and therefore where the macroscopic vapor flux is also maximal. A demonstration of this fact, using the spatial-averaging theorem, is given in Appendix A. Note that the assumption of saturated vapor at the ice surface, and therefore infinitely fast surface kinetics, has been regularly employed in studies about the diffusion of vapor in snow (e.g., Colbeck, 1993; Christon et al., 1994; Pinzer et al., 2012).

Even though this case corresponds to the maximal vapor flux, it can be shown that the macroscopic diffusion coefficient remains less than expected in free air, as pointed out by Giddings and LaChapelle (1962). This is due to the loss of diffusion space because of the ice, and we propose here to rederive the Giddings and LaChapelle (1962) demonstration using a more detailed framework. First, we assume that the thermal gradient is low enough so that the saturation vapor concentration dependence on temperature can be considered to be linear. For a thermal gradient of 100 K m^{−1} applied to a 1 cm sample, the deviation of vapor concentration from linear behavior is about 0.1 %, while the deviation of the derivative with respect to temperature is about 5 %. Moreover, this condition corresponds to the fact that the macroscopic vapor gradient should be constant over the sample, i.e., that the size of the sample is smaller than *L*_{macro}.

Under this assumption one can show that the vapor concentration is at saturation within the entire pore space. A demonstration is presented in Appendix B, and a similar conclusion was also reached by Yosida et al. (1955) and Pinzer et al. (2012). Consequently, the macroscopic vapor flux is expressed as

where *ϕ* is the snow porosity (not to be confused with the ice volume fraction), *V*_{a} is the volume of the pore space, ∇*T*_{a} is the microscopic temperature gradient in the air; we have used the chain rule $\mathrm{\nabla}{c}_{\mathrm{sat}}=\frac{\mathrm{d}{c}_{\mathrm{sat}}}{\mathrm{d}T}\mathrm{\nabla}{T}_{\mathrm{a}}$. As we considered that the saturation concentration of vapor does not deviate from a linear behavior, $\frac{\mathrm{d}{c}_{\mathrm{sat}}}{\mathrm{d}T}$ is taken as constant over the volume *V*_{a}. Thus

The precise relationship between the average microscopic thermal gradient in the air space and the macroscopic gradient ∇*T* depends on the particular snow microstructure (Calonne et al., 2011, 2014; Hansen and Foslien, 2015). However, Hansen and Foslien (2015) report that

where *V*_{i} is the volume of the ice space, and ∇*T*_{i} is the microscopic temperature gradient in the ice.

As snow is a transversely isotropic material with the vertical direction being the direction normal to the isotropy plane, one can expect for reason of symmetry that the average air and ice thermal gradients are aligned with the vertical macroscopic gradient. Moreover, the average air and ice thermal gradients are oriented in the same direction as the macroscopic gradient. Therefore, one has the inequality about the magnitudes of the air and macroscopic thermal gradients

which states that while the average thermal gradient in the air can be greater than the macroscopic thermal gradient, it cannot exceed it by a factor greater than 1∕*ϕ*. Intuitively, it states that the temperature drop in the pore space cannot exceed the temperature drop observed over the entire snow sample. One can show that the air thermal gradient is maximal in the special case of a microstructure composed of slabs perpendicular to the macroscopic temperature gradient. In this case the temperature gradient is almost entirely concentrated in the air, and furthermore Eq. (7) becomes an equality when the thermal conductivity of ice is assumed to be infinite.

Using the inequality of Eq. (7) in Eq. (5) leads to an inequality of the magnitude of the macroscopic flux

where $\mathrm{\nabla}C=\frac{\mathrm{d}{c}_{\mathrm{sat}}}{\mathrm{d}T}\mathrm{\nabla}T$ is the macroscopic vapor concentration gradient.

The macroscopic vapor flux is thus less than the vapor flux that would take place in free air, which can be similarly expressed by ${D}_{\mathrm{eff}}^{\mathrm{norm}}\le \mathrm{1}$. While the microscopic vapor flux in the pores is enhanced due to the enhancement of the microscopic temperature and vapor gradients, this effect is countered by the reduction in the space where vapor can diffuse. As the average air temperature gradient at the maximum is enhanced by a factor 1∕*ϕ*, while the reduction in pore space systematically decreases the macroscopic flux by a factor *ϕ*, the resulting macroscopic vapor flux cannot be greater than in free air. The equality ${D}_{\mathrm{eff}}^{\mathrm{norm}}=\mathrm{1}$ holds when the entire temperature gradient is concentrated in the pore space. However, since the thermal conductivity of ice is finite, the thermal gradient cannot be solely concentrated in the pore space, and thus one always has ${D}_{\mathrm{eff}}^{\mathrm{norm}}<\mathrm{1}$.

## 4.2 The slow-surface-kinetics case

The other limiting case is when the deposition and sublimation of water vapor at the ice grain surfaces is slow enough to be neglected. The diffusion of water vapor in snow then becomes equivalent to the diffusion of a gas in an inert porous structure. This problem has been extensively studied (e.g., Torquato and Haslach, 2002; Auriault et al., 2010), and in this case the effective diffusion coefficient is given by

where *τ* is defined as the tortuosity factor and is linked to the lengthening of the diffusion streamlines in the porous network. The tortuosity factor represents an impediment of diffusion and is thus less than or equal to unity. Moreover, *τ* depends solely on the structure of the porous medium and not on the specific diffusive species or the applied concentration gradient (Torquato and Haslach, 2002; Auriault et al., 2010). Under an assumption of slow surface kinetics, Calonne et al. (2014) report effective diffusion coefficients reduced from 20 % to 85 % compared to the free-air case, with lower diffusion coefficients corresponding to denser snow samples. Although we do not have a rigorous demonstration of this fact, it appears that the slow-kinetics assumption corresponds to the case where the macroscopic flux (and hence *D*_{eff}) is minimal for a given vapor concentration gradient. This proposition is empirically verified with numerical simulations in Sect. 5.

## 4.3 Comparison with previous works

We establish in Sect. 4.1 that even under the assumption of fast surface kinetics, the effective vapor diffusion coefficient in snow cannot be greater than that in free air. Yet several studies based on analytical and numerical models, which are not subjected to experimental errors, have reported opposite results. It thus appears important to elucidate why those previous results do not invalidate the demonstration made in Sect. 4.1 and the results of this work.

Colbeck (1993) proposed a theoretical model based on an idealized structure of disconnected and equally spaced ice spheres. In that model the vapor concentration is at saturation at the ice surface (i.e., surface kinetics is infinitely fast), and the vapor flux between two consecutive spheres can be analytically computed. In this case, the author concludes that the vapor diffusion coefficient is between 4 and 7 times greater than in air. However, as pointed out by Pinzer et al. (2012), Colbeck (1993) derives the diffusion coefficient in snow by computing the flux crossing a single plane between two spheres and not by averaging over the entire volume. As the plane between two spheres corresponds to a zone of maximal thermal gradient without any ice blockage, it is not surprising that the local microscopic vapor flux is severalfold that in free air. However, as is seen in Sect. 5.1, computing the macroscopic flux by performing a volume averaging of microscopic vapor fluxes over the entire microstructure significantly reduces the corresponding diffusion coefficient, down to a value below that of free air.

Christon et al. (1994) performed finite-element microscale simulations of vapor diffusion in snow under a thermal gradient using an idealized microstructure. They concluded that the vapor diffusion coefficient is between 1 and 2 times as large as that in air. Yet, in that study the macroscopic mass flux is not computed as a volume average but rather “as the weighted average of the mass flux rates over all of the exterior surfaces of the diffusion domain in order to capture the bulk vertical mass diffusion rate”. Here, the diffusion domain refers to the domain where vapor diffusion occurs, i.e., the pore space. This differs from volume averaging and leads to an overestimation of the macroscopic flux as the ice space is not included. As the loss of diffusion space due to the ice is neglected, the effective diffusion coefficient is overestimated by a factor of 1∕*ϕ*.

Similarly, Pinzer et al. (2012) performed finite-element microscale simulations of vapor diffusion, this time with microstructures measured by X-ray computed microtomography scanning. A diffusion coefficient slightly greater than that in free air is reported. Pinzer et al. (2012) noted that computing the mass flux crossing a single plane was insufficient for the reasons discussed in Sect. 3. To derive the macroscopic mass flux, Pinzer et al. (2012) computed the average mass flux in each plane and then averaged over all planes. However, it appears from the description of their methodology that the slice averaging was only performed in the pore space, not taking into account the reduction in macroscopic flux due to the presence of ice. As in the case of Christon et al. (1994), this would explain the diffusion coefficient higher than in free air. As is shown in Sect. 5, performing similar numerical simulations and computing the macroscale flux by total volume averaging leads to diffusion coefficients below that in free air.

Finally, Hansen and Foslien (2015) and Hansen (2019) proposed an analytical expression for the effective thermal conductivity of snow, taking into account the latent heat associated with the transport of water vapor. In their model, water vapor is at constant saturation in the pores (thus corresponding to the case of infinitely fast kinetics) and acts as an integral part of heat transfer by transporting latent heat between sublimation and deposition surfaces (as notably proposed by Yosida et al., 1955). One application of this effective-thermal-conductivity model is to allow the derivation of the vapor flux, which leads to the conclusion that the macroscopic vapor flux is greater than that in free air. To come to this conclusion, Hansen and Foslien (2015) determine the vapor flux by identifying the contribution of latent heat in their expression of the effective thermal conductivity. However, during the identification of the latent-heat contribution to the total energy flux, some of the heat conduction contribution of the ice is attributed to the latent-heat transport. This leads to an artificially increased vapor flux, and therefore an overestimated diffusion coefficient. A rederivation of the vapor flux with the thermal conductivity expression proposed by Hansen and Foslien (2015) is presented in Appendix C and leads to a macroscopic vapor flux below that in free air.

Most of the discrepancies between our results and those of the published literature thus reduce down to computations of the macroscopic fluxes that are inconsistent with fluxes expressed per unit surface of snow, as used in snow models and experimental studies. This leads to an overestimation of the value of the effective diffusion coefficient. Focusing on the magnitude of microscopic vapor fluxes as done by Colbeck (1993) or Christon et al. (1994) is of great interest for snow metamorphism as they govern the mass transfer between adjacent ice grains and the recrystallization rate. However, they do not correspond to the macroscopic mass flux expressed per unit surface of snow, as measured by Yosida et al. (1955) and subsequent experimental studies (e.g., Sokratov and Maeno, 2000). We reiterate that the macroscopic vapor flux responsible for the redistribution of mass at the macroscopic scale and which inspired the hand-to-hand delivery mechanism corresponds to the volume-averaged flux over the entire snow microstructure and must include the loss of diffusion space due to the ice.

In this section we present steady-state 3D numerical simulations of vapor diffusion in snow subjected to a macroscopic temperature gradient ∇*T* and a macroscopic vapor gradient ∇*C*. The macroscopic temperature gradient ∇*T* is obtained by imposing the top and bottom temperatures *T*^{top} and *T*^{bot}. The vapor concentrations in the pore space at the top and bottom of the sample are imposed to correspond to the saturation values for the top and bottom temperatures. We thus have

where *L*_{z} is the height of the sample considered. Conditions of zero heat and vapor normal fluxes are imposed on the other sides of the sample. For simplicity, we only consider the case of vertical temperature and vapor gradients, although the extension to the other directions is straightforward. Moreover, we do not take into account the impact of latent heat on the temperature field. At the microscopic level, adding latent-heat effects would act as an additional mechanism transporting heat from the warm sublimating surfaces towards the cold deposition surfaces. It would cool the sublimation surfaces and warm the deposition surfaces, decreasing the thermal gradient in the pore space. Therefore, taking latent-heat effects into account would not increase the effective vapor diffusion coefficient.

The thermal conductivities of the ice and the air *k*_{i} and *k*_{a} are set to 2.34 and 0.024 $\mathrm{W}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}$, respectively (Riche and Schneebeli, 2013), and the diffusion coefficient of vapor in air *D*_{0} is set to $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ m^{2} s^{−1} (Calonne et al., 2014). The vapor concentration is assumed to follow the Clausius–Clapeyron and ideal-gas laws, leading to

where $M=\mathrm{18}\times {\mathrm{10}}^{-\mathrm{3}}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{mol}}^{-\mathrm{1}}$ is the molar mass of water, $R=\mathrm{8.314}\phantom{\rule{0.125em}{0ex}}\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{mol}}^{-\mathrm{1}}$ is the ideal-gas constant, $\mathrm{\Delta}{H}_{\mathrm{s}}=\mathrm{51}\times {\mathrm{10}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{mol}}^{-\mathrm{1}}$ is the latent heat of sublimation of ice, *T*_{0}=273.15 K is a reference temperature, and *P*_{0}=611 Pa is the saturation pressure of vapor over ice at *T*_{0}. The different physical constants used in this article are tabulated in Appendix D with their references. All simulations are performed with an average temperature $({T}^{\mathrm{bot}}+{T}^{\mathrm{top}})/\mathrm{2}=\mathrm{258}\phantom{\rule{0.125em}{0ex}}\mathrm{K}$.

The heat and diffusion equations are solved using the finite-element method with the open-source software ElmerFEM (Malinen and Råback, 2013). We use the readily available ElmerFEM modules dedicated to the heat and diffusion equations, which are solved with iterative methods. We first solve the steady-state heat equation in order to obtain the temperature field in the entire microstructure. The steady-state vapor diffusion equation is then solved using the saturation concentration at the ice–pore interface resulting from the previously computed temperature field. In the case of simulations performed on measured snow microstructures, the tetrahedral meshes have been derived from X-ray computed microtomography images using the CGAL meshing library. The meshes have been refined to capture the ice–pore interface and contain between 18 and 50 million elements, depending on the snow sample. Moreover, in the case of snow samples the meshes have been partitioned into 20 sub-meshes, and the computations are performed using the parallel computing abilities of ElmerFEM. Under such conditions, a simulation typically takes a bit less than an hour to run. Finally, the outputs of the simulations are processed using the ParaView software to compute the volume averages.

As seen previously, the kinetics of the sublimation and deposition processes at the ice surface might significantly impact the magnitude of the macroscopic vapor flux. We recall that in general the boundary condition at the ice–air interface is given by the Hertz–Knudsen equation:

where ${v}_{\mathrm{kin}}\simeq \mathrm{140}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ at 258 K, and *α* is the sticking coefficient less than or equal to unity. In general *α* is not a constant and depends on the local vapor saturation as well as the crystallographic properties of the underlying ice crystal (Saito, 1996; Libbrecht and Rickerby, 2013).

For each microstructure, several simulations were performed with different values of *α* in order to assess the impact of the internal boundary conditions (IBCs) applied at the ice surface. We first performed simulations with constant *α* equal to 0, 10^{−5}, 10^{−4}, 10^{−3}, 10^{−2}, 10^{−1}, and 1. Simulations with constant *α* are referred to as linear-kinetics simulations in what follows. Among them, a special case is *α*=0, which corresponds to the diffusion of vapor in an inert porous medium. Moreover, we performed simulations similar to those of Christon et al. (1994) and Pinzer et al. (2012), where the Hertz–Knudsen boundary condition is replaced with the saturation of vapor at the ice surface, corresponding to the infinitely fast kinetics case. Finally, we performed simulations in which the dependence of *α* on the local vapor saturation is explicitly represented. For that we set $\mathit{\alpha}=\mathrm{exp}(-{\mathit{\sigma}}_{\mathrm{0}}/\mathit{\sigma})$, where $\mathit{\sigma}=(c-{c}_{\mathrm{sat}})/{c}_{\mathrm{sat}}$, and *σ*_{0}=0.01. Note that this expression was determined for the attachment of vapor to the basal and prismatic facets of ice crystals (Libbrecht and Rickerby, 2013) and might not properly apply for the entirety of ice surfaces in snowpacks. Indeed, this law has been derived using deposition measurement and might not apply for sublimating surfaces (Beckmann and Lacmann, 1982). Moreover, the presence of vicinal surfaces in snowpacks, where the proposed law does not apply, is likely (Legagneux and Domine, 2005). Therefore, the point of using such a law is to qualitatively study the potential impact of a dependence of *α* on the local vapor saturation rather than to produce quantitative results. Simulations using this law are referred to as non-linear-kinetics simulations. Finally, the macroscopic fluxes of the various simulations are computed by performing a total volume average, as defined in Sect. 3, and the effective diffusion coefficients are obtained by dividing these macroscopic fluxes by the macroscopic concentration gradients, i.e., ${D}_{\mathrm{eff}}=-\mathit{F}/\mathrm{\nabla}C$.

## 5.1 Idealized structure

We start with an idealized microstructure composed of disconnected ice spheres, similar to that used by Colbeck (1993). The structure is visible in Fig. 3. The domain is a cuboid with dimensions of $\mathrm{3.7}\times \mathrm{3.7}\times \mathrm{10}\phantom{\rule{0.125em}{0ex}}{\mathrm{mm}}^{\mathrm{3}}$ and three equidistant ice spheres with 3 mm diameters, which are vertically aligned at the center of the domain. The distance between two sphere centers is set to 3.3 mm. This microstructure is characterized by a porosity of 0.619 and a density of 349 kg m^{−3}.

The simulations were performed for the different IBCs described previously and for temperature gradients ranging from 5 to 200 K m^{−1}. The resulting normalized effective diffusion coefficients are displayed in Fig. 4.

We first analyze the 50 K m^{−1} temperature gradient simulations. Illustrations of the microscopic vapor fluxes for three IBCs, namely inert surfaces (*α*=0), $\mathit{\alpha}={\mathrm{10}}^{-\mathrm{4}}$, and infinitely fast surface kinetics, are displayed in Fig. 3. In the case of inert surfaces the vapor flux needs to go around the ice grains, which act as a blockage, leading to tortuous streamlines. In the case of infinitely fast surface kinetics, the vapor flux does not need to go around the ice grain and is rather moving from ice grain to ice grain, in agreement with the suggestion of Yosida et al. (1955) and the numerical simulations of Pinzer et al. (2012). Finally, the $\mathit{\alpha}={\mathrm{10}}^{-\mathrm{4}}$ case displays an intermediate behavior, with some of the vapor flux moving from ice grain to ice grain, while the rest bypasses the ice. This exemplifies that the microscopic vapor fluxes are strongly dependent on the kinetics of the vapor sublimation and deposition at the ice surface.

In the case of infinitely fast surface kinetics we find a normalized diffusion coefficient of 0.978, i.e., lower than in air, in agreement with the calculations of Sect. 4.1. Moreover, we computed the average air temperature gradient (in the pore space only) and found it to be 79.00 K m^{−1}. This is enhanced compared to the 50 K m^{−1} macroscopic gradient but still respects the inequality of Eq. (7). While the enhancement of the thermal gradient increases the microscopic vapor fluxes in the pores, it does not suffice to counter the loss of diffusion space, and the resulting macroscopic flux is lower than in free air.

To compare our results to the works of Colbeck (1993), Christon et al. (1994), and Pinzer et al. (2012), who worked under a similar assumption of infinitely fast kinetics, we used two alternate methods, different from total volume averaging, to compute the vapor flux. The first consists of averaging the microscopic vapor fluxes in the air space only, and we call the associated normalized diffusion coefficient ${D}_{\mathrm{air}}^{\mathrm{norm}}$. The second one consists of computing the flux crossing a horizontal plane placed between two spheres, and we call the associated diffusion coefficient ${D}_{\mathrm{plane}}^{\mathrm{norm}}$. As explained in Sect. 4.3, we believe that the first methodology is akin to works of Christon et al. (1994) and Pinzer et al. (2012), while the second was used by Colbeck (1993). Calculations yield a ${D}_{\mathrm{air}}^{\mathrm{norm}}$ of 1.580 and a ${D}_{\mathrm{plane}}^{\mathrm{norm}}$ of 2.986, consistent with the values reported by Christon et al. (1994), Pinzer et al. (2012), and Colbeck (1993). By not including the ice in the averaging or by selecting a peculiar plane where microscopic vapor fluxes are maximum, the macroscopic vapor flux is overestimated, leading to a diffusion coefficient greater than *D*_{0}.

The outcome of the other simulations performed with $\mathrm{\nabla}T=\mathrm{50}\phantom{\rule{0.125em}{0ex}}\mathrm{K}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}$ is reported in Fig. 4 and indicates that ${D}_{\mathrm{eff}}^{\mathrm{norm}}$ is maximal in the infinitely fast kinetics case, with a value of 0.978, and minimal in the inert-surface case, with a value of 0.512. Accordingly, the normalized effective diffusion coefficient increases with *α* and for the cases *α*=0.1 and *α*=1 differs by less than 0.3 % from the infinitely fast case. The use of the non-linear-surface-kinetics law leads to a normalized effective diffusion coefficient equal to 0.857, in between the inert (${D}_{\mathrm{eff}}^{\mathrm{norm}}=\mathrm{0.512}$) and infinitely fast kinetics (${D}_{\mathrm{eff}}^{\mathrm{norm}}=\mathrm{0.978}$) cases.

Similar observations can be made for the simulations performed with other temperature gradients. For the entire range of gradients tested, the infinitely fast kinetics and inert-surface cases correspond to the maximal and minimal macroscopic fluxes. Moreover, the associated effective diffusion coefficients are mostly independent of the macroscopic thermal or vapor gradients, suggesting that the effective diffusion coefficients could be intrinsic in these cases. Consistent results are observed for the simulations where *α* is constant. The obtained effective diffusion coefficients are mostly independent of the applied macroscopic gradient and are bounded by the infinitely fast kinetics and inert-surface cases. Note that the *α*=0.1 and *α*=1 cases are indistinguishable from the infinitely fast kinetics results in Fig. 4. Contrary to the rest of the simulations, the non-linear IBC yields effective diffusion coefficients that depend on the magnitude of the applied gradients. In this case, the macroscopic vapor flux and the macroscopic vapor concentration gradient are not proportionally linked by a single and well-defined material property. Furthermore, for low vapor and thermal gradients the non-linear case is close to the inert-surface case, while a transition towards the fast-kinetics case is observed for thermal gradients around 50 K m^{−1}. Again, even though the non-linear law used to express *α* as a function of local saturation does not necessarily accurately model water molecule attachment in real snowpacks, it illustrates the effects of a non-constant *α*.

## 5.2 Measured snow microstructures

Other numerical simulations of vapor diffusion have been performed, this time using measured snow microstructures instead of the idealized structure of Sect. 5.1. The microstructures were obtained by X-ray computed microtomography imaging of snow samples. In total six snow samples were analyzed, covering the snow types of decomposing and fragmented precipitation particles (DF), depth hoar (DH), rounded grains (RG), and melt forms (MF) (Fierz et al., 2009). The goal is not to provide effective diffusion coefficients on an exhaustive set of snow microstructural patterns but to illustrate the effects of the snow microstructure and surface kinetics on water vapor diffusion.

A close-up view showing the vapor streamlines inside the melt form sample is provided in Fig. 5. As with the idealized microstructure, in the inert-surface case vapor tends to go around the ice grains. In the infinitely fast kinetics case, vapor moves from ice grain to ice grain, as proposed by Yosida et al. (1955) and reported by Pinzer et al. (2012).

We start by analyzing the results of the simulations of the DF sample, characterized by a density of 125 kg m^{−3}. Similarly to Sect. 5.1, the simulations were performed by imposing external temperature and vapor gradients, with different selected IBCs characterizing the kinetics of the vapor sublimation and deposition process. The results are displayed in Fig. 6. As in the idealized case, the inert-surface, infinitely fast kinetics, and linear-kinetics cases yield normalized effective diffusion coefficients that are mostly independent of the applied gradients. Moreover, we observe that ${D}_{\mathrm{eff}}^{\mathrm{norm}}$ is minimal in the inert-surface case, with a value of 0.764, and maximal in the infinitely fast kinetics case, with a value of 0.982. As expected, the effective diffusion coefficient is systematically lower than that of air. The normalized effective diffusion coefficients in the linear-kinetics cases are distributed between the inert and infinitely fast values and increase with the value of *α*. For *α*=1, ${D}_{\mathrm{eff}}^{\mathrm{norm}}$ differs by less than 0.1 *%* from the infinitely fast kinetics case.

In contrast, the non-linear-kinetics case leads to a normalized effective diffusion coefficient that depends on the external gradients. As with the idealized disconnected-sphere structure of Sect. 5.1, we observe that for low gradients the non-linear case is close to the slow-kinetics simulations and transition towards faster kinetics with higher gradients. However, in the case of the DF sample this transition occurs more slowly and with higher temperature and vapor gradients.

Since the normalized effective diffusion coefficients appear to be independent of the external thermal and vapor gradient in the case of infinitely fast and linear surface kinetics, we only computed ${D}_{\mathrm{eff}}^{\mathrm{norm}}$ with a 50 K m^{−1} gradient for the five remaining snow samples. We also did not compute ${D}_{\mathrm{eff}}^{\mathrm{norm}}$ with non-linear surface kinetics (i.e., when alpha is not constant) as we are not confident in the validity of the chosen non-linear law for snow modeling. The resulting ${D}_{\mathrm{eff}}^{\mathrm{norm}}$ values are reported in Table 1 and displayed in Fig. 7 as a function of the sticking coefficient *α*. Again, ${D}_{\mathrm{eff}}^{\mathrm{norm}}$ is systematically minimal in the inert-surface case and maximal in the infinitely fast kinetics case. Figure 7 highlights that the normalized effective diffusion coefficient exhibits two different regimes depending on the value of *α*. The transition between the fast- and slow-surface-kinetics regimes occurs for values of *α* around 10^{−3}.

We observe that the effective diffusion coefficient is well correlated with density and show an almost systematic decrease in ${D}_{\mathrm{eff}}^{\mathrm{norm}}$ with increasing density for all values of *α*. The correlation between ${D}_{\mathrm{eff}}^{\mathrm{norm}}$ and the specific surface area is not so well marked, notably for the RG sample that shows a large value of specific surface area without any clear impact on ${D}_{\mathrm{eff}}^{\mathrm{norm}}$. That being said, our sample set is only composed of six samples, for which density and specific surface area are correlated. A detailed study of the influence of microstructural parameters on the effective diffusion coefficient would require a larger sample set, notably to be able to decipher the independent influence of specific surface area and density.

We have shown that the macroscopic vapor flux in snow is less than the flux in free air under the same water vapor gradient. This result is supported by a formal demonstration, inspired by the work of Giddings and LaChapelle (1962) as well as by numerical simulations of idealized and measured snow microstructures. While the interaction of water vapor with the ice structure results in a macroscopic flux larger than that of the inert-diffusion case, the macroscopic vapor flux cannot be enhanced compared to the free-air case. We have shown that most of the previous theoretical studies reporting enhanced macroscopic vapor flux compared to free air used faulty computations of the macroscopic vapor flux, which resulted in systematic overestimation.

As seen in this work, the sublimation and deposition fluxes at the ice surface play a great role in the final macroscopic flux. In particular we have shown that when the reaction is fast, i.e., *α* is large, the macroscopic fluxes can be close to those that would be observed in free air. Moreover, the dependence of *α* on the local vapor saturation might break the proportionality between the macroscopic vapor gradient and the macroscopic flux. In this case, it is no longer possible to define a single effective diffusion coefficient *D*_{eff} that proportionally relates the vapor flux to vapor gradient and that solely depends on the snow microstructure. In other words, with non-linear surface kinetics *D*_{eff} is not intrinsic. For all these reasons, it appears important to determine what the precise internal boundary conditions are that govern the sublimation and deposition of water vapor in snowpacks and in particular to determine whether the inert surfaces or infinitely fast kinetics case could accurately describe real snow. In the case of fast kinetics, one can have *D*_{eff}≥*ϕ**D*_{0} as the average microscopic vapor gradient can be greater than the macroscopic vapor gradient. In contrast, in the case of slow surface kinetics one has ${D}_{\mathrm{eff}}=\mathit{\varphi}\mathit{\tau}{D}_{\mathrm{0}}\le \mathit{\varphi}{D}_{\mathrm{0}}$ since *τ*≤1. An experimental distinction between fast and slow kinetics could thus be made by observing whether the quantity *D*_{eff}∕(*ϕ**D*_{0}) is greater than unity or not. Using the experimental results of Sokratov and Maeno (2000), which are the experimental results with the lowest reported diffusion coefficient, we observe that *D*_{eff}∕(*ϕ**D*_{0}) is almost always greater than unity, which supports the notion of fast rather than slow kinetics. This is consistent with the study of Krol and Löwe (2016), which reports that fast kinetics is consistent with their microtomography-based observation of the temperature gradient metamorphism of a snow sample. That being said, experimental determination of the macroscopic vapor fluxes is difficult, as exemplified by the large spread of reported values, and more observations would be needed to make a decisive conclusion on this point.

This work investigated the effective diffusion coefficient of vapor in snow with a phenomenological approach, where the diffusion coefficient is simply defined as the ratio of the macroscopic vapor flux to the vapor concentration gradient. A rigorous upscaling of the microscale equations to derive the equivalent macroscopic formulation would greatly benefit the understanding and modeling of the macroscopic vapor flux. Note that such an approach was used by Calonne et al. (2014) with the method of asymptotic-scale expansion but limited itself to small *α*. Applying a similar method to the case of non-negligible surface sublimation and deposition would lead to a proper definition of the macroscopic quantities, notably of the effective diffusion coefficient, and to the proper formulation of the equations governing the macroscopic scale. Furthermore, we assumed in this study that the macroscopic water vapor gradient is equal to the macroscopic gradient of saturated vapor, driven by the macroscopic thermal gradient. This assumption has been regularly made in the snow science community (Yosida et al., 1955; Colbeck, 1993; Sokratov and Maeno, 2000; Pinzer et al., 2012) and is supported by the idea that the ice in the snowpack tends to impose water vapor saturation at the macroscopic scale. It however remains possible that the macroscopic water concentration deviates from saturation, notably if the deposition and sublimation kinetics is slow. A rigorous upscaling method yielding the equations governing macroscopic water concentration would therefore also help quantifying if such a situation of non-saturation at the macroscopic scale is likely to occur in real snowpacks and indicate how the macroscopic vapor flux should be computed in such a case.

Finally, the fact that there is no macroscopic enhancement of the water vapor flux in snow suggests that most of the mass flux observed in subarctic and Arctic snow and which would necessitate effective diffusion coefficients several times higher than that of free air to be explained solely by diffusion (e.g., Sturm and Benson, 1997; Domine et al., 2016, 2018) could rather be due to convection. The importance of convective mass transport in subarctic snowpacks has notably been pointed out by Trabant and Benson (1972) and Sturm and Johnson (1991) and thus appears to be a good candidate to explain the high vapor movement in subarctic snowpacks. Currently, detailed snow physics models do not include the mechanism of convective mass transport (Lehning et al., 2002; Vionnet et al., 2012) and assume that all mass transport results from diffusion, sometimes using a diffusion coefficient larger than that in free air (e.g., Jafari et al., 2020). Further modeling efforts to include convective mass transport in detailed snow models could enhance their ability to model snowpack evolution.

This work investigated the macroscopic vapor fluxes that arise in snowpacks due to large-scale vapor gradients. We first considered the seminal work of Yosida et al. (1955) and their formulation of the hand-to-hand delivery mechanism, which was meant to explain the large vapor flux they measured. We argue that it is reasonable to assume that the concentration of the thermal gradient in the pore space would lead to strong vapor gradients between ice grains and drive the sublimation of water molecules from some grains and subsequent deposition on others. Yet, we disagree with the proposed idea that the process where one water molecule deposits on one side of an ice grain while another molecule sublimates on the other side is equivalent to a situation where the depositing molecule skipped the ice, virtually increasing the vapor flux.

We demonstrated that the specific internal boundary conditions governing the sublimation and deposition of water molecules have a significant impact on the macroscopic vapor flux. In particular, we showed that in the case of infinitely fast kinetics the macroscopic flux is enhanced compared to the slow-kinetics case but still cannot exceed the vapor flux that would happen in free air under an equivalent vapor gradient. This demonstration is confirmed by numerical simulations of both idealized and measured snow microstructures. The discrepancies with previous studies that report vapor fluxes greater than the free-air case originate from erroneous computations of how the macroscopic flux was obtained from the microscopic vapor fluxes at the pore scale. We argue that the method used in this article, i.e., volume averaging over an entire microstructure including the ice, is the only one consistent with the actual nature of the macroscopic water vapor flux.

The numerical simulations also indicate that the infinitely fast kinetics and inert ice surface cases, respectively, are the upper and lower limits for the vapor flux in snow. The use of more complex laws describing the sublimation and deposition of water molecules at the ice surface leads to flux values in between both previously mentioned cases. Moreover, the use of a non-constant attachment coefficient breaks the proportionality between the macroscopic vapor flux and the vapor gradient. In that case, it is no longer possible to define a constant and intrinsic effective diffusion coefficient, proportionally relating the macroscopic vapor flux to the macroscopic concentration gradient, independently of the magnitude of applied vapor concentration gradient.

The aim of this appendix is to demonstrate that the macroscopic vapor flux is maximal in the case of infinitely fast kinetics. For this we start by applying the spatial-averaging theorem (Whitaker, 1999) to the vapor concentration in the pores *c*

where 〈•〉 is an operator defined as $\langle \u2022\rangle =\frac{\mathrm{1}}{V}{\int}_{{V}_{\mathrm{a}}}\u2022\mathrm{d}V$, and the concentration *c* in the surface integral is the vapor concentration at the ice–pore interface. Multiplying by *D*_{0} and using the notation introduced in this article for the macroscopic vapor flux ** F**, we have

Moreover, using the Hertz–Knudsen equation we find that the concentration at the interface is

Equation (A2) can thus be written as

Applying the same spatial-averaging theorem to the saturation concentration *c*_{sat}, we have

Injecting Eq. (A5) in Eq. (A4) thus yields

As we assume that the macroscopic vapor concentration gradient equals the macroscopic saturation concentration gradient (as in Yosida et al., 1955; Colbeck, 1993; Sokratov and Maeno, 2000; Pinzer et al., 2012), we have that ∇〈*c*〉=∇〈*c*_{sat}〉. Thus

Let us now assume, without loss of generality, that the macroscopic vapor and thermal gradients are oriented downward. As seen in Fig. A1, surfaces that are characterized by a normal vector pointing upward are deposition surfaces. The product ∇*c*⋅** n** is therefore negative, and (∇

*c*⋅

**)**

*n***is a vector pointing downward. Similarly, surfaces that are characterized by a normal vector pointing downward are sublimation surfaces. The product ∇**

*n**c*⋅

**is thus positive, and the vector (∇**

*n**c*⋅

**)**

*n***is pointing downward. Therefore, for both types of surfaces (∇**

*n**c*⋅

**)**

*n***is pointing downward. The surface integral term in Eq. (A7) thus acts in opposition to −〈**

*n**D*

_{0}∇

*c*

_{sat}〉 and tends to reduce the macroscopic vapor flux. We thus have the inequality

We now show that this upper bound is reached in the infinitely fast kinetics case. Indeed, under the infinitely fast kinetics hypothesis the product *α**v*_{kin} can be treated as going to infinity. At the same time, the surface integral of Eq. (A7) remains bounded as the concentration gradient in the vicinity of the interface does not diverge. The surface integral thus vanishes, and the norm of the vapor flux is given by

that is to say that the upper bound of the macroscopic vapor flux is reached under infinitely fast kinetics. Moreover, note that we rederived that, in the infinitely fast kinetics case, the macroscopic vapor flux is given by the spatial average of the saturation vapor concentration in the pore space.

In the case of infinitely fast surface kinetics and assuming a linear relation between saturation concentration and temperature, the equations governing the vapor concentration are

where *A* and *B* are two constants characterizing the linear relationship between temperature and vapor concentration, and *T* is the temperature of the ice surface. Thanks to the linearity of the divergence and gradient operators and owing to the fact that ∇*B*=0, the equations can be reformulated to

where $\mathit{\theta}=(c-B)/A$, and we have used the fact that *D*_{0} is a non-zero constant to eliminate it from the first equation. Moreover let us recall that in the air temperature *T*_{a} is a solution of the following Laplace equation:

Systems of Eqs. (B2) and (B3) are identical, and since the solution of such a boundary value problem is unique, it follows that ${T}_{\mathrm{a}}=\mathit{\theta}=(c-B)/A$ over the entire pore space. It thus follows that $c=A{T}_{\mathrm{a}}+B={c}_{\mathrm{sat}}\left({T}_{\mathrm{a}}\right)$ in the pores.

Hansen and Foslien (2015) proposed that the heat flux *q*_{s} through a snow sample under a macroscopic thermal gradient ∇*T* be expressed as

where *q*_{tub} and *q*_{lam} are the heat fluxes through idealized snow structures corresponding, respectively, to a tubular structure and a lamellae structure, submitted to the same macroscopic thermal gradient ∇*T*, and *ϕ* is the porosity of the snow sample (not to be confused with the ice volume fraction). Concerning the tubular microstructure, one has

where *k*_{i} and *k*_{a} are the thermal conductivities of ice and air, and *L* is the latent heat of sublimation of ice. The contribution of the vapor flux is $\mathit{\varphi}L{D}_{\mathrm{0}}\frac{\mathrm{d}{c}_{\mathrm{sat}}}{\mathrm{d}T}\mathrm{\nabla}T$, and the vapor flux in the tubular microstructure is $\mathit{\varphi}{D}_{\mathrm{0}}\frac{\mathrm{d}{c}_{\mathrm{sat}}}{\mathrm{d}T}\mathrm{\nabla}T=\mathit{\varphi}{D}_{\mathrm{0}}\mathrm{\nabla}C$.

Similarly, concerning the lamellae microstructure one has

In their article, Hansen and Foslien (2015) then rewrite *q*_{lam} under the form

and identify the second term with the latent-heat flux. We however argue that Eq. (C4) is only one way among many to rewrite *q*_{lam} under the form $A\mathrm{\nabla}T+L\frac{\mathrm{d}{c}_{\mathrm{sat}}}{\mathrm{d}T}B\mathrm{\nabla}T$ and thus that the identification of the latent-heat flux with the second term of the decomposition is arbitrary.

To derive the latent-heat flux, we first start from Eqs. (79) and (80) of Hansen and Foslien (2015) and compute the thermal gradients in the ice and in the air, respectively, as

and

The heat flux *q*^{cond} through the sole process of conduction is thus given by

meaning that the latent-heat flux, which is the remaining contribution to *q*_{lam}, is ${q}_{\mathrm{lam}}-{q}^{\mathrm{cond}}=\frac{\mathit{\varphi}{k}_{\mathrm{i}}L{D}_{\mathrm{0}}\frac{\mathrm{d}{c}_{\mathrm{sat}}}{\mathrm{d}T}}{(\mathrm{1}-\mathit{\varphi})({k}_{\mathrm{a}}+L{D}_{\mathrm{0}}\frac{\mathrm{d}{c}_{\mathrm{sat}}}{\mathrm{d}T})+\mathit{\varphi}{k}_{\mathrm{i}}}\mathrm{\nabla}T$ and that the vapor flux is $\frac{\mathit{\varphi}{k}_{\mathrm{i}}{D}_{\mathrm{0}}}{(\mathrm{1}-\mathit{\varphi})({k}_{\mathrm{a}}+L{D}_{\mathrm{0}}\frac{\mathrm{d}{c}_{\mathrm{sat}}}{\mathrm{d}T})+\mathit{\varphi}{k}_{\mathrm{i}}}\mathrm{\nabla}C$. Note that the *ϕ* term in the numerator is not present in the original Hansen and Foslien (2015) demonstration, leading to an overestimation of the vapor flux.

Finally, the total vapor flux in the Hansen and Foslien (2015) model is computed as the weighted average of the tubular and lamellae vapor fluxes

and the expression in square bracket is therefore the effective vapor diffusion coefficient that one can show to be less than *D*_{0}.

The physical constants used in this article are listed in Table D1 with their units, numerical values, and references.

Calonne et al. (2014)Lide (2006)Lide (2006)Riche and Schneebeli (2013)Riche and Schneebeli (2013)The codes used for the simulations were developed with Python 3 and ElmerFEM. They will be provided upon request to the corresponding author.

FD designed the research with input from PH and KF. FD obtained funding. KF performed research and wrote the paper with input from FD and PH.

The authors declare that they have no conflict of interest.

This work contributes to the APT project (Acceleration of Permafrost Thaw), funded by the Climate Initiative program of the BNP Paribas Foundation. We thank Jacques Roulle for his help during the tomography imaging and cold-room work. Neige Calonne and Marie Dumont provided valuable input on the subject. Henning Löwe kindly indicated the usage of the spatial-averaging theorem. We thank Kevin Hammonds and Quirine Krol for their helpful review of the manuscript and Jürg Schweizer for editing the article.

This research has been supported by the BNP Paribas Foundation (Acceleration of Permafrost Thaw program).

This paper was edited by Jürg Schweizer and reviewed by Kevin Hammonds and Quirine Krol.

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

Auriault, J.-L., Boutin, C., and Geindreau, C.: Homogenization of coupled phenomena in heterogenous media, John Wiley & Sons, Hoboken, New Jersey, USA, 2010. a, b, c, d, e

Beckmann, W. and Lacmann, R.: Interface kinetics of the growth and evaporation of ice single crystals from the vapour phase: II. Measurements in a pure water vapour environment, J. Cryst. Growth, 58, 433–442, https://doi.org/10.1016/0022-0248(82)90292-5, 1982. a

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

Calonne, N., Geindreau, C., and Flin, F.: Macroscopic modeling for heat and water vapor transfer in dry snow by homogenization, J. Phys. Chem. B, 118, 13393–13403, https://doi.org/10.1021/jp5052535, 2014. a, b, c, d, e, f, g, h, i, j, k

Christon, M., Burns, P. J., and Sommerfeld, R. A.: Quasi-steady temperature gradient metamorphism in idealized, dry snow, Numer. Heat Transf., Part A, 25, 259–278, https://doi.org/10.1080/10407789408955948, 1994. a, b, c, d, e, f, g, h, i

Colbeck, S. C.: Theory of metamorphism of dry snow, J. Geophys. Res.-Oceans, 88, 5475–5482, https://doi.org/10.1029/JC088iC09p05475, 1983. a, b

Colbeck, S. C.: The vapor diffusion coefficient for snow, Water Resour. Res., 29, 109–115, https://doi.org/10.1029/92WR02301, 1993. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p

Domine, F., Barrere, M., and Sarrazin, D.: Seasonal evolution of the effective thermal conductivity of the snow and the soil in high Arctic herb tundra at Bylot Island, Canada, The Cryosphere, 10, 2573–2588, https://doi.org/10.5194/tc-10-2573-2016, 2016. a

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

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 Classificationi for Seasonal Snow on the Ground, IHP-VII Technical Documents in Hydrology No. 83, IACS Contribution No. 1, UNESCO-IHP, Paris, 2009. a, b

Gavriliev, R.: A Model for Calculating the Effective Diffusion Coefficient of Water Vapour in Snow, in: Ninth International Conference on Permafrost, edited by: Kane, D. L. and Hinkel, K. M., 505–501, Institute of Northern Engineering, University of Alaska Fairbanks, Fairbanks, Alaska, USA, 2008. a

Giddings, J. C. and LaChapelle, E.: The formation rate of depth hoar, J. Geophys. Res., 67, 2377–2383, https://doi.org/10.1029/JZ067i006p02377, 1962. a, b, c, d, e, f

Hagenmuller, P., Flin, F., Dumont, M., Tuzet, F., Peinke, I., Lapalus, P., Dufour, A., Roulle, J., Pézard, L., Voisin, D., Ando, E., Rolland du Roscoat, S., and Charrier, P.: Motion of dust particles in dry snow under temperature gradient metamorphism, The Cryosphere, 13, 2345–2359, https://doi.org/10.5194/tc-13-2345-2019, 2019. a

Hansen, A.: Revisiting the vapor diffusion coefficient in dry snow, The Cryosphere Discuss. [preprint], https://doi.org/10.5194/tc-2019-143, 2019. a, b, c

Hansen, A. C. and Foslien, W. E.: A macroscale mixture theory analysis of deposition and sublimation rates during heat and mass transfer in dry snow, The Cryosphere, 9, 1857–1878, https://doi.org/10.5194/tc-9-1857-2015, 2015. a, b, c, d, e, f, g, h, i, j, k, l

Jafari, M., Gouttevin, I., Couttet, M., Wever, N., Michel, A., Sharma, V., Rossmann, L., Maass, N., Nicolaus, M., and Lehning, M.: The Impact of Diffusive Water Vapor Transport on Snow Profiles in Deep and Shallow Snow Covers and on Sea Ice, Front. Earth Sci., 8, 249, https://doi.org/10.3389/feart.2020.00249, 2020. a

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

Legagneux, L. and Domine, F.: A mean field model of the decrease of the specific surface area of dry snow during isothermal metamorphism, J. Geophys. Res.-Earth Surf., 110, F04011, https://doi.org/10.1029/2004JF000181, 2005. a, b

Lehning, M., Bartelt, P., Brown, B., Fierz, C., and Satyawali, P.: A physical SNOWPACK model for the Swiss avalanche warning: Part II. Snow microstructure, Cold Reg. Sci. Tech., 35, 147–167, https://doi.org/10.1016/S0165-232X(02)00073-3, 2002. a

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.jcrysgro.2013.04.037, 2013. a, b, c

Lide, D. R.: CRC handbook of chemistry and physics, chap. Properties of ice and supercooled water, 6–5, CRC press, Taylor and Francis, Boca Raton, FL, 85th edn., 2006. a, b

Malinen, M. and Råback, P.: Elmer Finite Element Solver for Multiphysics and Multiscale Problems, in: Multiscale Modelling Methods for Applications in Materials Science, edited by: Kondov, I. and Sutmann, G., 101–113, Forschungszentrum Jülich GmbH, Jülich, Germany, 2013. a

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. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v

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. a, b, c

Saito, Y.: Statistical physics of crystal growth, World Scientific, Singapore, 1996. a, b

Shertzer, R. H. and Adams, E. E.: A Mass Diffusion Model for Dry Snow Utilizing a Fabric Tensor to Characterize Anisotropy, J Adv. Model. Earth Sys., 10, 881–890, https://doi.org/10.1002/2017MS001046, 2018. a, b, c, d

Sokratov, S. A. and Maeno, N.: Effective water vapor diffusion coefficient of snow under a temperature gradient, Water Resour. Res., 36, 1269–1276, https://doi.org/10.1029/2000WR900014, 2000. a, b, c, d, e, f, g

Sommerfeld, R. A., Friedman, I., and Nilles, M.: The fractionation of natural isotopes during temperature gradient metamorphism of snow, in: Seasonal Snowcovers: Physics, Chemistry, Hydrology, 95–105, D. Reidel Publishing Company, Dordrecht, the Netherlands, 1987. a, b

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

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

Torquato, S. and Haslach Jr., H.: Random heterogeneous materials: microstructure and macroscopic properties, Appl. Mech. Rev., 55, B62–B63, 2002. a, b

Trabant, D. and Benson, C.: Field experiments on the development of depth hoar, Geol. Soc. Am. Mem., 135, 309–322, 1972. a

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

Whitaker, S.: The method of volume averaging, vol. 13, Springer, the Netherlands, https://doi.org/10.1007/978-94-017-3389-2, 1999. a

Yosida, Z., Oura, H., Kuroiwa, D., Huzioka, T., Kojima, K., Aoki, S.-I., and Kinosita, S.: Physical Studies on Deposited Snow. I. Thermal Properties, Contributions from the Institute of Low Temperature Science, 7, 19–74, 1955. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac

- Abstract
- Introduction
- Does the hand-to-hand mechanism enhance macroscopic vapor diffusion?
- Defining the macroscopic vapor flux and the effective diffusion coefficient
- Bounding the effective diffusion coefficient of water vapor in snow
- Numerical modeling
- Discussion
- Conclusions
- Appendix A: Demonstration that the macroscopic vapor flux is maximal under infinitely fast kinetics
- Appendix B: Saturation of vapor in the infinitely fast surface kinetics case
- Appendix C: Vapor flux in the Hansen and Folsien (2015) thermal conductivity
- Appendix D: Physical constants
- Code availability
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Does the hand-to-hand mechanism enhance macroscopic vapor diffusion?
- Defining the macroscopic vapor flux and the effective diffusion coefficient
- Bounding the effective diffusion coefficient of water vapor in snow
- Numerical modeling
- Discussion
- Conclusions
- Appendix A: Demonstration that the macroscopic vapor flux is maximal under infinitely fast kinetics
- Appendix B: Saturation of vapor in the infinitely fast surface kinetics case
- Appendix C: Vapor flux in the Hansen and Folsien (2015) thermal conductivity
- Appendix D: Physical constants
- Code availability
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References