# Study of a temperature gradient metamorphism of snow from 3-D images: time evolution of microstructures, physical properties and their associated anisotropy

N. Calonne, F. Flin, C. Geindreau, B. Lesaffre, S. Rolland Du Roscoat

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

## To cite this version:

N. Calonne, F. Flin, C. Geindreau, B. Lesaffre, S. Rolland Du Roscoat. Study of a temperature gradient metamorphism of snow from 3-D images: time evolution of microstructures, physical properties and their associated anisotropy. The Cryosphere Discussions, 2014, 8 (1), pp.1407-1451. 10.5194/tc-8-2255-2014. hal-02021652

HAL Id: hal-02021652

https://hal.univ-grenoble-alpes.fr/hal-02021652

Submitted on 13 Dec 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

# Study of a temperature gradient metamorphism of snow from 3-D images: time evolution of microstructures, physical properties and their associated anisotropy

N. Calonne¹,²,³, F. Flin¹, C. Geindreau²,³, B. Lesaffre¹, and S. Rolland du Roscoat²,³

¹Météo-France – CNRS, CNRM – GAME UMR3589, CEN, 38400 Saint Martin d'Hères, France

²Univ. Grenoble Alpes, 3SR, 38000 Grenoble, France

³CNRS, 3SR, 38000 Grenoble, France

*Correspondence to:* N. Calonne (neige.calonne@meteo.fr) and F. Flin (frederic.flin@meteo.fr)

Received: 28 January 2014 – Published in The Cryosphere Discuss.: 28 February 2014

Revised: 3 September 2014 – Accepted: 26 September 2014 – Published: 5 December 2014

**Abstract.** We carried out a study to monitor the time evolution of microstructural and physical properties of snow during temperature gradient metamorphism: a snow slab was subjected to a constant temperature gradient in the vertical direction for 3 weeks in a cold room, and regularly sampled in order to obtain a series of three-dimensional (3-D) images using X-ray microtomography. A large set of properties was then computed from this series of 3-D images: density, specific surface area, correlation lengths, mean and Gaussian curvature distributions, air and ice tortuosities, effective thermal conductivity, and intrinsic permeability. Whenever possible, specific attention was paid to assess these properties along the vertical and horizontal directions, and an anisotropy coefficient defined as the ratio of the vertical over the horizontal values was deduced. The time evolution of these properties, as well as their anisotropy coefficients, was investigated, showing the development of a strong anisotropic behavior during the experiment. Most of the computed physical properties of snow were then compared with two analytical estimates (self-consistent estimates and dilute beds of spheroids) based on the snow density, and the size and anisotropy of the microstructure through the correlation lengths. These models, which require only basic microstructural information, offer rather good estimates of the properties and anisotropy coefficients for our experiment without any fitting parameters. Our results highlight the interplay between the microstructure and physical properties, showing that the physical properties of snow subjected to a temperature gradient cannot be described accurately using only

isotropic parameters such as the density and require more refined information. Furthermore, this study constitutes a detailed database on the evolution of snow properties under a temperature gradient, which can be used as a guideline and a validation tool for snow metamorphism models at the micro- or macroscale.

## 1 Introduction

Natural snowpacks are frequently subjected to temperature gradients induced by their environment. Due to temperature differences in the snowpack, the morphology of snow at the microscale, i.e., the snow microstructure, quickly evolves with time. This metamorphism, called temperature gradient (TG) metamorphism, is mainly characterized by the reorganization of ice along the gradient direction by sublimation of the warmest parts of the grains, water vapor transport across the air pores, and its deposition on the coldest zones of the ice matrix (Yosida et al., 1955; de Quervain, 1973; Colbeck, 1997; Flin and Brzoska, 2008). In terms of snow type, this leads to faceted crystals and depth hoar, which constitute often the weakest layers of the snowpack. Experimental and theoretical studies such as those of Yosida et al. (1955), de Quervain (1973), Akitaya (1974), Marbouty (1980), Colbeck (1983a, b), Fukuzawa and Akitaya (1993) and Satyawali et al. (2008) provide a good base of knowledge on TG metamorphism, with descriptions of the evolution of the snow grains mostly based on photographs. With the development
of X-ray microtomography for snow (Brzoska et al., 1999a; Schneebeli, 2000; Coléou et al., 2001; Lundy et al., 2002; Pinzer and Schneebeli, 2009a; Chen and Baker, 2010), very precise studies related to TG metamorphism are now available. Up to now, two different approaches have been used: the static approach, where the metamorphism of a homogeneous snow slab can be monitored by imaging different impregnated snow samples collected in the slab (Flin and Brzoska, 2008; Srivastava et al., 2010), and the dynamic or in vivo approach, which gives access to the grain to grain evolution of the same snow sample by time-lapse tomography (Schneebeli and Sokratov, 2004; Pinzer and Schneebeli, 2009b; Pinzer et al., 2012). They allow for a better understanding of the mechanisms involved and highlight the impact of snow microstructure on its physical and mechanical properties.

In particular, snow properties are often expressed as functions of snow density such as for the effective thermal conductivity (Yen, 1981; Sturm et al., 1997; Calonne et al., 2011; Löwe et al., 2013) or the intrinsic permeability (Shimizu, 1970; Jordan et al., 1999; Courville, 2010; Zermatten et al., 2011; Calonne et al., 2012), leading to simple parameterizations that can be used to estimate properties in snow-pack models, e.g., Crocus (Brun et al., 1989) and Snow-pack (Lehning et al., 1999). However, Schneebeli and Sokratov (2004) and Satyawali et al. (2008) have shown that during TG metamorphism, the effective thermal conductivity of snow evolves without significant changes in density, but only because of the ice/pore reorganization. Such studies suggest that there is a need to refine the parameterizations of snow properties, at least for snow subjected to temperature gradients. In addition, as recently shown for the effective thermal conductivity (Calonne et al., 2011; Shertzer and Adams, 2011; Riche and Schneebeli, 2013), the intrinsic permeability (Calonne et al., 2012), or the effective vapor diffusion (Calonne et al., 2014), this type of snow exhibits anisotropic behavior and requires more systematic investigations. Recently, Löwe et al. (2013) proposed a refined parameterization of the effective thermal conductivity tensor of snow based on anisotropic second-order bounds. Their results show the importance of taking into account the microstructural anisotropy for the estimation of the effective thermal conductivity during TG metamorphism.

We propose addressing these issues by studying the evolution of snow morphology together with several physical properties during a typical experiment of TG metamorphism. The main objective consists in better understanding the relationships between the snow microstructure and its properties. In this context, our paper focuses on the description of the time evolution of a snow slab of 294 kg m⁻³ subjected to a vertical temperature gradient of 43 K m⁻¹ in a cold room at -4 °C. The temperature and gradient values were chosen to observe a significant but not extreme evolution of the snow in a reasonable time (3 weeks of experiment). Moreover, these experimental conditions are in the range of conditions

frequently encountered by natural alpine snowpacks. Snow specimens were regularly sampled from the snow slab and, after treatment, scanned by X-ray microtomography to obtain a set of 3-D images showing the time evolution of the snow microstructure. Then, computations were performed on the 3-D images to estimate various geometrical and physical properties. Whenever possible, specific attention was paid to assess these properties in the x, y, and z directions; z being along the direction of gravity and of the macroscopic temperature gradient. In addition, following the approach of Löwe et al. (2013), we present two anisotropic analytical estimates for the determination of the physical properties of snow based on the knowledge of basic microstructural information (porosity, correlation lengths in the x, y, and z directions). This offers interesting possibilities for the improvement of the parameterizations of snow properties.

The new contributions of this study lie in the following points: (i) a wide range of snow properties (mean and Gaussian curvature distributions, directional correlation lengths, specific surface area, air and ice tortuosities, intrinsic permeability, effective thermal conductivity) are investigated during the same experiment; (ii) the time evolution of most properties computed in the x, y, and z directions is provided, allowing monitoring the anisotropy of properties with time; and (iii) the physical properties computed on 3-D images are compared with those determined by anisotropic analytical estimates based on basic microstructural properties.

## 2 Materials and methods

### 2.1 Experimental setup and 3-D images

Natural snow was collected at Chamrousse (1800 m, French Alps) on 22 February 2011 and stored at -20 °C for 2 weeks. This snow was then sieved in a cold room at -5 °C to obtain a horizontal snow slab of 100 cm length, 50 cm width, and 14 cm height, composed of rounded grains (RG; Fierz et al., 2009) at 300 ± 15 kg m⁻³ (result from macroscopic density measurements). The snow slab was confined at the base and the top between two copper plates whose temperature was controlled by a thermoregulated fluid circulation. The whole system was insulated with 8 cm thick polystyrene plates. An illustration of the experimental setup is given in Fig. 1. Isothermal conditions at -5 °C were first applied to the snow slab during 24 h. This aimed at sintering snow grains whose bonds may have been destroyed by sieving. During the following 3 weeks, the temperature of the cold room was held at -4 °C and the upper and lower copper plates were maintained at -7 and -1 °C, respectively, generating a steady vertical temperature gradient of 43 K m⁻¹ through the snow slab.

The snow slab was sampled using a cylindrical core drill approximately every 3 days over the 3 weeks, leading to seven samples in total at the end of the experiment. Macro
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Photograph in cold room of the apparatus designed to control and monitor temperatures at the top and bottom of a snow slab. The front and side vertical polystyrene plates were removed from the device for visualization purposes.</figcaption></figure>

photographs of snow particles were also taken to characterize
snow type. During the sampling operation, the temperature
of the cold room was temporarily held at -7 °C (temperature
of the upper copper plate) in order to minimize the change of
boundary conditions of the snow slab. The polystyrene plates
and the upper copper plate were then temporarily removed in
order to access the snow slab. The samples were taken in
the middle height of the layer and at a minimum distance
of 5 cm from edges and from regions already sampled. The
air gap created by the sampling was systematically refilled
with freshly sieved snow to prevent strong modifications of
the thermal field of the snow slab. Immediately after sam-
pling, each snow specimen was put in a plastic box and im-
pregnated with 1-chloronaphthalene. This organic product, in
liquid state above -15 °C, was poured along the box walls,
slowly filling the open pores of snow. Then, the sample was
frozen in an iso-octane bath cooled by dry ice (–78 °C) to al-
low for the solidification of the 1-chloronaphthalene. The im-
pregnation is required to stop the metamorphism of the snow
microstructure and consolidate the snow sample for further
machining processes. The absorption properties of the ice,
air and 1-chloronaphthalene ensure a good contrast between
these three components for the X-ray tomographic acquisi-
tion. Cylindrical snow cores were extracted from the sam-
ples by machining with a press drill which is mounted on a
lathe and operated in a cold room at -30 °C. Each snow core
was then glued to the upper part of a copper sample holder
by a droplet of 1-chloronaphthalene and sealed into a Plexi-
glas cap. The prepared samples were finally stored at -20 °C
until the tomographic acquisitions.

