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

microtomography. A large panel of properties was then computed from this series of 3-D images: density, specific surface area, correlation length, mean and Gaussian curvature distributions, air and ice tortuosities, effective thermal conductivity, and intrinsic permeability. Whenever possible, a specific attention was paid to assess these properties along the vertical and horizontal directions, and an anisotropy coefficient defined 10


Introduction
Natural snowpacks are frequently subjected to temperature gradients induced by their close environment. The morphology of the snow structure at the micro scale, or microstructure, evolves with time in a typical way which is called Temperature Gradient (TG) metamorphism and wherein one of the main features is the reorganization of 5 ice along the gradient direction by sublimation (resp. condensation) of ice (resp. water vapor) together with water vapor diffusion in the pore spaces (Yosida et al., 1955;Colbeck, 1983;Flin and Brzoska, 2008). In terms of snow type, this leads to facetted 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); Akitaya 10 (1974); Marbouty (1980); Colbeck (1983) and Satyawali et al. (2008) provide a good base of knowledge on the 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;Coléou et al., 2001;Lundy et al., 2002;Pinzer and Schneebeli, 2009;Chen and Baker, 2010), very precise studies related to 15 TG metamorphism are now available (Schneebeli and Sokratov, 2004;Kaempfer et al., 2005;Flin and Brzoska, 2008;Srivastava et al., 2010;Pinzer et al., 2012). They allow 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 20 as for the effective thermal conductivity (Yen, 1981;Sturm et al., 1997;Calonne et al., 2011) or the permeability (Shimizu, 1970;Jordan et al., 1999;Zermatten et al., 2011;Calonne et al., 2012), leading to simple parameterizations that can be used to estimate properties in snowpack models, e.g. Crocus (Brun et al., 1989) and Snowpack (Lehning et al., 1999). However, Marbouty (1980); 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/pores reorganization. Such studies suggest that there is a need to refine the 1409 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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;Löwe et al., 2013;Riche and Schneebeli, 2013) or the intrinsic permeability (Calonne et al., 2012), this type of snow exhibits specific anisotropic behaviors, and require more systematic investigations. 5 We propose to address 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 −3 subjected to a verti-10 cal temperature gradient of 43 K m −1 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 (three weeks of experiment). Moreover, these experimental conditions are in the range of conditions frequently undergone by natural alpine snowpacks. Snow cores were regularly sampled into the snow slab and scanned by X-ray most of these properties is provided, allowing to monitor the anisotropy of each quantity with time; and (iii) the physical properties computed on 3-D images are compared with those predicted by anisotropic analytical models based on basic microstructural properties.
2 Materials and methods 5

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 two weeks. This snow was then sieved in a cold room at −5 • C to obtain a horizontal snow slab of 100 cm×50 cm×14 cm composed of rounded grains (RG, Fierz et al., 2009) at 300 ± 15 kg m −3 (result from macroscopic density measure-10 ments). The snow slab was confined at the base and the top between two copper plates whose temperature was controlled by a thermo-regulated fluid circulation. The whole system was insulated by 8 cm thick polystyrene plates. An illustration of the experimental set-up is given in Fig. 1. Isothermal conditions at −5 • C were firstly applied to the snow slab during 24 h. This aimed at sintering snow grains whose bonds may have 15 been destroyed by sieving. During the following three weeks, the temperature of the cold room was held at −4 • C and the upper and lower copper plates were maintained at −1 • C and −7 • C, respectively, generating a steady vertical temperature gradient of 43 K m −1 through the snow slab.
The snow slab was sampled using a cylindrical core drill approximately every three 20 days over the three weeks, leading to seven samples in total at the end of the experiment. Macro 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 25 plate were punctually removed in order to access the snow slab. The samples were 1411 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 sampling, each snow sample was put in a plastic box and impregnated with 1-chloronaphthalene. This organic product, in liquid state 5 above −15 • C, was poured along the box walls, filling slowly the open porosities of snow. Then, the sample was frozen in an iso-octane bath cooled by dry ice (−78 • C) to allow the solidification of the 1-chloronaphthalene. The impregnation is required to stop the metamorphism of the snow microstructure and consolidate the snow sample for the further machining process. The absorption properties of the ice, air and 10 1-chloronaphthalene insure a correct contrast between these three components for the X-ray tomographic acquisition. Cylindrical snow cores were extracted from the samples by machining with a drill mounted on a lathe operating in a cold room at −30 • C. Each snow core was then stuck on the upper part of a copper sample holder by a droplet of 1-chloronaphthalene and sealed into a Plexiglas cap. The prepared samples were 15 finally stored at −20 • C until the tomographic acquisitions. Each core was scanned using the conical X-ray microtomograph of the 3S-R lab set with an acceleration voltage of 75 kV and a current intensity of 100 µA. As this microtomograph operates in an ambient temperature room, the snow core was placed in a specially designed cryogenic cell composed of a Peltier module, allowing to maintain 20 a regulated temperature of around −30 • C at the bottom part of the copper sample holder. A continuous dry and cold air circulation between the sample holder and the double-wall Plexiglas chambers of the cell prevents the condensation of water vapor on their sides. In addition, the heat generated by the Peltier module is dissipated by a water circulation. The whole system is able to rotate of 360 • during the acquisition. 25 Each tomographic acquisition lasted around 2 h during which 1200 radiographies were taken. Horizontal cross sections of the sample were reconstructed from radiographies using the DigiXCT 1 software. After image processing (see Flin et al., 2003;Lesaffre et al., 2004;Hagenmuller et al., 2013), we finally obtained seven 3-D images showing the microstructural evolution of the snow slab with time. The 3-D images have a resolution of 8.4 or 9.7 µm pixel −1 and a volume size of 5.9 3 , 9.2 3 or 9.7 3 mm 3 . Detailed information for each image is given in Table 1.

