3-D image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy

We used three-dimensional (3-D) images of snow microstructure to carry out numerical estimations of the full tensor of the intrinsic permeability of snow (K). This study was performed on 35 snow samples, spanning a wide range of seasonal snow types. For several snow samples, a significant anisotropy of permeability was detected and is consistent with that observed for the effective thermal conductivity obtained from the same samples. The anisotropy coefficient, defined as the ratio of the vertical over the horizontal components of K, ranges from 0.74 for a sample of decomposing precipitation particles collected in the field to 1.66 for a depth hoar specimen. Because the permeability is related to a characteristic length, we introduced a dimensionless tensor K∗=K/r2es, where the equivalent sphere radius of ice grains (res) is computed from the specific surface area of snow (SSA) and the ice density (ρi) as follows: res=3/(SSA×ρi). We define K and K∗ as the average of the diagonal components of K and K ∗, respectively. The 35 values of K∗ were fitted to snow density (ρs) and provide the following regression: K=(3.0±0.3) r2esexp((−0.0130±0.0003)ρs). We noted that the anisotropy of permeability does not affect significantly the proposed equation. This regression curve was applied to several independent datasets from the literature and compared to other existing regression curves or analytical models. The results show that it is probably the best currently available simple relationship linking the average value of permeability, K, to snow density and specific surface area.