Each core was scanned using the conical X-ray microt-
mograph of the 3SR lab set with an acceleration voltage of
75 kV and a current of 100 µA. As this microtomograph op-
erates in an ambient temperature room, the snow core was
placed in a specially designed cryogenic cell composed of

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Illustration of the cryogenic cell used during the tomographic acquisition.</figcaption></figure>

a Peltier module, which maintains a regulated temperature
of around -30 °C at the bottom part of the copper sample
holder. Figure 2 shows a schematic of this cell. A continuous
dry and cold air circulation between the sample holder and
the double-wall Plexiglas chambers of the cell prevents the
deposition or condensation of water vapor on their sides. In
addition, the heat generated by the Peltier module is dissi-
pated by water circulation. The whole system is able to ro-
tate 360° during the acquisition. Each tomographic acqui-
sition lasted around 2 h during which 1200 radiographs of
the entire impregnated snow sample were taken. Horizontal
cross sections of the sample were reconstructed from radio-
graphs using DigiXCT<fnref n="1" /> software. Image processing was then
applied to the grayscale, reconstructed images to obtain bi-
nary images representative of the ice-pore arrangement (Flin
et al., 2003; Lesaffre et al., 2004; Hagenmuller et al., 2013).
The method used consists of the following steps: (i) remov-
ing the remaining air bubbles and associating them with the
pore phase, (ii) smoothing and thresholding, and (iii) visual
verifications and 3-D post processing. One can refer to the
section “Three or more materials”, pp. 862–863, of Hagen-
muller et al. (2013) for detailed information on the exact pro-
cedure applied. We finally obtained seven binary 3-D images,
extracted in the middle of the whole reconstructed volumes
and showing the microstructural evolution of the snow slab
with time. The 3-D images have a voxel size of 8.4 or 9.7 µm
and a volume size of 5.9³, 9.2³ or 9.7³ mm³. Detailed infor-
mation for each image is given in Table 1.

## 2.2 Computation of structural properties

2.2.1 Density

Snow porosity φ (dimensionless), also called volume fraction
of air, was estimated from 3-D images using a standard voxel

<footnote n="1">¹DigiXCT: http://www.digisens3d.com/en/</footnote>
<table><caption>**Table 1.** Microstructural and physical properties computed from 3-D images. Snow types are given according to the international classification (Fierz et al., 2009).</caption><thead><tr><th></th><th></th><th>0A</th><th>1A</th><th>2A</th><th>3A</th><th>4A</th><th>5G</th><th>7G</th></tr></thead><tbody><tr><td>Name</td><td></td><td>RG</td><td>RG</td><td>FC</td><td>DH</td><td>DH</td><td>DH</td><td>DH</td></tr><tr><td>Snow type</td><td></td><td>RG</td><td>RG</td><td>FC</td><td>DH</td><td>DH</td><td>DH</td><td>DH</td></tr><tr><td>Length size of 3-D image (mm)</td><td></td><td>5.9</td><td>5.9</td><td>5.9</td><td>5.9</td><td>5.9</td><td>9.7</td><td>9.2</td></tr><tr><td>Voxel size of 3-D image (µm)</td><td></td><td>8.4</td><td>8.4</td><td>8.4</td><td>8.4</td><td>8.4</td><td>9.7</td><td>9.7</td></tr><tr><td>Time under temperature gradient (h)</td><td></td><td>0</td><td>73</td><td>144</td><td>217</td><td>313</td><td>409</td><td>500</td></tr><tr><td>Density ρ<sub>s</sub> (kg m<sup>-3</sup>)</td><td></td><td>315</td><td>275</td><td>283</td><td>275</td><td>315</td><td>286</td><td>310</td></tr><tr><td>Porosity φ (-)</td><td></td><td>0.657</td><td>0.700</td><td>0.692</td><td>0.700</td><td>0.656</td><td>0.688</td><td>0.622</td></tr><tr><td>Mean curvature H - upward surfaces</td><td>Average (mm<sup>-1</sup>)</td><td>4.6</td><td>4.3</td><td>3.4</td><td>2.3</td><td>1.0</td><td>1.4</td><td>0.8</td></tr><tr><td></td><td>Standard deviation (mm<sup>-1</sup>)</td><td>8.3</td><td>9.0</td><td>9.6</td><td>8.9</td><td>8.5</td><td>8.1</td><td>8.4</td></tr><tr><td>Mean curvature H - downward surfaces</td><td>Average (mm<sup>-1</sup>)</td><td>4.7</td><td>4.5</td><td>3.2</td><td>2.5</td><td>1.2</td><td>1.5</td><td>0.9</td></tr><tr><td></td><td>Standard deviation (mm<sup>-1</sup>)</td><td>8.4</td><td>9.2</td><td>9.1</td><td>9.6</td><td>9.3</td><td>8.8</td><td>9.2</td></tr><tr><td>Gaussian curvature K - all surfaces</td><td>Average (mm<sup>-2</sup>)</td><td>-81.4</td><td>-31.2</td><td>-24.0</td><td>-19.3</td><td>-18.8</td><td>-19.0</td><td>-18.3</td></tr><tr><td></td><td>Standard deviation (mm<sup>-2</sup>)</td><td>284.7</td><td>246.8</td><td>228.7</td><td>208.2</td><td>194.0</td><td>173.9</td><td>184.1</td></tr><tr><td>SSA (m<sup>2</sup> kg<sup>-1</sup>)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>x</td><td>27.4</td><td>22.8</td><td>20.8</td><td>18.2</td><td>15.2</td><td>15.0</td><td>13.5</td></tr><tr><td></td><td>y</td><td>26.5</td><td>22.6</td><td>20.5</td><td>18.2</td><td>15.2</td><td>14.8</td><td>13.4</td></tr><tr><td></td><td>z</td><td>29.2</td><td>24.7</td><td>21.0</td><td>18.2</td><td>14.6</td><td>14.9</td><td>13.3</td></tr><tr><td>SSA (m<sup>2</sup> kg<sup>-1</sup>)</td><td>(x + y + z)/3</td><td>27.7</td><td>23.4</td><td>20.8</td><td>18.2</td><td>15.0</td><td>14.9</td><td>13.4</td></tr><tr><td>Correlation length l<sub>c</sub> (µm)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>x</td><td>70</td><td>95</td><td>104</td><td>128</td><td>146</td><td>160</td><td>181</td></tr><tr><td></td><td>y</td><td>73</td><td>98</td><td>109</td><td>133</td><td>147</td><td>160</td><td>182</td></tr><tr><td></td><td>z</td><td>66</td><td>91</td><td>112</td><td>143</td><td>174</td><td>202</td><td>225</td></tr><tr><td>Air tortuosity τ<sub>a</sub> (-)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>x</td><td>0.73</td><td>0.77</td><td>0.73</td><td>0.71</td><td>0.66</td><td>0.66</td><td>0.63</td></tr><tr><td></td><td>y</td><td>0.74</td><td>0.77</td><td>0.74</td><td>0.71</td><td>0.65</td><td>0.66</td><td>0.63</td></tr><tr><td></td><td>z</td><td>0.71</td><td>0.76</td><td>0.76</td><td>0.76</td><td>0.73</td><td>0.73</td><td>0.71</td></tr><tr><td>Ice tortuosity τ<sub>i</sub> (-)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>x</td><td>0.21</td><td>0.09</td><td>0.08</td><td>0.08</td><td>0.14</td><td>0.12</td><td>0.16</td></tr><tr><td></td><td>y</td><td>0.22</td><td>0.09</td><td>0.09</td><td>0.08</td><td>0.13</td><td>0.12</td><td>0.15</td></tr><tr><td></td><td>z</td><td>0.19</td><td>0.12</td><td>0.16</td><td>0.15</td><td>0.25</td><td>0.21</td><td>0.27</td></tr><tr><td>Thermal conductivity k (W m<sup>-1</sup> K<sup>-1</sup>)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>x</td><td>0.22</td><td>0.14</td><td>0.14</td><td>0.14</td><td>0.19</td><td>0.16</td><td>0.20</td></tr><tr><td></td><td>y</td><td>0.23</td><td>0.14</td><td>0.14</td><td>0.14</td><td>0.18</td><td>0.17</td><td>0.19</td></tr><tr><td></td><td>z</td><td>0.21</td><td>0.15</td><td>0.18</td><td>0.18</td><td>0.25</td><td>0.21</td><td>0.26</td></tr><tr><td>Permeability K (× 10<sup>-9</sup> m<sup>2</sup>)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>x</td><td>0.76</td><td>1.75</td><td>1.95</td><td>2.98</td><td>2.71</td><td>3.77</td><td>3.92</td></tr><tr><td></td><td>y</td><td>0.78</td><td>1.74</td><td>1.95</td><td>2.89</td><td>2.61</td><td>3.70</td><td>3.95</td></tr><tr><td></td><td>z</td><td>0.70</td><td>1.64</td><td>2.02</td><td>3.08</td><td>3.18</td><td>4.39</td><td>4.84</td></tr></tbody></table>

counting algorithm. Snow density $\rho_s$ (in kg m$^{-3}$) was simply deduced from $\phi$ as $\rho_s = \rho_i(1 - \phi)$, where $\rho_i$ is the ice density equal to 917 kg m$^{-3}$.

### 2.2.2 Specific surface area SSA

The specific surface area estimates along the x, y and z directions, denoted by SSAx, SS Ay and SSAz (in m$^2$ kg$^{-1}$), were computed from 3-D images, using a stereologic method (Arakawa et al., 2009; Flin et al., 2011):

$$SSA_x = \frac{2N_x}{L\rho_s}, \quad SSA_y = \frac{2N_y}{L\rho_s}, \quad \text{and} \quad SSA_z = \frac{2N_z}{L\rho_s}, \quad (1)$$

where $N_x$, $N_y$, and $N_z$ are the total number of intersections between air and ice along parallel testing lines in the x, y, and z directions, respectively, through the entire volume, and L is the total length of the testing lines (in m). We recall that

the z direction corresponds to the direction of gravity and of the macroscopic temperature gradient. In the following, we use the vector $SSA = (SSA_x, SSA_y, SSA_z)$, where SSAz is called the vertical component while the average value of SSAx and SSAy, noted SSAxy, is called the horizontal component. The orientation of this vector in the (x, y, z) coordinate system is thus a way to estimate the degree of anisotropy of the snow surfaces. In addition, averaging the three components of $SSA$ yields a precise estimate of the usual scalar SSA (Berryman, 1998; Flin et al., 2011) provided x, y, and z are aligned with the potential anisotropy axes of the sample.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Two-point probability function $S_2(r)$ for the whole range of $r$ (left panel) and magnified (right panel) in the $x$ (blue), $y$ (green) and $z$ directions (red), obtained from the 3-D image referred to as 7G in Table 1. Solid and dashed lines correspond, respectively, to the values computed from Eq. (3) and to the results of the expression $S_2(r_β) = (φ - φ^2) \exp(-r_β/l_{cβ}) + φ^2$ with $β = (x, y, z)$.</figcaption></figure>