Density
Snow porosity φ (dimensionless), also called volume fraction of air, was estimated from 3-D images using a standard voxel counting algorithm. Snow density ρ s (in kg m −3 ) was simply deduced from φ such as ρ s = ρ i (1− φ) where ρ i is the ice density equal to 917 kg m −3 . 10

Specific surface area
The specific surface area in the x, y and z directions, noted SSA x , SSA y and SSA z (in m 2 kg −1 ), were estimated from 3-D images using a stereologic method (Arakawa et al., 2009;Flin et al., 2011) such as: 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, L is the total length of the testing lines (in m). We recall that the z direction corresponds to the direction of the gravity and of the macroscopic temperature gradient.

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 a: where r is a vector and the angular brackets denote the volume average. The two-point 10 probability function is also called the two-point correlation function or the autocorrelation function in the literature. For statistically homogeneous media, S 1 is simply equal to the porosity (φ) and S 2 depends on r. In general, S 2 has the following asymptotic properties (Torquato, 2002): Dashed lines in Fig. 2 show the two-point probability function of a snow volume after 500 h of metamorphism (referred as 7G in Table 1) computed along the x, y and z directions in blue, green and red, respectively. As proposed by Löwe et al. (2011Löwe et al. ( , 2013, 20 by fitting the S 2 (r) function along the coordinate axes β = (x, y, z) to an exponential Fig. 2), one obtains the correlation lengths l c β (in µm) in the x, y, and z directions noted l c x , l c y and l c z (called the vertical component), respectively. The average value of l c x and l c y is called the horizontal 1414 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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.

Mean and Gaussian curvatures
At a given point of a 3-D object, the surface is characterized by two principal curvatures 5 κ 1 and κ 2 , which correspond to the maximum and minimum values of curvature at this point. Negative, zero and positive values of κ 1 and κ 2 define concave, flat and convex surfaces of ice, respectively. The mean curvature C (in m −1 ) and Gaussian curvature G (in m −2 ) are defined as: 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 15 positive values represent saddle-shaped surfaces corresponding for snow to necks and connections between ice grains, flat or cylindrical surfaces, and dome-shaped surfaces (convex or concave), respectively. Note that the curvature value is inversely proportional to the typical size of the considered structure. Many techniques have been proposed to estimate mean and Gaussian curvatures 20 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 estimates obtained on a discrete grid. In our 25 approach, the mean (C) and Gaussian (G) curvatures are adaptively computed from 1415 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the largest relevant neighborhoods, limiting thus their digitization noise (see Flin et al., 2004;Brzoska et al., 2007;Wang et al., 2012 for details). In short, we use a variational approach for C and G (Sethian, 1999) where the mean curvature can be defined as the divergence of the normal vector field n(p) at point p: 5 and the Gaussian curvature, as in the following partial differential equation: 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 C (Eq. 8), the normal vector field n(p) could also have been expressed as a partial derivative on d . However, Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 C while decreasing the sensitivity of this formula to digitization effects (see Flin et al., 2004). For G, we simply use Eqs. (9)-(11) where derivatives are computed on local neighborhoods whose size is obtained by the adaptive analysis computed for 5 C (Wang et al., 2012).

Computations of tortuosity, effective thermal conductivity and permeability tensors
The full 3-D tensors of tortuosity τ (dimensionless), of effective thermal conductivity k (in W m −1 K −1 ) and of intrinsic permeability K (in m 2 ) were computed from 3-D images.

10
For that purpose, specific boundary value problems arising from the homogenization process (Auriault et al., 2009) have been numerically solved on Representative Elementary Volumes (REVs) extracted from 3-D images of snow by using the software Geodict 2 , based on a finite volume method (Thoemen et al., 2008). We define the REV of size l by Ω wherein Ω i and Ω a are the domains occupied by the ice and the air, 15 respectively, and where Γ denotes the common boundary.
To compute the effective thermal conductivity tensor k, the following boundary value problem was solved (Calonne et al., 2011): where I is the identity tensor, k a = 0.024 W m −1 K −1 and k i = 2.107 W m −1 K −1 stand for the thermal conductivity of air and ice at 271 K respectively, and the two periodic vectors t i and t a are unknown. The effective thermal conductivity tensor k is defined as: The tortuosity tensor of the air phase τ a (resp. of the ice phase τ i ) is obtained by solving the same boundary value problem assuming that k a = 1 and k i = 0 (resp. k a = 0 and k i = 1). These tensors are defined as: Let us remark that the air tortuosity is simply linked to the effective diffusion tensor 10 for the water vapor by D = φD m τ a , where D m is the molecular diffusion coefficient (in m 2 s −1 ) at the pore scale. If we assume that the porous medium is constituted of an equivalent tortuous capillary of total length l , in contrast with the REV length l , it can be shown that, by definition, 0 < τ a ∝ (l /l ) 2 < 1. Consequently, τ a (or τ i ) tends toward 0 or 1 when the air (or ice) structure is strongly tortuous or straight, respectively. The tortuosity is also often defined as τ f ∝ (l /l ) (Kaempfer et al., 2005) such as our tortuosity definition corresponds to the inverse of τ 2 f . The tensor K of intrinsic permeability was obtained by solving the following boundary value problem: where the fluid velocity v and the pressure fluctuationp (with p = 0) are the periodic unknowns for a given macroscopic gradient of pressure ∇p and where µ is the dynamic viscosity of air (in Pa s). It can be shown that v = −(1/µ)b∇p where b is a second order tensor (Auriault et al., 2009) and consequently, the permeability tensor is defined as: As the non-diagonal terms of the tensors τ i , τ 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 note z , xy and as the vertical component, the average of the two horizontal components 10 ( x and y ), and the average of the three components of any tensor = (τ i , τ a , k or K) respectively. Moreover, for the sake of simplicity, xy and are called the horizontal and the average components, respectively.

Computations of anisotropy coefficients
The anisotropy coefficient A( ) was computed for each of the microstructural and 15 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( ) = z / xy where = (1/SSA, l c , τ i , τ a , k or K). Thus, when the coefficient A( ) is close to 1, it is considered as isotropic, and when A( ) > 1 (resp. < 1), the property is anisotropic and larger along the vertical (resp. horizontal) direc- 20 tion. Note that, since SSA xy and SSA z characterize the vertical (resp. horizontal) surfaces, we study the anisotropy coefficient of the inverse of the specific surface area 1/SSA = (1/SSA x , 1/SSA y , 1/SSA z ) to be consistent with the other coefficients.

1419
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Re-adjustment in density
After sieving, the density of the snow slab exhibited slight spatial inhomogeneities (300 ± 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 5 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) and Calonne et al. (2012), respectively, as follows:

Analytical models based on ellipsoidal inclusions
We consider two analytical models based on ellipsoidal inclusions allowing to estimate the physical properties of snow including their anisotropy: the self consistent estimates and the dilute beds of ellipsoids. These models require basic microstructural information obtained from 3-D images, which are the volume fraction of each phase (φ, 1 − φ) 20 and the inclusion characteristics (size and anisotropy) given by the correlation lengths 1420 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2.6.1 Self consistent estimates: effective thermal conductivity and air tortuosity According to the self consistent method (Bruggeman, 1935;Hill, 1965;Budiansky, 1965;Willis, 1977;Torquato, 2002), the snow microstructure is considered as a macroscopically anisotropic composite of two different types of aligned isotropic ellipsoidal inclusions of the same shape, such as the semiaxes are defined as a = (l c x +l c y )/4 and 5 b = l c z /2, with volume fractions φ, 1−φ and conductivities k a , k i (see Fig. 3). Each type of inclusion is embedded in an homogeneous equivalent medium, i.e. an infinite matrix whose effective thermal conductivity k sc is the unknown to be calculated. The solution of equations for an isolated inclusion then gives an implicit relation which can be solved for this effective property (Eq. 25). For simple inclusion shapes, this method often yields 10 algebraic expressions of the effective properties. This self consistent scheme treats each phase symmetrically: Eq. (25) is invariant under the simultaneous interchange k a ↔ k i and φ ↔ 1 − φ. This particular property can be viewed as a way to capture the connectivity of both phases. By contrast, in the so called generalized self consistent method or three phase model (Hashin, 1962;Christensen and Lo, 1979;Boutin, 15 2000; Boutin and Geindreau, 2010), bi-composite inclusions are considered, typically a spherical particle (phase 1) surrounded by a concentric shell (phase 2), which allows to capture that the matrix (phase 2) and dispersed particles (phase 1) are connected and disconnected, respectively. In the present case, the standard self consistent estimate of the effective thermal 20 conductivity k sc verifies (Torquato, 2002): (25) where A is the depolarization tensor for an ellipsoid in a matrix with an effective thermal conductivity k sc and is defined in the (x, y, z) frame as: Let us remark that, for spheres, b/a = 1 and Q = 1/3; thus 5 the tensor A = I/3 and the tensor k sc is isotropic. In general, the inclusion shape implies that the tensor k sc is transverse isotropic. Thus, from Eqs. (25) and (26), the horizontal components of k sc are written: The self consistent estimate of the air tortuosity tensor τ sc a can be easily deduced from (29), and we have:  (Willis, 1977) are also reported. The corresponding microstructure of the lower bounds (resp. upper bounds) can be viewed as ellipsoidal inclusions of ice (resp. air) of same aspect ratio (b/a) dispersed within the air matrix (resp. ice matrix), as shown by the upper part of Fig. 3. As expected, in each direction, the self consistent estimate 20 lies between the bounds. At low volume fractions of ice (resp. at high volume fractions),

1422
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the self consistent estimates and the lower bounds (resp. the upper bounds) are very close. Finally, Fig. 3 clearly shows the anisotropy of the effective thermal conductivity induced by the anisotropy of the microstructure. In this particular case, A(k sc ) ranges from 1 to 2.8 in the whole range of ice volume fraction. 5 In that case, the snow is seen as a dilute dispersion of ellipsoids of ice in a matrix of air. As previously, the semiaxes of each ellipsoid are still 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 K el is written in the (x, y, z) frame:

Dilute beds of spheroids: permeability
15 and with 1423 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 estimation of the permeability tensor K di is proposed: This relation allows to recover a classical expression of the permeability as where h is a function of the porosity and d c is a characteristic length of the microstructure. Figure 4 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 bet-15 ter visualization of the facetted shapes, the images are presented "upside down": the top (respectively, base) of the images corresponds to the lowest (highest) and warmest (coolest) side of the physical sample. We observe qualitatively that the initial rounded grains become bigger and bigger and more angular and facetted with time. After about 200 h, depth hoar is obtained showing characteristic striations on the surface of grains 20 (see photographs in Fig. 4). At the end of the experiment, the ice structure is preferably arranged along the vertical direction, i.e. the direction of the temperature gradient, as shown by the cross-sections in particular.

Time evolution of microstructural and physical properties of snow
All the snow properties computed 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. 5. The blue, green and red symbols represent the x, y and z values of the considered property, respectively. One can observe that: -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 be-5 tween 275 and 315 kg m −3 are observed from one image to another (Fig. 5a). 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 values in the three directions of specific surface area decrease con-10 tinuously with time from 27.7 to 13.4 m 2 kg −1 . Between 0 and 144 h, the z values are slightly higher than the horizontal ones (29.2 against 26.9 m 2 kg −1 at 0 h). After 144 h, values in the three directions become very close from each other (Fig. 5b).
-The values of correlation length increase continuously during the experiment, evolving from 71 to 181 µm in the xy directions and from 68 to 228 µm in the 15 z direction (Fig. 5c).
-The values of air tortuosity (∼ 0.7) are around five times higher than those of ice tortuosity (∼ 0.15) (Fig. 5d). The overall evolution of both properties is low and can be divided in two stages: the ice (resp. air) tortuosity decreases (resp. increases) between 0 and 73 h and then slightly increases (resp. decreases) until the end 20 of the experiment. 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. 5e, exhibit the same variations as the snow density, showing the strong relationship between these two variables. 25 The values range between 0.14 and 0.26 W m −1 K −1 , all directions combined. . 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 as much significant. , all directions combined. 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. 5f). Both computed and readjusted values increase over the experiment and vertical values become 15 higher than the horizontal ones after 217 h. Figure 6 shows the time evolution of the anisotropy coefficient for six snow properties. In 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 vari-20 able to another. The largest evolution is shown by the coefficient of ice tortuosity A(τ 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(τ 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 with the Fig. 5. 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 the inverse of the specific surface area are the smallest and evolve during the experiment from 0.97 to 1.13 and from 0.92 to 1.01, respectively.
The distributions, expressed in terms of occurrence ratio, of the mean curvature of the upward (left) and downward (right) surfaces of ice are presented in Fig. 7. They   5 give the proportion in % of the ice surface area that exhibits a mean curvature located in a particular range of values. The time evolution of these distributions is shown by the plots of different colors. Initially (red color), the distribution of upward and downward surfaces are similar, with a peak of mean curvature located around 6 mm −1 and an occurrence ratio of ∼ 3 %. With time, the peak of mean curvature moves toward lower 10 curvature values and sharpens, meaning that ice structures tend to become larger and of monodisperse size. 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 −1 (occurrence ratio of 4.2 %) for the upward ones, and at 0 mm −1 (occurrence ratio of 4.8 %) for the downward ones.

15
Using the same representation, Fig. 8 shows the Gaussian curvature distributions computed from the entire ice-air interface. All the distributions are centered at 0 mm −2 , but become tighter over time, with maximum occurrence ratios evolving from 4.5 % to 17.3 % between 0 and 500 h, implying that the proportion of large, flat or cylindrical structures of ice increases. The initial distribution (red color) stands out from the oth-

Comparisons with analytical models
5 Figure 10 shows the comparison between the time evolution of computed values (dashed lines) and values given by analytical models (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 models are described in the z and xy directions (the vertical and the horizontal plane). In order to quantitatively compare models with 10 experiments, we use in the following E , the mean of relative differences, defined as where is a component or the anisotropy coefficient of τ a , k and K. E is given with the associated standard deviation. For the vertical values of properties, E τ a z = 6.9±5.4 %, E k z = −13.4 ± 13.2 % and E K z = −9.1 ± 11.7 %, whereas for the horizontal ones, E τ a xy 15 = 7.7±4.7 %, E k xy = −19.4 ± 9.0 % and E K xy = −8.5 ± 14.1 %.
Using the same representation, the time evolution of anisotropy coefficients from computed values and values given by models of the three above properties is shown in Fig. 11. For the air tortuosity, permeability and thermal conductivity, E A(τ a ) = −0.7 ± 0.9 %, E A(K) = −0.3 ± 2.9 % and E A(k) = 9.9 ± 27.9 %, respectively. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Evolution of the microstructure
From the mean curvature distributions (Fig. 7) and the 3-D images (Fig. 4), we observe that the snow microstructure undergoes a strong directional faceting under temperature gradient. At the beginning of the experiment, the top and base of ice grains are identigetically more stable for the condensation (resp. sublimation) of a water molecule (see TLK or Kossel crystal models in Markov, 1995;Mutaftschiev, 2001). Thus, the convex surfaces that are submitted to condensation undergo a layer by layer growth 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 15 of the shapes (see e.g. Knight, 1966;Flin and Brzoska, 2008, for more details). Our mean curvature results confirm the above considerations: during our metamorphism experiment, the top (resp. base) of the ice structures is warmer (resp. colder) than the surrounding air and sublimates (resp. condenses), inducing a rounding (resp. facetting) of the surface due to the crystalline properties of ice. 20 From the temporal evolution of the snow (e.g. Fig. 5), the metamorphism can be divided in two periods: (i) between 0 and 73 h, the microstructure, resulting from equitemperature 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 in favor to the growth of large structures (see Gaussian cur-25 vature distributions in Fig. 8). During this transitional state, the ice structure become more tortuous (minimum value of tortuosity over the whole experiment reached at 73 h in Fig. 5d), because of the decrease of the links between grains. (ii) After 73 h, the 1429 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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. Because the air tortuosity is about five times higher than the ice tortuosity (see Fig. 5d), our snow samples can be considered as constituted by a tortuous skele-5 ton 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 −3 a vertical ice tortuosity of τ 2 f = 4.4, which corresponds to 0.23 according to our tortuosity definition. In accordance with the observations of Schneebeli 10 and Sokratov (2004) and Satyawali et al. (2008), 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 re-built 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 coef-15 ficients 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).

Link with the physical properties
Heat conduction in snow is mostly due to the conduction of ice which conduces 100 20 times more than the air. Consequently, the effective thermal conductivity of snow is strongly linked to its density (ice proportion) and to the ice tortuosity (way for heat conduction), as illustrated in Fig. 5. 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 computations We observe two distinct stages in the evolution of the effective thermal conductivity (Fig. 5e): 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 (snow and temperature conditions were however different to that of our experiment). This behavior 5 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 temperature gradient limits the heat conduction in snow and (ii) the growth of large ice structures in a second time, in particular in the z direction, enhances the heat transfer.
As already mentioned, the intrinsic permeability is proportional to h(φ)τ a d 2 c , where 10 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., 2011). 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 15 that of the inverse of the specific surface area stands out by remaining close to one (see Fig. 6), our results seem to show that the anisotropy of the permeability is mainly related to that of the air tortuosity and the correlation length. 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 20 in Sect. 3.1, but undergoes the same continuous increase over the experiment that is observed for the characteristic length. Figure 5 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 verti-25 cal 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., 2011Calonne et al., , 1431 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2012) should be improved by introducing for instance new microstructural parameters which could reflect the anisotropy phenomena.

