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

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 microor macroscale.


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 snowpack models, e.g., Crocus (Brun et al., 1989) and Snowpack (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 −3 subjected to a vertical 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 (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.

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 −3 (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 −1 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 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 sampling, each snow specimen was put in a plastic box and impregnated 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 allow for the solidification of the 1-chloronaphthalene.The impregnation 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 acquisition.Cylindrical snow cores were extracted from the samples 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 Plexiglas cap.The prepared samples were finally stored at −20 • C until the tomographic acquisitions.
Each core was scanned using the conical X-ray microtomograph of the 3SR lab set with an acceleration voltage of 75 kV and a current 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, 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 dissipated by water circulation.The whole system is able to rotate 360 • during the acquisition.Each tomographic acquisition 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 radiographs using DigiXCT1 software.Image processing was then applied to the grayscale, reconstructed images to obtain binary 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) removing 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 Hagenmuller et al. (2013) for detailed information on the exact procedure 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 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 Table 1.Microstructural and physical properties computed from 3-D images.Snow types are given according to the international classification (Fierz et al., 2009) , where ρ i is the ice density equal to 917 kg m −3 .

Specific surface area SSA
The specific surface area estimates along the x, y and z directions, denoted by SSA x , SSA y and SSA z (in m 2 kg −1 ), were computed from 3-D images, using a stereologic method (Arakawa et al., 2009;Flin et al., 2011): 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 SSA z is called the vertical component while the average value of SSA x and SSA y , noted SSA xy , 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.1. Solid and dashed lines correspond, respectively, to the values computed from Eq. ( 3) and to the results of the expression

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: where x is a position vector within the sample.The one-and two-point probability functions for the air phase are then defined as where 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 (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 (φ) and S 2 depends on r.In general, S 2 has the following asymptotic properties (Torquato, 2002): 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. (2011Löwe et al. ( , 2013)), by fitting the S 2 (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.

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 −1 ) and Gaussian curvature K (in m −2 ) are used to define the surface geometry with 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 www.the-cryosphere.net/8/2255/2014/  The Cryosphere, 8, 2255-2274, 2014 for H and K (Sethian, 1999), where the mean curvature can be defined as the divergence of the normal vector field n(p) at point p and the Gaussian curvature as with and 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 H (Eq. 8), the normal vector field 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 H while decreasing the sensitivity of this formula to digitization effects (see Flin et al., 2004).For K, we simply use Eqs.( 9)-( 11) where local estimations K(p) are averaged on the neighborhoods obtained for the adaptive analysis of H (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 m2 ) 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 2 , based on a finite difference method (Thoemen et al., 2008).We define the REV with a side length l by wherein i and a are the domains occupied by the ice and the air, respectively, and where denotes the common boundary.
To compute the effective thermal conductivity tensor k, the following boundary value problem was solved (Auriault et al., 2009;Calonne et al., 2011Calonne et al., , 2014)): where I is the identity tensor, n is the outward vector normal to the ice surface and the two periodic vectors t i and 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 ∇T applied on the REV.Finally, as in Calonne et al. (2011), stand for the thermal conductivity of air and ice at 271 K, respectively.The effective thermal conductivity tensor k is defined as The tortuosity tensor of the air phase τ a and of the ice phase τ i are obtained by solving the same above boundary value problem (Eqs.12-16) assuming that k a = 1 and k i = 0 for τ a , and k a = 0 and k i = 1 for τ i .These tensors are defined as Let us remark that for snow the air tortuosity is simply linked to the effective diffusion tensor for the water vapor by D = φ D m τ 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 < τ a ∝ (l/l ) 2 < 1 (Bear, 1972).Consequently, τ a tends toward 0 or 1 when the air structure is highly tortuous or straight, respectively (the same is true for τ i and the ice structure).The tortuosity is also often defined as τ f ∝ (l / l) (Kaempfer et al., 2005), so that 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 (Calonne et al., 2012): where v and p (with p = 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 ∇ p. µ 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 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 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 denote as z , xy and the vertical component, the average of the two horizontal components ( x and y ), and the average of the three components of any tensor (i.e., τ 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 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 = l SSA , l c , τ i , τ a , k or K.The property is considered isotropic if it exhibits a coefficient A( ) 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.

Readjustment 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 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. (2011Calonne et al. ( , 2012)), respectively, as follows: with ρ s as the computed snow density, ρ 294 = 294 kg m −3 , k fit (ρ s ) = 2.5 × 10 −6 ρ 2 s − 1.23 × 10 −4 ρ s + 0.024 and K fit (ρ s ) = 3.0 × r 2 es exp(−0.0130× ρ s ) where the equivalent sphere radius r es = 3/(SSA × ρ i ).In this way, we obtained readjusted thermal conductivity (k r x , k r y , k r z ) and permeability (K r x , K r y , K r z ) values for a density of 294 kg m −3 .

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 (φ, 1 − φ) 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 ).

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 (φ, 1 − φ) 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):  matrix with an effective thermal conductivity k sc .When k sc is transverse isotropic, the depolarization tensor A is defined in the (x, y, z) frame as (Giraud et al., 2007;Kushch and Sevostianov, 2014) with where γ is linked to the aspect ratio of the ellipsoid and anisotropy ratio of k sc by γ = (a/b) × (k sc xy /k sc z ).Thus, from Eqs. ( 25) and ( 26), the horizontal and vertical components of k sc are written as follows: where The self-consistent estimate of the air tortuosity tensor τ sc a can be easily deduced from Eq. ( 25) with k a = 1 and 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 ↔ k i and φ ↔ 1 − φ, meaning that each phase is treated symmetrically.
As an illustration, Fig. 4 shows the behavior of k sc x = k sc y and k sc z 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(k sc ) is defined as k sc x /k sc z and consequently, from Eqs. ( 29) and (30), is a function of the ratio k i /k a , the porosity φ, and the aspect ratio b/a.In the particular case of Fig. 4, k i /k a 100 and b/a = 1.45; so A(k sc ) depends only on the porosity and ranges from 1 to 1.6 in the whole range of the ice volume fraction.

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 K el is written in the (x, y, z) frame as where and where χ a and χ b are linked to the aspect ratio of the ellipsoid as χ 2 a = −χ 2 b = (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 K di is proposed such as where K el xy = K el x = K el y and K el z are the diagonal components of K 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) where h is a function of the porosity and d c is a characteristic length of the microstructure.

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.9m 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 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,   referred to as "readjusted", show values of thermal conductivity after readjustment at a density of 294 kg m −3 .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 −9 and 4.84 × 10 −9 m 2 .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(τ 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 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 −1 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 −1 (occurrence ratio of 4.2 %) for the upward ones, and at 0 mm −1 (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 −2 , but become sharper over time: the maximum  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 −2 and large between around −500 and −50 mm −2 , 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 −2 , 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.

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( ) 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.3 ± 4.5 %, E(k z ) = −22.6 ± 4.32 % and E(K z ) = −9.6 ± 10.9 %, whereas for the horizontal   ones, E(τ a xy ) = 8.1 ± 5.2 %, E(k xy ) = −15.6 ± 6.0 % and E(K xy ) = −8.1 ± 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(A(τ a )) = −1.6 ± 0.9 %, E(A(k)) = −7.9± 7.4 %, and E(A(K)) = −1.1 ± 4.2 %, respectively.

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 K m −1 at −3 • 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, andFlin andBrzoska, 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 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 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 −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 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).

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 twostage evolution of this property in roughly similar conditions (average temperature of −8 • C, temperature gradient of 100 K m −1 , snow density of 268 kg m −3 ).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(φ)τ a d 2 c , 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., 2011Calonne et al., , 2012) ) should be improved by introducing for instance new microstructural parameters which could reflect the anisotropy phenomena.

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 improvements 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 ellipsoid.These computations are not straightforward.Finally, other relevant estimates of these effective properties can be probably obtained by introducing (i) other microstructural information in bounds or exact contrast expansions by computing three-or four-point probability functions on 3-D images (Torquato, 1991(Torquato, , 1998(Torquato, , 2002;;Roberts and Garboczi, 2002) or (ii) refined variables allowing us to describe the evolution of grain shapes and contacts between grains (mean and Gaussian curvatures, fabric tensor) during the metamorphism.

Conclusions
An experiment of TG metamorphism was performed in a cold room in order to monitor the evolution of microstructural 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 −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 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 correlation 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 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 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 proposed analytical estimates, whose results were compared to those obtained by numerical computations, offer good estimations of the physical properties and anisotropy coefficients for our time series, without applying any fitting parameters.
In summary, 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 outlooks for the development of microand macroscale snow models.
The Supplement related to this article is available online at doi:10.5194/tc-8-2255-2014-supplement.

Figure 1 .
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.

Figure 2 .
Figure 2. Illustration of the cryogenic cell used during the tomographic acquisition.

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

Figure 4 .
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.

Figure 5 .
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 −1 .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.

Figure 6 .Figure 7 .
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.

Figure 8 .
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 −1 wide.

Figure 9 .
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 −2 wide.

Figure 10 .
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 −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 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.

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

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