### 2.2.3 Two-point probability function and correlation lengths

At a given time, within 3-D images of snow, we can define the following characteristic function of the air phase:

$$I^a(\mathbf{x}) = \begin{cases} 1 & \text{if } \mathbf{x} \text{ lies in the air phase} \\ 0 & \text{if } \mathbf{x} \text{ lies in the ice phase,} \end{cases}$$

where $\mathbf{x}$ is a position vector within the sample. The one- and two-point probability functions for the air phase are then defined as

$$S_1 = \langle I^a(\mathbf{x}) \rangle, \qquad (2)$$

$$S_2(\mathbf{r}) = \langle I^a(\mathbf{x}) I^a(\mathbf{x} + \mathbf{r}) \rangle, \qquad (3)$$

where $\mathbf{r}$ is a vector oriented in the $x$, $y$ or $z$ direction of the image and the angular brackets denote the volume average. $S_2(\mathbf{r})$ is also called the two-point correlation function or the autocorrelation function. For statistically homogeneous media, $S_1$ is simply equal to the porosity ($\phi$) and $S_2$ depends on $\mathbf{r}$. In general, $S_2$ has the following asymptotic properties (Torquato, 2002):

$$S_2(\mathbf{r} = 0) = S_1 = \phi, \qquad (4)$$

$$\lim_{\mathbf{r} \to \infty} S_2(\mathbf{r}) = \phi^2. \qquad (5)$$

As an illustration, dashed lines in Fig. 3 show the two-point probability function computed over a 3-D image of a snow sample after 500 h of metamorphism (referred as image 7G in Table 1) along the $x$, $y$ and $z$ directions in blue, green and red, respectively. As proposed by Löwe et al. (2011, 2013), by fitting the $S_2(\mathbf{r})$ function along the coordinate axes $β = (x, y, z)$ to an exponential $S_2(r_β) = (φ - φ^2) \exp(-r_β/l_{cβ}) + φ^2$ (solid lines in Fig. 3), one obtains a correlation length $l_{cβ}$ (in µm) in the $x$, $y$, and $z$ directions noted $l_{c_x}$, $l_{c_y}$ and $l_{c_z}$, respectively. $l_{c_z}$ is called the vertical component. The average value of $l_{c_x}$ and $l_{c_y}$ is called the horizontal component and noted $l_{c_{xy}}$. The vector $l_c = (l_{c_x}, l_{c_y}, l_{c_z})$ is often used to characterize

the typical sizes of the heterogeneities in the microstructure, i.e., to define the characteristic lengths of an ice grain and a pore without distinction.

### 2.2.4 Mean and Gaussian curvatures

At a given point on the surface of a 3-D object, the shape is characterized by two principal curvatures, $κ_1$ and $κ_2$, which correspond to the maximum and minimum values of normal curvature at this point. Negative, zero and positive values of $κ_1$ and $κ_2$ define concave, flat and convex lines on the ice surface, respectively. The mean curvature $H$ (in m⁻¹) and Gaussian curvature $K$ (in m⁻²) are used to define the surface geometry with

$$\mathcal{H} = \frac{\kappa_1 + \kappa_2}{2}, \qquad (6)$$

$$\mathcal{K} = \kappa_1 \kappa_2. \qquad (7)$$

The signs of the mean and Gaussian curvatures characterize the surface shape. For the mean curvature, negative, zero and positive values correspond to concave, flat and convex surfaces of ice, respectively. For the Gaussian curvature, negative, zero, and positive values represent saddle-shaped surfaces (typically bonds between ice grains), flat or cylindrical surfaces, and dome-shaped surfaces (convex or concave), respectively.

Many techniques have been proposed to estimate mean and Gaussian curvatures on either triangular or digital surfaces (Brzoska et al., 1999b; Nishikawa et al., 2001; Rieger et al., 2002; Zhang et al., 2002; Ogawa et al., 2006; Pottmann et al., 2009, e.g.,). Curvature estimations usually imply specific accuracy issues since these estimators are particularly sensitive to noise and digitization effects. This is mainly due to the fact that curvatures are second-order derivatives obtained on a discrete grid. In our approach, the mean ($H$) and Gaussian ($K$) curvatures are adaptively computed from the largest relevant neighborhoods, limiting their digitization noise (see Flin et al., 2004; Brzoska et al., 2007; Wang et al., 2012 for details). In short, we rely on the following definition
for $\mathcal{H}$ and $\mathcal{K}$ (Sethian, 1999), where the mean curvature can be defined as the divergence of the normal vector field $\mathbf{n}(p)$ at point $p$

$$ \mathcal{H}(p) = \frac{\nabla \cdot \mathbf{n}(p)}{2}, \qquad (8) $$

and the Gaussian curvature as

$$ \mathcal{K}(p) = \frac{M(p)}{N(p)}, \qquad (9) $$

with

$$ \begin{aligned} M(p) = & d_{,x}^2 (d_{,yy} d_{,zz} - d_{,yz}^2) \\ & + d_{,y}^2 (d_{,xx} d_{,zz} - d_{,xz}^2) \\ & + d_{,z}^2 (d_{,xx} d_{,yy} - d_{,xy}^2) \\ & + 2d_{,x} d_{,y} (d_{,xz} d_{,yz} - d_{,xy} d_{,zz}) \\ & + 2d_{,y} d_{,z} (d_{,xy} d_{,xz} - d_{,yz} d_{,xx}) \\ & + 2d_{,x} d_{,z} (d_{,xy} d_{,yz} - d_{,xz} d_{,yy}) \end{aligned} \qquad (10) $$

and

$$ N(p) = (d_{,x}^2 + d_{,y}^2 + d_{,z}^2)^2, \qquad (11) $$

where $d$ is the signed distance map at $p$ and $d_{,x}$, $d_{,y}$ and $d_{,z}$ denote the partial derivatives of $d$ along the $x$, $y$ and $z$ coordinates, respectively. For $\mathcal{H}$ (Eq. 8), the normal vector field $\mathbf{n}(p)$ could also have been expressed as a partial derivative of $d$. However, we use a specific normal vector estimation as proposed by Flin et al. (2005). Such an approach is based on an adaptive computation of the normal vector field using volumetric information obtained from the signed distance map. This gives us a precise estimation of $\mathcal{H}$ while decreasing the sensitivity of this formula to digitization effects (see Flin et al., 2004). For $\mathcal{K}$, we simply use Eqs. (9)–(11) where local estimations $\mathcal{K}(p)$ are averaged on the neighborhoods obtained for the adaptive analysis of $\mathcal{H}$ (Wang et al., 2012).

## 2.3 Computations of tortuosity, effective thermal conductivity and permeability tensors

The full 3-D tensors of tortuosity $\boldsymbol{\tau}$ (dimensionless) of effective thermal conductivity $\mathbf{k}$ (in W m$^{-1}$ K$^{-1}$) and of intrinsic permeability $\mathbf{K}$ (in m$^2$) were computed from 3-D images. For that purpose, specific boundary value problems arising from the homogenization process (Auriault et al., 2009; Calonne et al., 2014) have been numerically solved on representative elementary volumes (REVs) extracted from 3-D images of snow by using the software Geodict<fnref n="2" />, based on a finite difference method (Thoemen et al., 2008). We define the REV with a side length $l$ by $\Omega$ wherein $\Omega_i$ and $\Omega_a$ are the domains occupied by the ice and the air, respectively, and where $\Gamma$ denotes the common boundary.

To compute the effective thermal conductivity tensor $\mathbf{k}$, the following boundary value problem was solved (Auriault et al., 2009; Calonne et al., 2011, 2014):

$$ \nabla \cdot (k_i (\nabla \mathbf{t}_i + \mathbf{I})) = \mathbf{0} \quad \text{within } \Omega_i, \qquad (12) $$

$$ \nabla \cdot (k_a (\nabla \mathbf{t}_a + \mathbf{I})) = \mathbf{0} \quad \text{within } \Omega_a, \qquad (13) $$

$$ \mathbf{t}_i - \mathbf{t}_a = \mathbf{0} \quad \text{on } \Gamma, \qquad (14) $$

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

$$ \frac{1}{|\Omega|} \int_{\Omega} (\mathbf{t}_a + \mathbf{t}_i) \, d\Omega = \mathbf{0}, \qquad (16) $$

where $\mathbf{I}$ is the identity tensor, $\mathbf{n}$ is the outward vector normal to the ice surface and the two periodic vectors $\mathbf{t}_i$ and $\mathbf{t}_a$ are unknown. These two vectors characterize the fluctuations of the temperature field in the ice and air phase which are induced by a given macroscopic gradient of temperature $\nabla T$ applied on the REV. Finally, as in Calonne et al. (2011), $k_a = 0.024 \text{ W m}^{-1} \text{ K}^{-1}$ and $k_i = 2.107 \text{ W m}^{-1} \text{ K}^{-1}$ stand for the thermal conductivity of air and ice at 271 K, respectively. The effective thermal conductivity tensor $\mathbf{k}$ is defined as

$$ \mathbf{k} = \frac{1}{|\Omega|} \left( \int_{\Omega_a} k_a (\nabla \mathbf{t}_a + \mathbf{I}) \, d\Omega + \int_{\Omega_i} k_i (\nabla \mathbf{t}_i + \mathbf{I}) \, d\Omega \right). \qquad (17) $$

The tortuosity tensor of the air phase $\boldsymbol{\tau}_a$ and of the ice phase $\boldsymbol{\tau}_i$ are obtained by solving the same above boundary value problem (Eqs. 12–16) assuming that $k_a = 1$ and $k_i = 0$ for $\boldsymbol{\tau}_a$, and $k_a = 0$ and $k_i = 1$ for $\boldsymbol{\tau}_i$. These tensors are defined as

$$ \boldsymbol{\tau}_a = \frac{1}{|\Omega_a|} \int_{\Omega_a} (\nabla \mathbf{t}_a + \mathbf{I}) \, d\Omega, \quad \boldsymbol{\tau}_i = \frac{1}{|\Omega_i|} \int_{\Omega_i} (\nabla \mathbf{t}_i + \mathbf{I}) \, d\Omega. \qquad (18) $$