Modelling physical properties
The computed effective properties of snow were compared with analytical models based on basic information on microstructure such as the snow density and the 5 anisotropy of the grain shape through the correlation lengths. Even if these models fail to reflect all the details of the microstructure (bonds between grains ...), our results clearly show that they allow to capture in overall the evolution of the considered properties with a mean of relative differences ranging from −19.4 % to 7.7 % (see Sect. 3.1). These models do not include any fitting parameter, by contrast with the model of Löwe  Fig. 3) was adjusted to the computed data of thermal conductivity by introducing two parameters. These parameters are necessary since the lower bounds equation does not take into account the connectivity of both phases, in contrast with the self consistent estimates. In both of the presently proposed models, the time evolution of physical properties during the meta-15 morphism mainly depends on one parameter Q which is associated to the shape heterogeneity (pore/grain). They could be improved by introducing (i) other information on the microstructure by computing three-or four-point probability functions on 3-D images (Torquato, 2002) or (ii) refined variables allowing to describe the evolution of grain shapes and contacts between grains (mean and Gaussian curvatures, fabric ten-20 sor) during the metamorphism.

Conclusions
An experiment of temperature gradient (TG) metamorphism was realized in cold room in order to monitor the evolution of microstructural and physical properties of snow over time: seven snow samples were collected at regular intervals of time over three weeks into a snow slab subjected to a vertical temperature gradient of 43 K m −1 . Each sample was scanned by X-ray tomography to obtain a series of 3-D images showing the microstructural evolution of the snow slab during the temperature gradient metamorphism. Numerical computations were performed on these 3-D images to estimate various geometrical and physical properties of snow, such as the directional correla-5 tion 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)  These results highlight the strong interplay between microstructure and physical properties of snow, and confirm that the density alone, or any isotropic quantity, is not sufficient to describe the time evolution of physical properties during a TG metamorphism. To solve this problem, we applied analytical anisotropic models (self consistent estimates and dilute bed of spheroids) using microstructural parameters (directional 25 correlation lengths) that reflect the general shape of heterogeneities (size, anisotropy). The proposed analytical models, whose results were compared to those obtained by numerical computations, offer good estimates of the physical properties and anisotropy coefficients, without any fitting parameters.

1433
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | In overall, this study presents numerical tools to quantitatively 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-parameter to access directly with physical measurements during such a process. Used as a guideline or a validation tool, this database offers new 5 outlooks for the development of micro and macro-scale snow models.     24 Fig. 2. Two-point probability function S 2 (r ) for the whole range of r (left) and zoomed (right) in the x (blue), y (green) and z directions (red), obtained from the 3-D image referred 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).    For a better visualization of the facetted shapes, the images are presented "upside down". The arrows in blue, green and red correspond to the x-, y-and z-directions of images, respectively. 30 Fig. 9. The initial (left) and final (right) 3-D images of the experiment where colors represent the Gaussian curvature of the surfaces, ranging from −781 to +781 mm −2 . 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 facetted shapes, the images are presented "upside down". The arrows in blue, green and red correspond to the x, y and z directions of images, respectively.

1449
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig. 10. 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 models in the xy and z directions (solid lines in blue and red, respectively).

1450
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig. 11. 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 models (solid lines).