Introduction
The intrinsic permeability (K) is an important physical property of snow. Defined in a tensorial way, K (m 2 ) links the pressure gradient ∇p (Pa m −1 ) and the discharge per unit area q (m s −1 ) through the Darcy's law q = −(1/µ)K∇p, where µ is the dynamic viscosity of the fluid (kg m −1 s −1 ). The intrinsic permeability is of particular interest to quantify transport properties of snow and firn, such as wind pumping processes (Colbeck, 1989(Colbeck, , 1997, air convection (Akitaya, 1974;Powers et al., 1985;Brun and Touvier, 1987;Sturm and Johnson, 1991;Albert et al., 2004) and liquid water flow (Colbeck, 1975(Colbeck, , 1976Waldner et al., 2004;Katsushima et al., 2009;Yamaguchi et al., 2010). It is a key variable for a wide number of applications such as atmospheric and firn chemistry (Freitag et al., 2002;Neumann, 2003;Grannas et al., 2007;Domine et al., 2008) or snow-firn metamorphism (Albert, 2002;Hörhold et al., 2009). Permeability has also been proposed as a means of characterization for quantitative snow classification (Bader, 1939;Jordan et al., 1999;Arakawa et al., 2009). Indeed, most physical properties of snow, including its permeability, are linked to the geometrical arrangement of ice, air, water vapor and sometimes liquid water, referred to as its microstructure.
The first quantitative investigations of snow permeability are attributed to Bader (1939). Shimizu (1970), Sommerfeld and Rocchio (1993), Jordan et al. (1999) and Arakawa et al. (2010) realized extensive experiments on seasonal snow and proposed parameterizations depending on grain size and density. However, these studies are characterized by a significant scatter in the results, which seems mainly due to experimental issues (Sommerfeld and Rocchio, 1993) and to sometimes ambiguous definitions of grain size (Jordan et al., 1999). In addition existing field permeameters require snow samples to be several cm large (cylinder diameter on the order of 10 to 20 cm typically, and cylinder height from several to tens of cm) (Shimizu, 1970;Albert et al., 2000;Arakawa et al., 2009) which generally exceeds the size within which snow properties can reasonably be considered homogeneous (Schneebeli et al., 1999;Pielmeier and Schneebeli, 2003;Schneebeli, 2006, 2010;Marshall and Johnson, 2009).
More recently, the availability of 3-D images of snow and firn from X-ray tomography (Brzoska et al., 1999;Coléou et al., 2001;Freitag et al., 2004;Schneebeli and Sokratov, 2004;Kaempfer et al., 2005;Chen and Baker, 2010) opened the way to numerical simulations, applying to snow the methods developed for porous media in general (Spanne et al., 1994;Ferreol and Rothman, 1995;Martys and Chen, 1996;Arns et al., 2001Arns et al., , 2004. In particular, such techniques allow to estimate the permeability for snow samples much smaller than those used for experimental measurements, thereby reducing if not eliminating issues associated with the intrinsic heterogeneity of the samples probed. Recently, Freitag et al. (2002) and Courville et al. (2010) computed one component of the permeability of firn samples using Lattice-Boltzman modeling and Zermatten et al. (2011) performed direct porelevel simulations on five snow samples. While these numerical computations methods are particularly promising, the obtained results show only a qualitative agreement with previous experimental studies and theoretical models, and do not take into account the anisotropy of the permeability.
We carried out numerical estimations of the full 3-D tensor of intrinsic permeability (K) on 35 3-D images of seasonal snow obtained from microtomography. Computations were performed with the software Geodict (Thoemen et al., 2008;Koivu et al., 2009;Calonne et al., 2011) and were based on the periodic homogenization method (Auriault et al., 2009). The main objective of this paper was to elaborate a parameterization of the snow permeability from other variables measurable in the field. For this purpose, we studied the relationship between the computed permeability, snow density and grain size, defined here as the equivalent sphere radius, at the scale of the representative elementary volume (REV). This relationship obtained using our data was compared to existing literature datasets as well as other equations from theoretical models and regressions. In addition, we focused on the anisotropy of permeability, available from the computed 3-D tensor of permeability.

Snow samples
Numerical computations were performed on 35 tomographic images obtained from previous experiments or field sampling (e.g. Calonne et al., 2011), spanning most types of seasonal snow, i.e. precipitation particles (PP), decomposing and fragmented precipitation particles (DF), rounded grains (RG), faceted crystals (FC), depth hoar (DH) and melt forms (MF), according to the International Classification for Seasonal Snow on the Ground (ICSSG) (Fierz et al., 2009).
Two-thirds of the snow samples come from controlled cold-room experiments: a first series (PP, DF and RG) was obtained by subjecting deposited natural snow to isothermal conditions at 271 K (Flin et al., 2004). Another RG sample was obtained under similar conditions, but after sieving. A second series (RG, FC and DH) was obtained under a temperature gradient of 43 K m −1 at 269 K and corresponds to various stages of metamorphism of the initial sieved snow (RG) into FC then DH. Two other similar experiments, with a gradient of 16 and 100 K m −1 at 268 and 270 K, provided two samples (FC and DH) (Coléou et al., 2001;Flin and Brzoska, 2008). A series of MF samples was obtained by grain coarsening in water-saturated snow using the method of Raymond and Tusima (1979) followed by draining their liquid water content (Coléou et al., 2001;Flin et al., 2011).
The remaining snow samples were directly collected in the field: 10 snow specimens (PP, DF and RG) were sampled at increasing depths in the snowpack of the Girose glacier (Écrins, French Alps) . Another PP sample was collected at Col de Porte (Chartreuse, French Alps).
All 3-D images considered in this study are cubic, have an edge size ranging from 2.51 to 9.66 mm, and a resolution between 4.91 micrometer (µm) and 10 µm, depending on the snow type. Additional information are indicated in Table 1 of the Supplement. Auriault (1991Auriault ( , 2011 has shown that physical phenomena in random and in periodic medium can be modeled by a similar equivalent continuous macroscopic description, provided that the condition of separation of scales is satisfied. This fundamental condition may be expressed as ε ≪ 1, where ε represents the separation of scale and is defined as ε = l/L, in which 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 of the material and the physical phenomena. The REV constitutes the smallest fraction of the sample volume from which a variable representative of the whole can be determined. In practice, periodic boundary conditions are thus widely used to compute effective properties of random media (Kanit et Bernard et al., 2005;Koivu et al., 2009;Calonne et al., 2011). Assuming for snow that ε ≪ 1, we applied the homogenization method for periodic structures to derive the Darcy's law from the physics at the pore scale (Ene and Sanchez-Palencia, 1975), as described in the following. Let us consider a rigid porous matrix, which is periodic with period , i.e. the size of the REV, and is fully saturated by an incompressible Newtonian fluid of density ρ and viscosity µ. s and f are the domains occupied by the solid and the fluid, respectively. The common boundary of s and f is denoted Ŵ. The porosity is defined as φ = | f |/| |. At the pore scale, the steady state flow is described by

Computations of snow permeability
where v and p are the velocity and the pressure of the fluid, respectively. Using the homogenization method, it can be shown that the corresponding macroscopic description strongly depends on the order of magnitude of the pore Reynolds number, R el = |ρ(v · ∇)v|/|µ△v| = O(ρv c l/µ) where v c is a characteristic value of fluid velocity. When R el < O(1), the macroscopic flow is described by Darcy's law, otherwise non-linearities appear at the macroscopic scale (Mei and Auriault, 1991). In what follows, it is assumed that R el < O(1), which is typically the case if we consider airflow through snow and firn induced by moderate winds (< 6 m s −1 ) over the snowpack surface (Albert, 2002). Under such conditions, the macroscopic description is written: where v = q represents the Darcy's velocity and K is the intrinsic permeability tensor of the porous media. This tensor, which is symmetric and positive, is defined as K ij = k ij . The second order tensor k is the solution of the following boundary value problem over the REV, where v = −(1/µ)k∇p, the pressure fluctuationp (with p = 0) are the periodic unknowns and ∇p is a given macroscopic gradient of pressure. The components of the permeability tensor K were estimated by solving the above boundary value problem (Eqs. 3a-3c) on a REV extracted from tomographic images (see Sect. 3.6), using the software Geodict (http://www.geodict.de). The boundary value problem is solved by using the finite difference method. Within this method a staggered grid (voxel) is used. The values of velocity and pressure are defined at center points of the faces and volume of the cubic unit cells, respectively.
The partial differential equations for incompressible Stokes flow (Eqs. 3a-3c) are solved by using the FFF-Stokes solver based on fast Fourier transform. Periodic boundary conditions are applied on the external boundaries of each volume (see Wiegmann, 2007 for more details).
Before carrying out the computations, voxels that are part of the network of interconnected pores were detected by image analysis, allowing to determine the ratio between closed and open porosity and to check that it was very small (less than 0.004 for all the samples). This means that it is correct to consider that air can flow through the whole porosity of the REVs for all of our samples.
In the following, the non-diagonal terms of the tensor K, about 50 times lower than diagonal terms, are not presented (the x-, y-and z-axes of 3-D images correspond to the principal directions of the microstructure, z being along the direction of the gravity). We note K x , K y and K z the diagonal term of the permeability tensor computed in the x-, yand z-direction, respectively, K the average value of the three terms, and K xy the average value of K x and K y . In the following, we mostly use K z , K xy and K, i.e. the vertical component, the average of the two horizontal components, and the average value of the three components of the permeability tensor, respectively.

Dimensionless permeability
The intrinsic permeability is strongly linked to a characteristic length of the microstructure of the medium considered (Boutin and Geindreau, 2010). Because the dimension of the permeability is a square length, K is often normalized by a characteristic length to the square, leading to a dimensionless tensor that we note K * .
Among various existing length metrics for snow (Fierz et al., 2009), we focus on the equivalent sphere radius of snow (r es , in m) (Sommerfeld and Rocchio, 1993;Luciano and Albert, 2002), also called the optical radius, opticalequivalent grain size or OGS (e.g. Grenfell and Warren, 1999;Painter et al., 2006;Fierz et al., 2009;Brucker et al., 2011). It is a characteristic length of the ice grains at the microscopic scale, which corresponds to the radius of a monodisperse collection of spheres having the same specific surface area (SSA) value than the sample considered. The snow SSA is defined as the total surface area of the air-ice interface per unit mass and can be quantitatively estimated by various means experimentally (e.g. Matzl and Schneebeli, 2006;Domine et al., 2008;Gallet et al., 2009;Arnaud et al., 2011) and numerically using 3-D images (e.g. Flin et al., 2011). The equivalent sphere radius and snow SSA are related by the following equation: where ρ i = 917 kg m −3 is the ice density. In the following, K * thus corresponds to K/r 2 es and we keep the notation K * , www.the-cryosphere.net/6/939/2012/ The Cryosphere, 6, 939-951, 2012 K * z and K * xy for the average, vertical and horizontal components of K * , respectively.
One could also choose to normalize K by a characteristic length corresponding to the pore space and use the hydraulic radius (r h ), commonly applied to flow through pipes and open channels. This normalization is equivalent to using r es , because these two radii are linked by the simple relationship: r h = r es (1 − φ)/(3φ) (Bear, 1972), where φ represents snow porosity. The granulometric analysis is another possible way to estimate a characteristic length (Zermatten et al., 2011). However, this approach is less convenient since it cannot be performed in the field, but only using 3-D images.
The specific surface area (SSA, in m 2 kg −1 ) was computed from 3-D images using a stereological method as described by Flin et al. (2011) where this quantity is obtained by averaging SSA estimations computed along 3 orthogonal directions (x, y, z). The equivalent sphere radius was then computed from SSA using Eq. (4).
The anisotropy coefficient of permeability, A(K), was computed from numerical estimations of the permeability tensor such as A(K) = K z /K xy .
The full tensor of the effective thermal conductivity (k eff , in W m −1 K −1 ) of the 35 snow samples considered was computed using the periodic homogenization. The effective thermal conductivity of most snow samples presented here were already described and included in a previous study (Calonne et al., 2011). Computations were extended to a few samples additionally considered in the present study. The coefficient of anisotropy of the effective thermal conductivity A(k eff ) was defined similarly to that of intrinsic permeability, by ratioing the vertical and horizontal components of the tensor.

Representative elementary volume
The representative elementary volume (REV) of our samples was estimated by calculating values of a given variable from several sub-volumes of increasing sizes within the same sample. The size of the REV was assumed to be reached once values did not vary significantly when the size of the sub-volumes of computation increased. Note that the REV size depends on the variable and on the sample studied: REVs with respect to permeability are generally equal to or larger than those for other variables, such as density, SSA or the effective thermal conductivity (Kanit et al., 2003;Rolland du Roscoat et al., 2007). Thus, a special attention must be paid to ensure that computations are carried out on a sufficiently large volume.

Overview of the numerical calculations
The average value of the three components of the intrinsic permeability of snow, noted K, ranges from 4 × 10 −10 to 6 × 10 −9 m 2 for the 35 samples considered in this study. This range of values is consistent with previous experimental and numerical estimates of snow permeability (e.g. Sommerfeld and Rocchio, 1993;Albert et al., 2000;Luciano and Albert, 2002;Domine et al., 2008;Arakawa et al., 2009;Courville et al., 2010;Zermatten et al., 2011). Density and SSA span from 103 to 544 kg m −3 and from 4 to 56 m 2 kg −1 , respectively. For each sample, detailed values are provided in Table 2 of the Supplement. Figure 1 provides an overview of the above results for the 35 samples of this study, showing the vertical (K * z ) and horizontal (K * xy ) components of the dimensionless permeability vs. snow density. PP samples exhibit the largest values of K * (1.05 for ρ s = 103 kg m −3 ), while the lowermost values are obtained for MF samples (2.16 × 10 −3 for ρ s = 544 kg m −3 ). The figure clearly shows that K * z , K * xy and thus K * decrease with increasing ρ s .

Anisotropy
As shown by Fig. 1, values of the vertical (K * z ) and horizontal (K * xy ) components of the dimensionless permeability are not identical, and some samples exhibit significant differences. Indeed, the anisotropy coefficients of permeability (A(K)) range from 0.74 for a DF sample collected in the field, to 1.66 for a particularly evolved DH sample obtained in cold room. We note that this range of A(K) values is consistent with the values between 0.75 and 1.9 measured by Luciano and Albert (2002) on firn collected from 1 to 13 m depth at Summit, Greenland.
For our snow specimen exhibiting a A(K) value of 1.66, Fig. 2 shows the air flow through the sample induced by a pressure drop along the vertical (left) and horizontal (right) direction. The flow, constrained by the snow microstructure, is clearly higher in the vertical than in the horizontal direction, leading to a high K z value. Note that the value of 1.66 is of the same order of magnitude as the analytical coefficient of anisotropy for a network of vertical cylinders (Boutin, 2000). Figure 3 shows the relationship between the anisotropy coefficient computed for the intrinsic permeability (A(K)) and for the effective thermal conductivity (A(k eff )) of snow. The good relationship between these two variables indicates that the microstructure influences both variables in a similar manner. Nevertheless, some discrepancies can be observed depending on density, snow type and microstructure. In addition, Fig. 3 shows that the anisotropy coefficients, and especially A(k eff ), enable the differentiation of the DH and FC samples (highest values) from other snow types. This result suggests that the anisotropy coefficients may be helpful for The Cryosphere, 6, 939-951, 2012 www.the-cryosphere.net/6/939/2012/ Fig. 1. Dimensionless permeability vs. snow density. "T" symbols indicate the values obtained by our numerical computations. Tips and horizontal bars of the "T" shapes represent the vertical (K * z ) and horizontal (K * xy ) components of K * , respectively. Colors correspond to the ICSSG (Fierz et al., 2009). Analytical models, numerical computations and fits are also plotted. quantitative classification of snow. However, due to the thin layered nature of the snowpack (e.g. Schneebeli et al., 1999;Schneebeli, 2006, 2010;Marshall and Johnson, 2009), macroscopic measurements seem particularly challenging to access such information except in special cases where the investigated material is sufficiently homogeneous (e.g. Luciano and Albert, 2002).

Regression analysis
The numerical data presented in Sect. 3.1 were used to build a regression curve allowing to infer the intrinsic permeability from snow density and equivalent sphere radius. Following the pioneering work of Shimizu (1970), a regression of this form was seeked: Other mathematical forms of equations were tested and none of them proved better than Eq. (5). For the 35 samples studied, the parameters a and b in Eq. (5) were calculated using the nonlinear least-squares Marquardt-Levenberg fitting algorithm, from the numerical estimates of permeability, r es and ρ s . Note that these computations were performed separately for each diagonal component of K, as well as for the average value K. The fitting algorithm used provides estimates of the uncertainty pertaining to the parameters a and b, in the form of "asymptotic standard errors". These constitute an indication of the fit's accuracy. Table 1 shows an overview of the fit parameters a and b computed using the values of K x , K y , K z and K, with the associated "asymptotic standard errors". For the whole computation, the "asymptotic standard errors" values are small (on the order of 10 % and 2 % for a and b, respectively), indicating a strong correlation between permeability, density and r es . Moreover, the table indicates that the parameters obtained from the four different values of K are insignificantly different from each other. The regression is thus not affected by the anisotropy of K presented in the previous section. Based on our numerical estimates of K, r es and ρ s , which span a wide range of snow types, we propose the following regression to infer an average permeability value K from r es and ρ s :  This regression corresponds to the black solid line in Fig. 1. Providing a regression which reflects the anisotropy of permeability seems currently challenging and further investigations are needed. Nevertheless, we believe that it would require an additional variable representing the snow microstructure.

Test of the obtained regression curve against literature data
This section compares our numerical data and the obtained regression curve (Eq. 6) with K, r es and ρ s data from the literature. Sommerfeld and Rocchio (1993), Arakawa et al.  (a and b) of Eq. (5) for the three components of K, as well as for their mean value K. "Asymptotic standard errors" (±) resulting from the fitting algorithm are given both in absolute units and as percentage values.  Figure 4 provides a general view of the values of dimensionless permeability vs. density from the four datasets described above (colored symbols), as well as from our 35 computed values and the associated regression curve (in black). Data from the four studies are overall consistent with our computations, showing a similar relationship of dimensionless permeability with ρ s , even if some discrepancies can be observed.
The predicted values of K provided by the regression curve (Eq. 6) were compared to observed permeability data. Figure 5 displays the mean and standard deviation of the relative residuals of the regression curve against our own permeability data as well as the above-mentioned datasets. It indicates that our regression curve manages to estimate the intrinsic permeability of snow with a small relative bias on the order of 20 % maximum within a standard deviation on the order of 40 % maximum. This result shows that the numerical computations of Courville et al. (2010), Zermatten et al. (2011) and this study are in good agreement, although all three studies use different numerical methods and boundary conditions on the external faces of the sample. The only major deviation from a good performance of Eq. (6) against experimental and numerical data is encountered with the dataset of Arakawa et al. (2009), which exhibits an overall positive relative bias of 71 % with a standard deviation of 85 %, also seen in Fig. 4 (Fierz et al., 2009), excepting stars which are used when firn is considered or when the snow type is not specified in the paper. Analytical models and fits are also plotted.
this discrepancy is (are) not understood. This could be due to the difficulty of making reliable and reproducible measurements of permeability. A first source of error may be the physical damage of the sample caused during its sampling or its handling (Sommerfeld and Rocchio, 1993). Shimizu (1970) and Sommerfeld and Rocchio (1993) also pointed out a bias linked to possible condensation/sublimation of the snow microstructure induced by the airflow imposed through the snow sample during the measurement. Sample heterogeneity could also be invoked. Based on the mathematical equation of our proposed regression curve (Eq. 6), accounting for a 10 % uncertainty on the parameter a of the equation (here we neglect the uncertainty on the parameter b of the equation -see Table 1), and assuming that r es and ρ s both carry a measurement uncertainty on the order of 10 % (Matzl and Schneebeli, 2006;Painter et al., 2006;Gallet et al., 2009;Conger and McClung, 2009;Arnaud et al., 2011), the propagation of these relative errors in terms of K adds up to about 50 %. This experimental error is of the same order of magnitude as the minimum and maximum deviation found when applying the regression curve to independent data (black points in Fig. 5).

Comparisons of the obtained regression curve to models and fits
This section compares the permeability estimates provided by our regression (Eq. 6) and by various classical equations proposed in the literature. We used the following regressions, analytical formulas or numerical computations (note that Figs. 1 and 4 show the corresponding curves in terms of dimensionless permeability: K * = K/r 2 es ): i. We recall the regression fit (Eq. 6) proposed in this study (referred to as "Calonne" in the relevant figures): ii. The well-known Shimizu's fit (Shimizu, 1970) is expressed as iii. The self-consistent (SC) estimate was obtained assuming that the porous medium consists of a bicomposite spherical pattern made of an internal spherical grain and an external fluid shell that ensures fluid connectivity. The porous medium is defined using the most basic information, i.e. the porosity and the grain size. Using the SC method, Boutin (2000) showed that the above porous medium leads to the following estimate: iv. In the Carman-Kozeny (CK) model, the medium is treated as a bundle of capillarity tubes of equal length. By solving the Stokes equations simultaneously for all the channels passing through a cross-section normal to the flow in the porous medium, the permeability is written as where SSA I is the specific surface area per unit of ice volume (SSA I = SSA×ρ i , in m −1 ) and c a coefficient which characterizes the geometry of the channels in the model. Empirically, c is found to be equal to 0.2 for many types of porous media (Bear, 1972). The CK equation can thus be expressed as v. Numerical values of permeability (K sphere ) were computed using the periodic homogenization method and a finite element software, for a simple cubic packing of spheres (Boutin and Geindreau, 2010). Note that at high density, where spheres interpenetrate, the SSA is computed from the effective surface formed by the sphere's assembly.
Figure 5 allows to investigate the performance of the above equations by comparing their estimates (K Cal , K Shi , K SC , K CK and K sphere ) against the independently observed permeability data presented in Sect. 3.4. For each expression, the result is given in terms of average and standard deviation of the relative residuals (in color on Fig. 5). The whole set of the relative residual values is also shown by the black points. Note that in the case of the simple cubic packing of spheres, the corresponding K values, normalized by r es and then expressed as a function of density, were interpolated to be able to compute K sphere estimates for all values of density and r es . Figure 5 shows that the Carman-Kozeny (CK) model, selfconsistent (SC) model, and the numerical computations for a simple cubic packing of spheres, predict the numerical values obtained in this study with a relative bias of the order of 20 % maximum within a factor 2 at most (maximum relative deviation of 40 %). These expressions are thus consistent with our proposed regression curve. In contrast, the regression fit proposed by Shimizu (1970) behaves poorly against our numerical values, showing a negative relative bias of −36 % within a relative deviation of 45 % (see Fig. 5). This situation is also true for other datasets, where the fit proposed by Shimizu (1970) (Fierz et al., 2009). For a given sample, each symbol represents a particular volume from which K, SSA and ρ s were computed. The largest volume corresponds to the total volume of the sample. Analytical models and fits are also plotted.
21 % (± 23 %) and −57 % (± 31 %)). As pointed out by several studies (Sommerfeld and Rocchio, 1993;Jordan et al., 1999), we believe that this is largely due to the method used by Shimizu (1970) to estimate grain size. From snow cross-sections, he computed a grain diameter D = 2 r es = √ 6ρ s /πρ i n, where n, the number of ice grains appearing in the cross section per a unit area, is estimated by a counting process described as "a complicated-shaped grain having m remarkable constrictions was counted as m + 1 grains". This visual method is subjective and may lead to erroneous estimates of n and therefore of K * .
Over the whole ensemble of tested datasets, Fig. 5 indicates that the best estimates are obtained using our regression (Eq. 6), showing very small average values of relative residuals, from +13 % for the numerical study of Courville et al. (2010) to −19 % for Zermatten et al. (2011). Moreover, the set of values of relative residuals generally does not exceed ± 50 %. Again, the predicted values for the dataset reported by Arakawa et al. (2009) do not correspond to this situation, as discussed in Sect. 3.4.
The predicted values using a simple packing of spheres are also consistent with experimental and numerical estimates from snow samples. In particular, this latter model behaves better than the CK and SC approach for the MF samples, where these two models fail to reflect the snow microstructure and overestimate the permeability, as shown in Fig. 1. These results can be explained by the fact that at low density the airflow around a snow particle is little affected by flow around its neighbors. In contrast, at higher density, snow particles are close together and the flow around one of them disturbs the flow around the others. This last phenomenon is not captured by analytical models.
Finally, Zermatten et al. (2011) indicate significant differences between the CK model and their values, while Fig. 5 shows an excellent agreement between both (the CK model predicts their observed data with a relative bias of 0.8 % within a relative deviation of 27 %): they apparently plotted an erroneous CK equation, using the specific surface area per unit of snow volume (SSA V , in m −1 ) in the expression of K = φ 3 /(5(1 − φ) 2 SSA 2 I ) instead of SSA I .

Representative elementary volume
Permeability estimations of the REV were performed on one sample of each snow type and their edge size range from 2.5 mm for the PP sample to 5.5 mm for the MF sample, which corresponds to volumes smaller or equal to those of our 3-D images (see Supplement for further details). For three snow samples (RG, FC and MF snow types), computations of K, SSA and ρ s were performed on the total available 3-D image and on sub-volumes of different sizes. Figure 6 shows the obtained values of K * as a function of ρ s . By increasing the volume size of calculation (from volume 1 to 4 or 5 in Fig. 6), the K * value is closer and closer to the result based on the previous volume, which is consistent with the definition of a REV. We can also see that if K, SSA and/or ρ s are computed on a volume smaller than their REV, the relationship between K * and ρ s remains consistent with the regression curve (Eq. 6) proposed in this study. This observation confirms the robustness of the relationship between K, SSA and ρ s .

Conclusions
The intrinsic permeability tensor K was computed on 35 tomographic images of various snow types by solving numerically a specific boundary value problem arising from the homogenization process on a representative elementary volume (REV). The equivalent sphere radius (r es ) was used as a characteristic length of the microstructure to reduce K to a dimensionless tensor K * . A regression equation using the 35 computed values of mean permeability (K), density (ρ s ) and equivalent sphere radius was proposed, such as: K = (3.0 ± 0.3) r 2 es exp ((−0.0130 ± 0.0003) ρ s ) and compared to existing literature data. Our main conclusions are summarized below: 1. The intrinsic permeability of snow can be anisotropic depending on the snow microstructure. The anisotropy coefficients of permeability range from 0.74 for a sample of decomposing precipitation particles collected in the field, to 1.66 for a particularly evolved depth hoar specimen, and are consistent with the anisotropy coefficient of the effective thermal conductivity. It appears that the use of these coefficients could be helpful for the quantitative classification of snow, as it may enable to distinguish the depth hoar and the facetted crystals from other snow types.
2. Permeability, density and equivalent sphere radius, directly related to the specific surface area, are strongly correlated. The equivalent sphere radius is thus a relevant characteristic length for permeability, which in addition is rather easy to determine in the field or from 3-D images. However, this strong correlation between snow permeability, density and SSA precludes consid-ering them as independent variables for the sake of objective snow classification.
3. The anisotropy of permeability does not affect the regression curve performed on the 35 snow samples. Indeed, very similar regression curves have been computed by using each diagonal component of K as well as the average of these three terms. Thus, the proposed regression allows estimating an average permeability value only.
4. Our numerical computations of permeability from 3-D images are consistent with datasets from other experimental and numerical studies. Moreover, our regression succeeds in estimating the permeability data of previous studies with a small relative bias on the order of 20 % maximum, within a factor ≈ 2 maximum relative deviation (40 %), except for the dataset of Arakawa et al. (2009) for which the agreement is lower (positive relative bias of 71 % with a standard deviation of 85 %). By comparing with other equations from theoretical models, fits and numerical data available in literature, our regression appears to be the best currently available simple relationship linking the average value of permeability to snow density and specific surface area. In particular, the well known fit proposed by Shimizu (1970) seems to significantly underestimate the permeability of snow for most of the tested datasets, most probably due to an inconsistency between the equivalent sphere radius derived from the specific surface area of snow and the empirically-defined snow grain size used by this author.