Let us remark that for snow the air tortuosity is simply linked to the effective diffusion tensor for the water vapor by $\mathbf{D} = \phi D_m \boldsymbol{\tau}_a$, where $D_m$ (in m$^2$ s$^{-1}$) is the molecular diffusion coefficient of the vapor in air at the pore scale. If we assume that the porous medium consists of an equivalent tortuous capillary of total length $l'$, in contrast with the REV length $l$, it can be shown that, by definition, $0 < \tau_a \propto (l/l')^2 < 1$ (Bear, 1972). Consequently, $\tau_a$ tends toward 0 or 1 when the air structure is highly tortuous or straight, respectively (the same is true for $\tau_i$ and the ice structure). The tortuosity is also often defined as $\tau_f \propto (l'/l)$ (Kaempfer et al., 2005), so that our tortuosity definition corresponds to the inverse of $\tau_f^2$.

The tensor $\mathbf{K}$ of intrinsic permeability was obtained by solving the following boundary value problem (Calonne et al., 2012):

$$ \mu \Delta \mathbf{v} - \nabla \tilde{p} - \nabla p = \mathbf{0} \quad \text{within } \Omega_a, \qquad (19) $$

$$ \nabla \cdot \mathbf{v} = 0 \quad \text{within } \Omega_a, \qquad (20) $$

$$ \mathbf{v} = \mathbf{0} \quad \text{on } \Gamma, \qquad (21) $$

<footnote n="2">²Geodict: http://www.geodict.de</footnote>
where $v$ and $\tilde{p}$ (with $\langle \tilde{p} \rangle = 0$) are the periodic unknowns which represent, respectively, the fluid velocity and the pressure fluctuation in a REV induced by a given macroscopic gradient of pressure $\nabla p$. $\mu$ is the dynamic viscosity of air (in Pa s). It can be shown that $v = -(1/\mu) b \nabla p$, where $b$ is a second order tensor which characterizes the variation of the fluid velocity at the pore scale over a REV induced by a given macroscopic gradient of pressure (Auriault et al., 2009). Consequently, the permeability tensor is defined as

$$ K = \frac{1}{|\Omega|} \int_{\Omega_a} bd\Omega. \qquad (22) $$

As the non-diagonal terms of the tensors $\tau_i$, $\tau_a$, $k$ and $K$ are negligible compared to the diagonal terms (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 and of the temperature gradient), we only focus on the latter ones. In the following, we denote as $\star_x$, $\star_{xy}$ and $\star$ the vertical component, the average of the two horizontal components ($\star_x$ and $\star_{xy}$), and the average of the three components of any tensor $\star$ (i.e., $\tau_i$, $\tau_a$, $k$ or $K$), respectively. Moreover, for the sake of simplicity, $\star_{xy}$ and $\star$ are called the horizontal and the average components, respectively.

## 2.4 Computations of anisotropy coefficients

The anisotropy coefficient $A(\star)$ was computed for each of the microstructural and physical properties mentioned above, except for density and curvatures. This coefficient is defined as the ratio between the vertical component over the horizontal one, such as $A(\star) = \star_z/\star_{xy}$, where $\star = l_{SSA}$, $l_c$, $\tau_i$, $\tau_a$, $k$ or $K$. The property is considered isotropic if it exhibits a coefficient $A(\star)$ close to 1, otherwise the property is anisotropic. Note that, since SSA$_{xy}$ characterizes the vertical surfaces while SSA$_{z}$ describes the horizontal ones, we study the anisotropy coefficient of the vector $l_{SSA} = (1/SSA_x, 1/SSA_y, 1/SSA_z)$ to be consistent with the other coefficients.

## 2.5 Readjustment in density

After sieving, the density of the snow slab exhibited slight spatial inhomogeneities ($300 \pm 15$ kg m$^{-3}$, from macroscopic measurements with a corer). In order to focus only on the evolution of snow properties driven by the temperature gradient, readjusted values of effective thermal conductivity and permeability, $k^r$ and $K^r$, were computed as if the density was homogeneous in the snow slab and equal to $294$ kg m$^{-3}$ (average of the density values computed from 3-D images) using the regression proposed in Calonne et al. (2011, 2012), respectively, as follows:

$$ k^r = \frac{k \times k^{\text{fit}}(\rho_s)}{k^{\text{fit}}(\rho_{294})}, \qquad (23) $$

$$ K^r = \frac{K \times K^{\text{fit}}(\rho_s)}{K^{\text{fit}}(\rho_{294})}, \qquad (24) $$

with $\rho_s$ as the computed snow density, $\rho_{294} = 294$ kg m$^{-3}$, $k^{\text{fit}}(\rho_s) = 2.5 \times 10^{-6} \rho_s^2 - 1.23 \times 10^{-4} \rho_s + 0.024$ and $K^{\text{fit}}(\rho_s) = 3.0 \times r_{es}^2 \exp(-0.0130 \times \rho_s)$ where the equivalent sphere radius $r_{es} = 3/(SSA \times \rho_i)$. In this way, we obtained readjusted thermal conductivity ($K_x^r$, $K_y^r$, $K_z^r$) and permeability ($K_x^r$, $K_y^r$, $K_z^r$) values for a density of $294$ kg m$^{-3}$.

## 2.6 Analytical estimates based on ellipsoidal inclusions

We used two analytical estimates based on ellipsoidal inclusions to estimate the physical properties of snow in the $x$, $y$ and $z$ directions: the self-consistent estimates and the dilute beds of spheroids. These estimates require basic microstructural information, which are the volume fraction of each phase ($\phi$, $1-\phi$) and the inclusion aspect ratio and size. These latter were obtained from the 3-D images, and we chose to describe the inclusion characteristics using the correlation lengths ($l_{c_x}$, $l_{c_y}$, $l_{c_z}$).

### 2.6.1 Self-consistent estimates: effective thermal conductivity and air tortuosity

The snow microstructure is considered here as a macroscopically anisotropic composite, which corresponds to an assemblage of isotropic ellipsoidal inclusions of air and ice with a major axis collinear with the $z$ direction, of same aspect ratio $a/b$, with volume fractions ($\phi$, $1-\phi$) and thermal conductivities ($k_a$, $k_i$), (see Fig. 4). Since the correlation lengths ($l_{c_x}$, $l_{c_y}$, $l_{c_z}$) characterize the typical sizes of the heterogeneities (air and ice without distinction), it seems reasonable, in a first order of approximation, to assume that the aspect ratio is given by $a/b = (l_{c_x} + l_{c_y})/(2l_{c_z})$. According to the self-consistent scheme (Bruggeman, 1935; Hill, 1965; Budiansky, 1965; Willis, 1977; Torquato, 2002), each type of inclusion is successively embedded in a homogeneous equivalent medium, i.e., an infinite matrix whose effective thermal conductivity $k^{sc}$ is the unknown to be calculated, which is a way to capture the connectivity of both phases. The solution of equations for an isolated inclusion then gives an implicit relation which can be solved for this effective property. In the present case, the self-consistent estimate of the effective thermal conductivity $k^{sc}$ verifies the following implicit relation (Torquato, 2002):

$$ \phi (k_a \mathbf{I} - \mathbf{k}^{\text{sc}}) (\mathbf{I} + \mathbf{A}[\mathbf{k}^{\text{sc}}])^{-1} (k_a \mathbf{I} - \mathbf{k}^{\text{sc}}) \mathbf{I} + (1-\phi) (k_i \mathbf{I} - \mathbf{k}^{\text{sc}}) (\mathbf{I} + \mathbf{A}[\mathbf{k}^{\text{sc}}])^{-1} (k_i \mathbf{I} - \mathbf{k}^{\text{sc}}) \mathbf{I} = 0, \qquad (25) $$

with $k_i = 2.107$ Wm$^{-1}$ K$^{-1}$ and $k_a = 0.024$ Wm$^{-1}$ K$^{-1}$, and where $A$ is the depolarization tensor for an ellipsoid in a
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Schematic representation of the microstructure corresponding to the two-point bounds and the self-consistent scheme. Effective thermal conductivity versus ice volume fraction when $b/a = 1.45$: self-consistent estimates (black curves), two-point lower (blue curves) and upper (red curves) bounds.</figcaption></figure>

matrix with an effective thermal conductivity $\mathbf{k}^{\text{sc}}$. When $\mathbf{k}^{\text{sc}}$ is transverse isotropic, the depolarization tensor $\mathbf{A}$ is defined in the $(x, y, z)$ frame as (Giraud et al., 2007; Kushch and Sevostianov, 2014)

$$ \mathbf{A} = \begin{pmatrix} Q & 0 & 0 \\ 0 & Q & 0 \\ 0 & 0 & 1 - 2Q \end{pmatrix}, \qquad (26) $$

with

$$ Q = \frac{1}{2(1 - \gamma^2)} \left( 1 - \frac{\gamma^2}{\sqrt{\gamma^2 - 1}} \tan^{-1} \left( \sqrt{\gamma^2 - 1} \right) \right) \quad (27) $$
if $\gamma \ge 1$,

$$ Q = \frac{1}{2(1 - \gamma^2)} \left( 1 - \frac{\gamma^2}{2\sqrt{1 - \gamma^2}} \ln \left( \frac{1 + \sqrt{1 - \gamma^2}}{1 - \sqrt{1 - \gamma^2}} \right) \right) \quad (28) $$
if $\gamma \le 1$,

where $\gamma$ is linked to the aspect ratio of the ellipsoid and anisotropy ratio of $\mathbf{k}^{\text{sc}}$ by $\gamma = (a/b) \times (k_{xy}^{\text{sc}}/k_z^{\text{sc}})$.

Thus, from Eqs. (25) and (26), the horizontal and vertical components of $\mathbf{k}^{\text{sc}}$ are written as follows:

$$ k_x^{\text{sc}} = k_y^{\text{sc}} = k_{xy}^{\text{sc}} = \frac{-(k_a(\phi - Q) + k_i((1 - \phi) - Q)) - \sqrt{\Delta_{xy}}}{2(Q - 1)}, \quad (29) $$

$$ k_z^{\text{sc}} = \frac{-(k_a(\phi - (1 - 2Q)) + k_i((1 - \phi) - (1 - 2Q))) - \sqrt{\Delta_z}}{2((1 - 2Q) - 1)}, \quad (30) $$

where

$$ \Delta_{xy} = (k_a(\phi - Q) + k_i((1 - \phi) - Q))^2 - 4(Q - 1)Qk_ik_a, \quad (31) $$

$$ \Delta_z = (k_a(\phi - (1 - 2Q)) + k_i((1 - \phi) - (1 - 2Q)))^2 - 4(1 - 2Q) - 1(1 - 2Q)k_ik_a. \quad (32) $$

The self-consistent estimate of the air tortuosity tensor $\boldsymbol{\tau}_a^{\text{sc}}$ can be easily deduced from Eq. (25) with $k_a = 1$ and $k_i = 0$:

$$ \boldsymbol{\tau}_a^{\text{sc}} = \phi^{-1} \mathbf{k}^{\text{sc}} (k_a = 1, k_i = 0), \quad (33) $$

with $\tau_a^{\text{scx}} = \tau_a^{\text{scy}} = \tau_a^{\text{scxy}} = \phi^{-1} k_{xy}^{\text{sc}} (k_a = 1, k_i = 0)$ and $\tau_a^{\text{scz}} = \phi^{-1} k_{z}^{\text{sc}} (k_a = 1, k_i = 0)$. Let us remark that in the particular case described above, the depolarization tensor of Eq. (26) is equal in both phases and Eq. (25) is invariant under the simultaneous interchanges $k_a \leftrightarrow k_i$ and $\phi \leftrightarrow 1 - \phi$, meaning that each phase is treated symmetrically.

As an illustration, Fig. 4 shows the behavior of $k_x^{\text{sc}} = k_y^{\text{sc}}$ and $k_z^{\text{sc}}$ vs. the ice volume fraction for $b/a = 1.45$. In this figure, the two-point bounds for anisotropic composites (Willis, 1977) are also shown. The corresponding microstructure of the lower bounds can be viewed as ellipsoidal inclusions of ice of the same aspect ratio ($b/a$) dispersed within the air matrix, as shown in the upper-left part of Fig. 4. Inversely, for the lower bounds, the microstructure is seen as ellipsoidal inclusions of air of the same aspect ratio ($b/a$) dispersed within the ice matrix (see upper-right part of Fig. 4). As expected, in each direction, the self-consistent estimate lies between the bounds: at low volume fractions of ice, the self-consistent estimates and the lower bounds are very close; conversely, at high volume fractions of ice, the self-consistent estimates are quite similar to the upper bounds. Finally, Fig. 4 clearly shows the anisotropy of the effective thermal conductivity induced by the anisotropy of the microstructure. The anisotropy coefficient $A(\mathbf{k}^{\text{sc}})$ is defined as $k_x^{\text{sc}}/k_z^{\text{sc}}$ and consequently, from Eqs. (29) and (30), is a function of the ratio $k_i/k_a$, the porosity $\phi$, and the aspect ratio $b/a$. In the particular case of Fig. 4, $k_i/k_a \simeq 100$ and $b/a = 1.45$; so $A(\mathbf{k}^{\text{sc}})$ depends only on the porosity and ranges from 1 to 1.6 in the whole range of the ice volume fraction.

## 2.6.2 Dilute beds of spheroids: permeability

The snow is seen as a dilute dispersion of ellipsoids of ice in a matrix of air. The semiaxes of each ellipsoid are defined
as $a = (l_{c_x} + l_{c_y})/4$ and $b = l_{c_z}/2$. It can be shown (Torquato, 2002) that the permeability tensor estimate $\mathbf{K}^{\text{el}}$ is written in the $(x, y, z)$ frame as

$$ \mathbf{K}^{\text{el}} = \frac{2a^2}{9(1-\phi)} \begin{pmatrix} f(b/a) & 0 & 0 \\ 0 & f(b/a) & 0 \\ 0 & 0 & g(b/a) \end{pmatrix}, \quad (34) $$

where

$$ f(b/a) = \frac{3}{16\chi_b^3} \left( (3\chi_b^2 - 1) \ln\left(\frac{1+\chi_b}{1-\chi_b}\right) + 2\chi_b \right) \quad (35) $$
if $b/a \ge 1$

$$ f(b/a) = \frac{3}{8\chi_a^3} \left( (1 + 3\chi_a^2) \tan^{-1}(\chi_a) - \chi_a \right) \quad (36) $$
if $b/a \le 1$,

and

$$ g(b/a) = \frac{3}{8\chi_b^3} \left( (1 + \chi_b^2) \ln\left(\frac{1+\chi_b}{1-\chi_b}\right) - 2\chi_b \right) \quad (37) $$
if $b/a \ge 1$,

$$ g(b/a) = \frac{3}{4\chi_a^3} \left( (\chi_a^2 - 1) \tan^{-1}(\chi_a) + \chi_a \right) \quad (38) $$
if $b/a \le 1$,

where $\chi_a$ and $\chi_b$ are linked to the aspect ratio of the ellipsoid as $\chi_a^2 = -\chi_b^2 = (a/b)^2 - 1$. By definition, this estimation of the permeability does not depend on the spatial arrangement of the ellipsoids and consequently does not capture the real tortuosity of the porous media induced by the connectivity of both air and ice phases. In order to overcome this problem, the following permeability tensor estimate $\mathbf{K}^{\text{di}}$ is proposed such as

$$ K_x^{\text{di}} = K_y^{\text{di}} = K_{xy}^{\text{di}} = \tau_{a_{xy}}^{\text{sc}} K_{xy}^{\text{el}}, \quad K_z^{\text{di}} = \tau_{a_z}^{\text{sc}} K_z^{\text{el}}, \quad (39) $$

where $K_{xy}^{\text{el}} = K_x^{\text{el}} = K_y^{\text{el}}$ and $K_z^{\text{el}}$ are the diagonal components of $\mathbf{K}^{\text{el}}$ in the $(x, y, z)$ frame (see Eq. 34). This relation allows for the recovery of an expression of the permeability similar to the one of Carman-Kozeny (Bear, 1972): $K^{\text{di}} \propto h(\phi) \tau_a d_c^2 = h(\phi) d_c^2 / \tau_f^2$, where $h$ is a function of the porosity and $d_c$ is a characteristic length of the microstructure.

# 3 Results

## 3.1 Time evolution of microstructural and physical properties of snow

Figure 5 illustrates the time evolution of the snow microstructure during the experiment of temperature gradient metamorphism. 3-D images obtained from X-ray tomography are presented together with the corresponding vertical cross section and photograph. The color coding of the 3-D

images corresponds to the mean curvature field. For a better visualization of the faceted shapes, the images are presented “upside down”: the top of the images corresponds to the lowest and warmest side of the physical sample. We observe qualitatively that the initial rounded grains become bigger and more angular and faceted with time. After about 200 h, depth hoar is obtained showing characteristic striations on the surface of grains (see photographs in Fig. 5). At the end of the experiment, the ice structure is preferentially arranged along the vertical direction, i.e., the direction of the temperature gradient, as shown by the cross sections.

All the snow properties computed based on the 3-D images are summarized in Table 1. The time evolution of snow density, specific surface area, correlation length, tortuosity of ice and air, thermal conductivity and permeability are depicted in Fig. 6. The blue, green and red symbols represent the $x$, $y$ and $z$ values of the considered property, respectively. One can observe the following:

- The snow density shows no significant evolution with time and the average value over the experiment is 294 kg m$^{-3}$ (porosity of 0.68). In detail, low variations between 275 and 315 kg m$^{-3}$ are observed from one image to another (Fig. 6a). As explained in Sect. 2.5, these variations reflect the spatial heterogeneity initially present in the sieved snow layer, and not a real-time evolution generated by the temperature gradient conditions.

- The average value of the SSA estimates in the three directions decreases continuously with time from 27.7 to 13.4 m$^2$ kg$^{-1}$. Between 0 and 144 h, the $z$ estimates are slightly higher than the horizontal ones (29.2 vs. 26.9 m$^2$ kg$^{-1}$ at 0 h). After 144 h, values in the three directions become very close to each other (Fig. 6b).

- The values of correlation length increase continuously during the experiment, evolving from 71 to 181 µm in the $x$ and $y$ directions and from 68 to 228 µm in the $z$ direction (Fig. 6c).

- The values of air tortuosity (~0.7) are around 5 times higher than those of ice tortuosity (~0.15) (Fig. 6d). The overall evolution of both properties is low and can be divided in two stages: the ice tortuosity decreases between 0 and 73 h and then slightly increases until the end of the experiment, while the air tortuosity shows the opposite trends during these two equal periods. During the second stage, the $z$ values stand out and become increasingly higher than the horizontal ones for both phases.

- The raw values of the effective thermal conductivity, which are referred as “computed” and depicted by the dashed lines in Fig. 6e, exhibit the same variations as the snow density, showing the strong relationship between these two variables. The values range between 0.14 and 0.26 W m$^{-1}$ K$^{-1}$. The solid lines in Fig. 6e,
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Microstructure evolution during the TG metamorphism. For each stage of the evolution, the following views are given: (i) 3-D images of the snow samples where colors represent the mean curvature of the surfaces, ranging from -36 to +36 mm⁻¹. Convexities, flat shapes and concavities are shown in red, yellow and green, respectively. Images have a size of 3 mm × 3 mm × 3 mm. For a better visualization of the faceted shapes, the images are presented "upside down". The arrows in blue, green and red correspond to the x, y and z directions of the images, respectively. (ii) Vertical cross sections from 3-D images of 5.5 mm × 5.5 mm × 5.5 mm, where ice is in white and air in black. (iii) Photographs of snow grains.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Time evolution of microstructural and physical properties of snow during the whole temperature gradient experiment. Values in the *x*, *y* and *z* directions are given in blue, green and red, respectively.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Time evolution of the anisotropy coefficient of the whole computed properties.</figcaption></figure>

referred to as “readjusted”, show values of thermal conductivity after readjustment at a density of 294 kg m⁻³. As explained in Sect. 2.5, this readjustment is one way to estimate the thermal conductivity without taking into account the influence of density variations. The evolution of readjusted values is thus smoother than the one of computed values, but is as significant. Computed and readjusted values decrease from 0 to 73 h and then continuously increase until the end of the metamorphism, showing an evolution in two stages similar to the one of tortuosities. After 73 h, the z values become much higher than the x and y values with time.

- The computed values of permeability range between 0.70 × 10⁻⁹ and 4.84 × 10⁻⁹ m². They exhibit an opposite evolution to the one of snow density, but this dependence seems less pronounced than for the thermal conductivity. Indeed, despite the influence of density, a significant evolution is still observed during the metamorphism (dashed lines in Fig. 6f). Both computed and readjusted values increase over the experiment and vertical values become higher than the horizontal ones after 217 h.

Figure 7 shows the time evolution of the anisotropy coefficient for five snow properties. Overall, the evolution is the same for all properties: the coefficient increases from values lower than one at the beginning of the experiment to values greater than one at the end. In detail, the magnitude of the coefficients is strongly different from one variable to another. The largest evolution is shown by the coefficient of ice tortuosity $A(\tau_i)$, which increases from 0.86 to 1.90 between 0 and 144 h and then slightly decreases to reach 1.70 at the end of the metamorphism. The anisotropy coefficient of the thermal conductivity $A(k)$ increases from 0.90 to 1.34. The largest increases of $A(\tau_i)$ and $A(k)$ are observed between 0 and 144 h, showing that changes are almost concentrated during

this period for both of the properties, as already pointed out above and shown in Fig. 6. The correlation length and permeability show similar values of anisotropy coefficients, which evolve from 0.96 to 1.25 and from 0.91 to 1.23, respectively. Finally, anisotropy coefficients of air tortuosity and of $l_{SSA}$ are the smallest and evolve during the experiment from 0.97 to 1.13 and from 0.92 to 1.01, respectively.

The distributions of mean curvature of the upward (left panel) and downward (right panel) surfaces of ice are presented in Fig. 8. They are expressed in terms of occurrence ratio, which gives the ratio in percentage of the ice surface area that exhibits a mean curvature located in a particular range of values over the total ice surface area. The time evolution of these distributions is shown by the plots of different colors. The area averaged and the standard deviation of the mean curvature values are given in Table 1. Initially (red color), the distribution of upward and downward surfaces are similar, with a peak of mean curvature located around 6 mm⁻¹ and an occurrence ratio of ~3%. With time, the area-averaged mean curvature decreases gradually, meaning that ice structures tend to become larger. At the end of the experiment (purple color), the upward and downward surfaces exhibit clearly distinct distributions: the peak of mean curvature is now located at 1 mm⁻¹ (occurrence ratio of 4.2%) for the upward ones, and at 0 mm⁻¹ (occurrence ratio of 4.8%) for the downward ones.

Using the same representation, Fig. 9 shows the Gaussian curvature distributions computed from the entire ice-air interface. A log scale was used to allow for a better visualization of the curves. Here, we do not discern between the upward and downward surfaces of ice because the Gaussian curvature distributions of these two cases are similar. Table 1 provides the time evolution of the area-averaged and standard deviation values. All the distributions are centered at 0 mm⁻², but become sharper over time: the maximum
<figure><img src="image_12.png" /><figcaption>Figure 8. Time evolution of the mean curvature distribution computed from the upward (left panel) and downward (right panel) surfaces of the 3-D images. Each curvature class is 0.5 mm⁻¹ wide.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 9. Time evolution of the Gaussian curvature distribution computed from the whole surface of the 3-D images. Each curvature class is 10 mm⁻² wide.</figcaption></figure>

occurrence ratio evolves from 4.5 to 17.3% between 0 and 500 h and the standard deviation continuously decreases. This implies that the proportion of large, flat or cylindrical structures of ice increases. The initial distribution (red color) stands out from the others: it exhibits occurrence ratios which are small between around −50 and 50 mm⁻² and large between around −500 and −50 mm⁻², compared to those of other distributions. The most important changes in Gaussian curvatures of ice surfaces appear thus early after the beginning of the temperature gradient, where high negative values

corresponding to small saddle-shaped structures disappear in favor of low values close to 0 mm⁻², which reflect large and less curved structures. Qualitatively, the same observations can be done from the 3-D images of snow at 0 (left panel) and 500 h (right panel) of metamorphism represented in Fig. 10 by looking at the color maps of Gaussian curvature. In the initial image, we see a lot of green color (negative Gaussian curvature) at bonds between grains, while the final image exhibits mainly yellow-based colors (nearly zero Gaussian curvatures) representing flat and large structures. The seven images of snow represented with the color map of the Gaussian curvature, as well as those of the mean curvature, are available in the auxiliary materials.

## 3.2 Comparisons with analytical estimates

Figure 11 shows the comparison between the time evolution of computed values (dashed lines) and values given by analytical estimates (solid lines) of air tortuosity, effective thermal conductivity and permeability. Computed values are given in the x, y and z directions while the ones given by estimates are described in the z and xy directions. In order to compare estimates and experiments quantitatively, we use the mean of relative differences $E(\star)$ defined as

$$E(\star) = \frac{\star_{\text{from estimate}} - \star_{\text{computed}}}{\star_{\text{computed}}}, \qquad (40)$$

where $\star$ is a component or the anisotropy coefficient of $\tau_a$, k and K. $E(\star)$ is given with the associated standard deviation. For the vertical values of properties, $E(\tau_a) = 6.3 \pm 4.5\%$, $E(k_z) = -22.6 \pm 4.32\%$ and $E(K_z) = -9.6 \pm 10.9\%$, whereas for the horizontal
<figure><img src="image_14.png" /><figcaption>**Figure 10.** The initial (left panel) and final (right panel) 3-D images of the experiment where colors represent the Gaussian curvature of the surfaces, ranging from -781 to +781 mm⁻². Dome-shaped (concave or convex), flat or cylindrical and saddle-shaped surfaces are shown in red, yellow and green, respectively. Images have a size of 3 mm × 3 mm × 3 mm. For a better visualization of the faceted shapes, the images are presented "upside down". The arrows in blue, green and red correspond to the x, y and z directions of the images, respectively.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Time evolution of air tortuosity, thermal conductivity and permeability during the whole experiment. Comparison between values computed from 3-D images in the x, y and z directions (dashed lines in blue, green and red, respectively) and values given by analytical estimates in the xy and z directions (solid lines in blue and red, respectively).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Time evolution of the anisotropy coefficients of air tortuosity, thermal conductivity and permeability. Comparison between coefficients deduced from values computed on 3-D images (dashed lines) and values given by analytical estimates (solid lines).</figcaption></figure>

ones, $E(\tau_{xy}) = 8.1 \pm 5.2\%$, $E(k_{xy}) = -15.6 \pm 6.0\%$ and $E(K_{xy}) = -8.1 \pm 14.7\%$.

Using the same representation, the time evolution of anisotropy coefficients from computed values and values given by estimates of the three properties above is shown in Fig. 12. For the air tortuosity, permeability and thermal conductivity, $E(\mathcal{A}(\tau_a)) = -1.6 \pm 0.9\%$, $E(\mathcal{A}(\mathbf{k})) = -7.9 \pm 7.4\%$, and $E(\mathcal{A}(\mathbf{K})) = -1.1 \pm 4.2\%$, respectively.

# 4 Discussion

## 4.1 Evolution of the microstructure

From the mean curvature distributions (Fig. 8) and the 3-D images (Fig. 5), we observe that the snow microstructure undergoes a strong directional faceting under a temperature gradient. At the beginning of the experiment, the top and base of the ice grains are identical in shape and exhibit mainly convex surfaces, characteristic of rounded grains. With time, these regions become strongly different, showing mostly convex surfaces at the top and faceted surfaces at the bottom of the grains, which is typical of faceted crystals and depth hoar. This observation is in agreement with the results of Flin et al. (2006) and Flin et al. (2007), who show asymmetries in the mean curvature distributions of snow samples submitted to low TGs (3 and $16 \text{ K m}^{-1}$ at $-3 \text{ }^{\circ}\text{C}$) during 3 weeks. This asymmetric behavior can be attributed to the growth and decay properties of the ice crystal at the molecular level, where some sites are energetically more stable for the deposition or sublimation of a water molecule (see TLK or Kossel crystal models in Markov, 1995; Mutaftschiev, 2001). Thus, the convex surfaces that undergo vapor deposition grow by a layer by layer process resulting in the production of facets. Conversely, the sublimation of convex crystals preferentially

leads to the generation of kink and step sites, which results in the rounding of the shapes (see e.g., Knight, 1966, and Flin and Brzoska, 2008, for more details). Our mean curvature results confirm the above considerations: during the experiment, the upward-oriented grain surfaces are warmer than the surrounding air and sublimate, inducing a rounding of the interface. Symmetrically, the downward-oriented surfaces of the ice microstructure, colder than air, undergo vapor deposition and generate facets.

From the temporal evolution of the snow (e.g., Fig. 6), the metamorphism can be divided in two periods: (i) between 0 and 73 h, the microstructure, resulting from equi-temperature conditions, recrystallizes toward a new pattern of structure to adapt to the temperature gradient conditions. In particular, the small but numerous connections between grains sublimate and favor the growth of large structures (see Gaussian curvature distributions in Fig. 9). During this transitional state, the ice structure becomes more tortuous (minimum value of tortuosity over the whole experiment reached at 73 h in Fig. 6d) because of the decrease of the links between grains. (ii) After 73 h, the structure evolves in the continuity of the pattern developed during the first stage (consolidation): grains and connections become increasingly larger especially in the gradient direction, and the ice network becomes gradually less tortuous. Our observations confirm the previous studies of Bradley et al. (1977), Schneebeli et al. (1999) and Schneebeli and Sokratov (2004) concerning the evolution of bonds under TG.

Because the air tortuosity is about 5 times higher than the ice tortuosity (see Fig. 6d), our snow samples can be considered as being comprised of a tortuous skeleton of ice surrounded by large channels of air in the three directions. Our initial z value of ice tortuosity of 0.19 (referred as 0A in Table 1) is quantitatively in good agreement with the work of Kaempfer et al. (2005), who computed from a 3-D image
of rounded grains at 268 kg m⁻³ a vertical ice tortuosity of τᵢ² = 4.4, which corresponds to 0.23 according to our tortuosity definition. In accordance with the observations of Marbouty (1980), Schneebeli and Sokratov (2004), Satyawali et al. (2008) and Pinzer et al. (2012), the snow density remains constant over the experiment, while grains and pores grow continuously as shown by correlation lengths and curvature distributions. It means that the microstructure is rebuilt without a supply or loss of mass. Moreover, the snow becomes more and more anisotropic with a structure elongated in the vertical direction, as illustrated by the anisotropy coefficients of the correlation length and of the ice tortuosity. From a thermodynamic point of view, this microstructural anisotropy can be seen as a way to decrease the vertical temperature gradient by enhancing the heat flow in that direction (Staron et al., 2014).

## 4.2 Link with the physical properties

Heat conduction in snow is mostly due to the conduction of ice which conducts 100 times better than the air. Consequently, the effective thermal conductivity of snow is strongly linked to its density and to the ice tortuosity, as illustrated in Fig. 6. Anisotropy of thermal conductivity is lower than that of ice tortuosity because the air phase taken into account in the computation of conduction reduces this anisotropy. Air conduction cannot be neglected in the computation of snow effective conductivity even if the ratio $k_i/k_a$ is about 88 at 271 K, as already shown by e.g., Calonne et al. (2011).

We observe two distinct stages in the evolution of the effective thermal conductivity (Fig. 6e): a decrease between 0 and 73 h, followed by an increase until the end of the experiment. This observation is in agreement with the result of Schneebeli and Sokratov (2004), who also noticed a two-stage evolution of this property in roughly similar conditions (average temperature of -8 °C, temperature gradient of 100 K m⁻¹, snow density of 268 kg m⁻³). This behavior can be related to the evolution of the microstructure explained in Sect. 4.1 since (i) the destruction of connections between grains during the first period of the temperature gradient limits the heat conduction in snow and (ii) the growth of large ice structures in the second period, in particular in the z direction, enhances the heat transfer.

As already mentioned, the intrinsic permeability is proportional to $h(\phi)\tau_a d_c^2$, where $d_c$ is a characteristic length of the microstructure which could be linked to the correlation length or to the inverse of the specific surface area (Calonne et al., 2012). Consequently, the anisotropy of the permeability depends on both anisotropies of the air tortuosity and of the chosen characteristic length. Since the anisotropy coefficient of permeability, air tortuosity, and correlation length present similar time evolutions, while that of $l_{SSA}$ stands out by remaining close to one (see Fig. 7), our results seem to show that the anisotropy of the permeability is mainly related

to that of the air tortuosity and of the correlation length for snow. The fact that the permeability is related to the density and a characteristic length can explain that this property is less influenced by the density variations than the thermal conductivity as described in Sect. 3.1, but undergoes the same continuous increase over the experiment that is observed for the characteristic length.

Figure 6 shows an overall increase of the effective thermal conductivity and permeability at a constant density, underlying the role of the microstructure rearrangement (Schneebeli and Sokratov, 2004; Satyawali et al., 2008). Moreover, a significant vertical anisotropy is observed for both variables with coefficients until ~ 1.3. Thus, for snow subjected to a temperature gradient, the first approximation of expressing the physical properties depending on the density and on a single characteristic length for the permeability (e.g., Shimizu, 1970; Yen, 1981; Sturm et al., 1997; Calonne et al., 2011, 2012) should be improved by introducing for instance new microstructural parameters which could reflect the anisotropy phenomena.

## 4.3 Estimates of physical properties

The computed effective properties of snow were compared with analytical estimates based on basic information of the microstructure such as the snow density and the anisotropy of heterogeneities (air and ice) through the correlation lengths. Even if the analytical estimates fail to reflect all the details of the microstructure (e.g., bonds between grains), our results clearly show that they capture the overall evolution of the considered properties with a mean of relative differences ranging from -22.6 to 6.3 % (see Sect. 3.1). These estimates do not include any fitting parameter; in contrast to the parameterization of Löwe et al. (2013) where the lower bounds equation (see Fig. 4) was adjusted to the computed data of thermal conductivity by introducing two parameters. The fitting parameters are necessary since the lower bounds equation does not take into account the connectivity of both phases, in contrast to the self-consistent estimates. We compared the self-consistent estimates and the parameterization of Löwe et al. (2013) quantitatively; overall, both estimates are in agreement for the present data.

In both of the presently proposed estimates, the time evolution of physical properties during the metamorphism mainly depends on one parameter $Q$ which is associated with the shape heterogeneity (pore/grain). In our case, the use of the correlation length to define $Q$ seems relevant. Of course, future comparisons between these estimates and other sets of data are needed in order to validate this approach in a large density range and other values of temperature gradients. Let us remark that these self-consistent estimates can be improved by introducing a finer description of the snow microstructure, i.e., by considering the ice skeleton or the pore network as a complex assemblage of several classes of ellipsoidal inclusions with a given orientation, a given aspect
ratio and a given volume fraction. However, these improve-
ments will require developing specific algorithms in order to
compute, from the 3-D images, the geometric pore size/shape
distribution using an anisotropic structural element like an el-
lipsoid. These computations are not straightforward. Finally,
other relevant estimates of these effective properties can be
probably obtained by introducing (i) other microstructural in-
formation in bounds or exact contrast expansions by comput-
ing three- or four-point probability functions on 3-D images
(Torquato, 1991, 1998, 2002; Roberts and Garboczi, 2002) or
(ii) refined variables allowing us to describe the evolution of
grain shapes and contacts between grains (mean and Gaus-
sian curvatures, fabric tensor) during the metamorphism.

**5 Conclusions**

An experiment of TG metamorphism was performed in a
cold room in order to monitor the evolution of microstruc-
tural and physical properties of snow over time: seven snow
samples were collected, at regular time intervals over 3
weeks, from a snow slab subjected to a vertical temperature
gradient of 43 K m⁻¹. Each sample was scanned by X-ray to-
mography to obtain a series of 3-D images showing the mi-
crostructural evolution of the snow slab during temperature
gradient metamorphism. Numerical computations were per-
formed on these 3-D images to estimate various geometrical
and physical properties of snow, such as the directional corre-
lation lengths, the mean and Gaussian curvatures, the specific
surface area, the air and ice tortuosities, the effective thermal
conductivity and the intrinsic permeability. When possible,
such quantities were evaluated in the x, y and z directions to
monitor the evolution of their anisotropy.

The main results concerning the TG experiment can be summarized as follows: (i) as shown by many other studies, density remained almost constant during the whole experiment. (ii) The grains were continuously faceting at their bottom parts (deposition) while the upper parts underwent rounding due to ice sublimation. (iii) Overall, grain growth and neck growth were observed during the metamorphism. (iv) The intrinsic permeability, linked mainly to the air phase, increased continuously. (v) The snow microstructure evolved in two stages: a short period of strong modifications of the ice structure due to the TG initiation (0–73 h), where the ice tortuosity and the thermal conductivity decreased, followed by a stage of consolidation of this new structure (73–500 h), where the above properties increased gradually. (vi) The anisotropy coefficients of all properties increased during the metamorphism, with larger values in the gradient direction.

These results highlight the strong interplay between the
microstructure and physical properties of snow, and con-
firm that the density alone, or any isotropic quantity, is not
sufficient to describe the time evolution of physical prop-
erties during a TG metamorphism. To solve this problem,
we applied analytical anisotropic estimates (self-consistent

estimates and dilute beds of spheroids) using microstructural
parameters (directional correlation lengths) that reflect the
general shape of heterogeneities (size, anisotropy). The pro-
posed analytical estimates, whose results were compared to
those obtained by numerical computations, offer good esti-
mations of the physical properties and anisotropy coefficients
for our time series, without applying any fitting parameters.

In summary, this study presents numerical tools to quan-
titatively monitor snow properties using 3-D images and
provides a detailed database describing snow under a TG
metamorphism. In particular, it provides a quantification of
snow anisotropy, which is a key – but challenging – param-
eter to access directly with physical measurements during
such a process. Used as a guideline or a validation tool, this
database offers new outlooks for the development of micro-
and macroscale snow models.

**The Supplement related to this article is available online**
**at doi:10.5194/tc-8-2255-2014-supplement.**

*Acknowledgements.* We thank Henning Löwe, Martin Schneebeli, and three anonymous reviewers for their contributions in the improvement of this paper. Funding by Météo-France, INSU-LEFE and DigitalSnow (ANR-11-BS02-009-03) is acknowledged. We thank P. Charrier and J. Desrues of the 3SR laboratory, where the 3-D images were obtained. We are also grateful to J. Roulle, J.-M. Panel, P. Puglièse, C. Carmagnola and S. Morin of the CNRM-GAME for their support during the experiment. N. Calonne thanks S. Morin for his co-supervision of her master's internship, during which the TG experiment has been realized. CNRM-GAME/CEN is part of the Labex OSUG@2020 (ANR10 LABX56). The laboratory 3SR is part of the LabEx Tec 21 (Investissements d'Avenir – grant agreement ANR11LABX0030).

Edited by: M. Schneebeli

**References**

Akitaya, E.: Studies of depth hoar, Low. Temp. Sci. Ser. A, 26, 1–67, 1974.

Arakawa, H., Izumi, K., Kawashima, K., and Kawamura, T.: Study on quantitative classification of seasonal snow using specific surface area and intrinsic permeability, Cold. Reg. Sci. Technol., 59, 163–168, doi:10.1016/j.coldregions.2009.07.004, 2009.

Auriault, J.-L., Boutin, C., and Geindreau, C.: Homogenization of coupled phenomena in heterogenous media, Wiley-ISTE, London, 2009.

Bear, J.: Dynamics of Fluids in Porous Media, American Elsevier Publishing Company, Inc., Dover, 1972.

Berryman, J. G.: Planar spatial correlations, anisotropy, and specific surface area of stationary random porous media, J. Appl. Phys., 83, 1685–1693, doi:10.1063/1.366885, 1998.

Bradley, C., Brown, R., and Williams, T.: Gradient metamorphism, zonal weakening of the snow-pack and avalanche initiation, J. Glaciol., 19, 335–342, 1977.
Bruggeman, D.: The calculation of various physical constants of heterogeneous substances, I. The dielectric constants and conductivities of mixtures composed of isotropic substances, Ann. Phys., 24, 636–679, 1935.

Brun, E., Martin, E., Simon, V., Gendre, C., and Coléou, C.: An energy and mass model of snow cover suitable for operational avalanche forecasting, J. Glaciol., 35, 333–342, 1989.

Brzoska, J.-B., Coléou, C., Lesaffre, B., Borel, S., Brissaud, O., Ludwig, W., Boller, E., and Baruchel, J.: 3D visualization of snow samples by microtomography at low temperature, ESRF Newsletter, 32, 22–23, 1999a.

Brzoska, J.-B., Lesaffre, B., Coléou, C., Xu, K., and Pieritz, R. A.: Computation of 3D curvatures on a wet snow sample, Eur. Phys. J.-Appl. Phys., 7, 45–57, doi:10.1051/epjap:1999198, 1999b.

Brzoska, J.-B., Flin, F., and Ogawa, N.: Using Gaussian curvature for the 3D segmentation of snow grains from microtomographic data, in: Physics and Chemistry of Ice, edited by: Kuhs, W., Special Publication Nr. 311, RSC Publishing, Cambridge, UK, Proceedings of the 11th International Conference on the Physics and Chemistry of Ice held at Bremerhaven, 23–28 July 2006, Bremerhaven, Germany, 125–132, 2007.

Budiansky, B.: On the elastic moduli of some heterogeneous materials, J. Mech. Phys. Solids, 13, 223–227, 1965.

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, doi:10.1029/2011GL049234, 2011.

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

Calonne, N., Geindreau, C., and Flin, F.: Macroscopic modeling for heat and water vapor transfer in dry snow by homogenization, J. Phys. Chem. B, doi:10.1021/jp5052535, in press, 2014.

Chen, S. and Baker, I.: Evolution of individual snowflakes during metamorphism, J. Geophys. Res., 115, D21114, doi:10.1029/2010JD014132, 2010.

Colbeck, S. C.: Ice crystal morphology and growth rates at low supersaturations and high temperatures, J. Appl. Phys., 54, 2677–2682, 1983a.

Colbeck, S. C.: Theory of metamorphism of dry snow, J. Geophys. Res., 88, 5475–5482, 1983b.

Colbeck, S. C.: A review of sintering in seasonal snow, CRREL, Hanover, NH, USA, 1997.

Coléou, C., Lesaffre, B., Brzoska, J.-B., Ludwig, W., and Boller, E.: Three-dimensional snow images by X-ray microtomography, Ann. Glaciol., 32, 75–81, doi:10.3189/172756401781819418, 2001.

Courville, Z., Hörhold, M., Hopkins, M., and Albert, M.: Lattice–Boltzmann modeling of the air permeability of polar firn, J. Geophys. Res., 115, F04030, doi:10.1029/2009JF001549, 2010.

de Quervain, M. R.: Snow structure, heat and mass flux through snow, Proceeding of Symposium at Banff 1972, The Role of Snow and Ice Hydrology, IAHS Publ., 107, 203–226, 1973.

Fierz, C., Armstrong, R. L., Durand, Y., Etchevers, P., Greene, E., McClung, D. M., Nishimura, K., Satyawali, P. K., and Sokravov, S. A.: The international classification for seasonal snow on the ground, IHP-VII Technical Documents in Hydrology no. 83, IACS Contribution no. 1, UNESCO-IHP, Paris, France, 2009.

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

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

Flin, F., Brzoska, J.-B., Lesaffre, B., Coléou, C., and Pieritz, R. A.: Three-dimensional geometric measurements of snow microstructural evolution under isothermal conditions, Ann. Glaciol., 38, 39–44, doi:10.3189/172756404781814942, 2004.

Flin, F., Brzoska, J.-B., Coeurjolly, D., Pieritz, R. A., Lesaffre, B., Coléou, C., Lamboley, P., Teytaud, O., Vignoles, G. L., and Delesse, J.-F.: Adaptive estimation of normals and surface area for discrete 3-D objects: application to snow binary data from X-ray tomography, IEEE T. Image Process., 14, 585–596, doi:10.1109/TIP.2005.846021, 2005.

Flin, F., Brzoska, J.-B., Pieritz, R. A., Lesaffre, B., Coléou, C., and Furukawa, Y.: A tentative model for the temperature gradient snow metamorphism and its first validations on X-ray microtomographic data, edited by: S. Aoki, A., Kagoshima, Y., and Suzuki, Y., Proceedings of the 8th International Conference on X-ray Microscopy, XRM2005, Himeji, Japan, 306–308, 2006.

Flin, F., Brzoska, J.-B., Pieritz, R. A., Lesaffre, B., Coléou, C., and Furukawa, Y.: The temperature gradient metamorphism of snow: Model and first validations using X-ray microtomographic images, Proceedings of the 11th International Conference on the Physics and Chemistry of Ice held at Bremerhaven, 23–28 July 2006, Bremerhaven, Germany, 181–189, 2007.

Flin, F., Lesaffre, B., Dufour, A., Gillibert, L., Hasan, A., Rolland du Roscoat, S., Cabanes, S., and Pugliese, P.: On the computations of Specific Surface Area and Specific Grain Contact Area from snow 3D images, in: P. C. I., edited by: Furukawa, Y., Hokkaido University Press, Sapporo, JP, proceedings of the 12th International Conference on the Physics and Chemistry of Ice held at Sapporo, 5–10 September 2010, Sapporo, Japan, 321–328, 2011.

Fukuzawa, T. and Akitaya, E.: Depth-hoar crystal growth in the surface layer under high temperature gradient, Ann. Glaciol., 18, 39–45, 1993.

Giraud, A., Gruescu, C., Do, D. P., Homand, F., and Kondo, D.: Effective thermal conductivity of transversely isotropic media with arbitrary oriented ellipsoidal inhomogeneities, Int. J. Solids Struct., 44, 2627–2647, 2007.

Hagenmuller, P., Chambon, G., Lesaffre, B., Flin, F., and Naaim, M.: Energy-based binary segmentation of snow microtomographic images, J. Glaciol., 59, 859–873, doi:10.3189/2013JoG13J035, 2013.

Hill, R.: A self consistent mechanics of composite materials, J. Mech. Phys. Solids, 13, 213–222, 1965.

Jordan, R., Andreas, E., and Makshtas, A.: Heat budget of snow-covered sea ice at North Pole 4, J. Geophys. Res., 104, 7785–7806, 1999.
Kaempfer, T., Schneebeli, M., and Sokratov, S.: A microstructural approach to model heat transfer in snow, Geophys. Res. Lett., 32, L21503, doi:10.1029/2005GL023873, 2005.

Knight, C. A.: Formation of crystallographic etch pits on ice and its application to the study of hailstones, J. Appl. Meteorol., 5, 710–714, 1966.

Kushch, V. I. and Sevostianov, I.: Dipole moments, property contribution tensors and effective conductivity of anisotropic particulate composites, Int. J. Eng. Sci., 74, 15–34, doi:10.1016/j.ijengsci.2013.08.002, 2014.

Lehning, M., Bartelt, P., Brown, B., Russi, T., Stöckli, U., and Zim-merli, M.: SNOWPACK model calculations for avalanche warn-ing based upon a new network of weather and snow stations, Cold Reg. Sci. Technol., 30, 145–157, 1999.

Lesaffre, B., Broska, J.-B., Coléou, C., Flin, F., and Pieritz, R.: Im-ages tridimensionnelles de neige: des prélèvements in situ aux fichiers de données volumiques. Application à une expérience de métamorphose d'isothermie, Tech. Rep., CNRM-GAME/CEN, Météo-France - CNRS, Grenoble, France, 2004.

Löwe, H., Spiegel, J., and Schneebeli, M.: Interfacial and structural relaxations of snow under isothermal conditions, J. Glaciol., 57, 499–510, doi:10.3189/002214311796905569, 2011.

Löwe, H., Riche, F., and Schneebeli, M.: A general treatment of snow microstructure exemplified by an improved rela-tion for thermal conductivity, The Cryosphere, 7, 1473–1480, doi:10.5194/tc-7-1473-2013, 2013.

Lundy, C. C., Edens, M. Q., and Brown, R. L.: Measurement of snow density and microstructure using computed tomography, J. Glaciol., 48, 312–316, doi:10.3189/172756502781831485, 2002.

Marbouty, D.: An experimental study of temperature-gradient metamorphism, J. Glaciol., 26, 303–312, 1980.

Markov, I. V.: Crystal growth for beginners, World Scientific Pub-lishing Co. Pte. Ltd, Singapore, 1995.

Mutaftschiev, B.: The atomistic nature of crystal growth, Series in Material Science, Springer, Berlin, 2001.

Nishikawa, Y., Koga, T., Hashimoto, T., and Jinnai, H.: Measure-ments of interfacial curvatures of bicontinuous structure from three-dimensional digital images, 2. A sectioning and fitting method, Langmuir, 17, 3254–3265, 2001.

Ogawa, N., Flin, F., and Brzoska, J.-B.: Representation of two cur-vatures of surface and its application to snow physics, Memoirs of the Hokkaido Institute of Technology, 34, 81–87, 2006.

Pinzer, B. and Schneebeli, M.: Breeding snow: an instrumented sample holder for simultaneous tomographic and thermal studies, Meas. Sci. Technol., 20, 095705, doi:10.1088/0957-0233/20/9/095705, 2009a.

Pinzer, B. and Schneebeli, M.: Snow metamorphism under al-ternating temperature gradients: Morphology and recrystal-lization in surface snow, Geophys. Res. Lett., 36, L23503, doi:10.1029/2009GL039618, 2009b.

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, doi:10.5194/tc-6-1141-2012, 2012.

Pottmann, H., Wallner, J., Huang, Q.-X., and Yang, Y.-L.: Integral invariants for robust geometry processing, Comput. Aided Geom. D., 26, 37–60, 2009.

Riche, F. and Schneebeli, M.: Thermal conductivity of snow mea-sured by three independent methods and anisotropy considera-tions, The Cryosphere, 7, 217–227, doi:10.5194/tc-7-217-2013, 2013.

Rieger, B., Timmermans, F. J., and van Vliet, L. J.: Estimation of curvature on surfaces in 3D grey-value images, in: Proc ASCI 2002, 8th annual conf. of the advanced school for com-puting and imaging, Delft, the Netherlands, 170–177, 2002.

Roberts, A. P., and Garboczi, E. J.: Computation of the linear elas-tic properties of random porous materials with a wide variety of microstructure, P. Roy. Soc. Lond., 458, 1033–1054, 2002.

Satyawali, P. K., Singh, A. K., Dewali, S. K., Kumar, P., and Ku-mar, V.: Time dependence of snow microstructure and associated effective thermal conductivity, Ann. Glaciol., 49, 43–50, 2008.

Schneebeli, M.: Three-dimensional snow: what snow really looks like, in: Proceedings of the International Snow Science Work-shop, 1–6 October 2000, Big Sky, Montana, edited by: Montana, S., American Avalanche Association, Big Sky, Montana, USA, 407–408, 2000.

Schneebeli, M. and Sokratov, S. A.: Tomography of tem-perature gradient metamorphism of snow and associated changes in heat conductivity, Hydrol. Process., 18, 3655–3665, doi:10.1002/hyp.5800, 2004.

Schneebeli, M., Pielmeier, C., and Johnson, J. B.: Measuring snow microstructure and hardness using a high resolution penetrometer, Cold Reg. Sci Technol., 30, 101–114, 1999.

Sethian, J. A.: Level set methods and fast marching methods: Evolv-ing interfaces in computational geometry, fluid mechanics, com-puter vision, and materials science, Vol. 3, Cambridge University Press, 1999.

Shertzer, R. H. and Adams, E. E.: Anisotropic thermal conductiv-ity model for dry snow, Cold Reg. Sci. Technol., 69, 122–128, doi:10.1016/j.coldregions.2011.09.005, 2011.

Shimizu, H.: Air permeability of deposited snow, Institute of Low Temperature Science, Hokkaido University, Contrib. Inst. Low Temp. Sci., A22, 1–32, 1970.

Srivastava, P., Mahajan, P., Satyawali, P., and Kumar, V.: Obser-vation of temperature gradient metamorphism in snow by X-ray computed microtomography: measurement of microstruc-ture parameters and simulation of linear elastic properties, Ann. Glaciol., 51, 73–82, 2010.

Staron, P. J., Adams, E. E., and Miller, D. A.: Nonequilibrium thermodynamics of kinetic metamorphism in snow, Cold Reg. Sci. Technol., 97, 60–71, doi:10.1016/j.coldregions.2013.10.007, 2014.

Sturm, M., Holmgren, J., König, M., and Morris, K.: The thermal conductivity of seasonal snow, J. Glaciol., 43, 26–41, 1997.

Thoemen, H., Walther, T., and Wiegmann, A.: 3D simulation of macroscopic heat and mass transfer properties from the mi-crostructure of wood fibre networks, Compos. Sci. Technol., 68, 608–616, doi:10.1016/j.compscitech.2007.10.014, 2008.

Torquato, S.: Random heterogeneous media: microstructure and improved bounds on the effective properties, Appl. Mech. Rev., 44, 37–76, 1991.

Torquato, S.: Effective stiffness tensor of composite media: ii. application to isotropic dispersions, J. Mech. Phys. Solids, 46, 1411–1440, 1998.

Torquato, S.: Random heterogeneous materials: microstructure and macroscopic properties, Springer, New York, USA, 2002.
Wang, X., Gillibert, L., Flin, F., and Coeurjolly, D.: Curvature-driven volumetric segmentation of binary shapes: An application to snow microstructure analysis, in: Proceedings of the 21st International Conference on Pattern Recognition, ICPR 2012, 11–15 November 2012, Tsukuba, Japan, 742–745, 2012.

Willis, J. R.: Bounds and self-consistent estimates for the overall properties of anisotropic composites, J. Mech. Phys. Solids, 25, 185–202, 1977.

Yen, Y. C.: Review of thermal properties of snow, ice and sea ice, Tech. Rep. 81-10, CRREL, Hanover, NH, USA, 1981.

Yosida, Z., Oura, H., Kuroiwa, D., Huzioka, T., Kojima, K., Aoki, S., and Kinosita, S.: Physical studies on deposited snow: I Thermal properties, Tech. Rep. 7, Institute of Low Temperature Science, Hokkaido University, Sapporo, Japan, 1955.

Zermatten, E., Haussener, S., Schneebeli, M., and Steinfeld, A.: Tomography-based determination of permeability and Dupuit-Forchheimer coefficient of characteristic snow samples, J. Glaciol., 57, 811–816, doi:10.3189/002214311798043799, 2011.

Zhang, Y., Paik, J., Koschan, A., Abidi, M. A., and Gorsich, D.: Simple and efficient algorithm for part decomposition of 3-D tri-angulated models based on curvature analysis, in: vol. 3, Proceedings of the 2002 International Conference on Image Processing, Rochester, New York, USA, III-273–III-276, 2002.