2.1. Sampling

A 0.5 m × 1m slab of recent snow was first collected in the field (on16 January 2002) 15 hours after a snowfall (exterior temperature –1°C, slab thickness 12 cm). To prevent sublimation and temperature gradient effects, the slab was stored in a closed Styrodur^® (extruded rigid polystyrene foam) box inside the cold room maintained at –2±0.2°C. Till the end of the experiment, all manipulations were done in the cold room at this temperature. A 3 cm wide core was sampled at increasing time intervals, ranging from 24 hours at the beginning to 1 week at the end of the experiment. The collected slab was constituted of three layers whose thicknesses were measured at each sampling. All samples were taken in the middle layer of the slab, always >10 cm away from already sampled regions.

2.2. X-ray absorption microtomography

Once sampled, each core was impregnated with 1-chloronaphthalene and machined into the shape of a 9 × 9 mm cylinder for microtomographic acquisition. Three-dimensional images of snow samples were obtained at the ID19 beamline of the European Synchrotron Radiation Facility (ESRF) at Grenoble, France, by X-ray absorption microtomography (Reference Coléou, Lesaffre, Brzoska, Ludwig and BollerColéou and others, 2001) using a specially designed refrigerated cell (Reference BrzoskaBrzoska and others, 1999b). All the images were obtained at 18–20 keV, with a voxel (volume element) size of 4.91 μm. The grey-level images obtained, reconstructed at the ESRF, were contoured using a semi-automatic procedure. More information about core preparation and image processing can be found in Reference Flin, Brzoska, Lesaffre, Coléou and R. A.Flin and others (2003). The 3-D images obtained were 600 voxels wide, which is a considerable quantity of volumetric data. To allow reasonable processing times for our normal and curvature algorithms, the resolution of the 3-D images was reduced by a factor of 2 in the three axes. Thus, in all the 3-D snow images presented in this paper, one voxel corresponds to 9.82 μm.

3.1. Porosity

The porosity P is the ratio of the pore volume by the total volume of the considered sample. The evolution of this parameter during the metamorphism is directly related to the packing of the snow. It can be obtained simply by counting the voxels belonging to the pore phase in a chosen region of interestD. Let N_D be the number of voxels p in D; we have:

(1)

where f(p) = 1 for p belonging to the pore, 0 otherwise.

3.2. Specific surface area

The specific surface area (SSA) of a snow sample is defined by the total surface area of the air/ice interface per mass unity of the considered sample. The SSA is a useful parameter for describing the metamorphism of a snow layer, as it indicates the potential of this snow to undergo physical evolution. This parameter can be computed by estimating the surface area S_D included in a region of interest D. Stereological methods can be used to obtain such a surface estimation, and give good results for section planes that contain a large number of grains (Reference SerraSerra, 1982; Reference ChermantChermant, 1992). In Reference BrzoskaBrzoska and others (2001), an original method was proposed and compared to other non-stereological methods. Here the algorithm we used for the S_D determination is briefly reviewed.

The unit outward normal vectors are adaptively calculated on each surface voxel p 2 D by an original algorithm (Reference Flin, Brzoska, Lesaffre, Coléou and LamboleyFlin and others, 2001; Reference Coeurjolly, Flin, Teytaud and TougneCoeurjolly and others, 2003).

For each surface voxel p, a weight g(p) is computed as follows:

(2)

where n_x(p), n_y(p) and n_z(p) are the projections of along the three coordinate directions of the voxel grid.

S_D is obtained, in voxel units, by summing all the contributions g(p) in D.

Then, the SSA in physical units can be calculated as follows:

(3)

where ρ is the ice density and l₀ the size of one voxel.

3.3. Anisotropy

By “anisotropy”, we mean the anisotropy of the normal vector field describing the ice surface. An estimation of the anisotropy in snow samples is of great importance for understanding processes during metamorphism and packing. We propose here a tridimensional anisotropy estimator inspired by the work of Reference SerraSerra (1982) but using the knowledge of the normal vector in each surface voxel. Let us define α_x(p), α_y(p) and α_z(p) the three coordinate angles of

Note that the distributions are only defined between 0 and π rad. As this algorithm consists in interpreting volumetric data by angles with one direction of the space, the distributions obtained are not straightforward.

For example, let us consider the distributions of a sphere. As this object is fundamentally isotropic, the distributions are invariant by rotation and we have:

(4)

with α the angle between and a chosen direction of the space. The number of points corresponding to a class of angle α is equal to the length of the parallel line defined by α. This value amounts to 2πR sin α, R being the radius of the sphere (Fig. 1). The surface area of a sphere being 4πR², we have:

Fig. 1. For a sphere of radius R;A_D(α) = (100 sin α)/2R%. This distribution has a maximal value for π/2 (equatorial line).

(5)

with A_D(α) expressed in percentage of the total surface voxels belonging to D.

Note that A_D(α) will have a maximal value at π/2 (equatorial line) and minima at 0 and π rad (poles).

For more complex objects, one can detect anisotropies by comparing the distributions obtained in the three directions of the space and the distribution of its equivalent sphere (see section 4.3).

3.4. Mean curvature distribution

The 3-D mean curvature C of a surface corresponds to the quantity δA/(2δV ), where δA is the incremental change in the element’s area when it is normally displaced by local addition of material of volume δV. As can be seen from Kelvin’s equation (Reference AdamsonAdamson, 1990), the mean curvature is directly involved in vapour exchanges, and the determination of this quantity at each point of the surface is essential for monitoring the snow metamorphism.

Several methods were proposed for computing mean curvatures (see, e.g., Reference Bullard, Garboczi, Carter and FullerBullard and others, 1995; Reference Rieger, Timmermans and vanVlietRieger and others, 2002). In a preceding paper (Reference Brzoska, Lesaffre, Coléou, Xu and PieritzBrzoska and others, 1999a), we presented such an algorithm, in which C(p), the curvature on each voxel p, was computed as follows:

(6)

where the values 1/r₁(p) and 1/r₂(p) of the 2-D curvature were obtained on two orthogonal planes that contain by fitting the discrete curves by parabolas. This method gives fairly good results but presents some drawbacks:

In this paper, we propose an adaptive algorithm that avoids the above disadvantages. This method involves the following steps:

The background distance map of the ice phase M (_ice) is first generated. In other words, we label all voxels in the object with the distance to the closest background voxel. Many algorithms exist to compute such a map: some of them use chamfer metrics to approximate the Euclidean metric (Reference BorgeforsBorgefors, 1986; Reference VerwerVerwer, 1991); other methods compute the exact Euclidean distance transform (Reference HirataHirata, 1996; Reference Meijster, Roerdink and HesselinkMeijster and others, 2000). We chose a second-neighbour chamfer distance d_{5-7-9-11-12-15} due to its simplicity and its good accuracy (Reference VerwerVerwer, 1991).

The background distance map of the pore phase M(_pore) is then constructed. Note that the same nearest background points used for the determination of M(_ice) should be used for M(_pore).

A distance map on the whole image M(_image) is obtained by combining the two background distance maps as follows: M(_image) = M(_ice) -M(_pore).

We compute the gradient map of M(image) using a classical Prewitt first-neighbour mask derivative filter (Reference PrewittPrewitt, 1970).

In each surface voxel p, the divergence of the unit-normal vector is estimated by averaging the first-neighbour divergences of the background distance gradient in an appropriate surface neighbourhood around p. Note that this neighbourhood is obtained by applying the same angular and symmetry criteria as for the determination of normals (Reference Flin, Brzoska, Lesaffre, Coléou and LamboleyFlin and others, 2001; Reference Coeurjolly, Flin, Teytaud and TougneCoeurjolly and others, 2003) and can be directly deduced from this computation.

C(p) is then given by the following relation (Reference SethianSethian, 1996; Reference BullardBullard, 1997b):

(7)

where d is the dimensionality of the considered R_d (d = 3 in this paper).

The curvature distribution H_D(C) on D is evaluated by counting the number of voxels corresponding to each class of curvature. On each histogram obtained, this number is expressed in percentage of the total surface voxels belonging to D.

4.1. Porosity evolution

The porosity was computed for ten tomographic images and is plotted vs time in Figure 2a. To estimate the soundness of the results obtained, the inverse of density was computed from the porosity value. This quantity was compared to the thickness variation of the snow layers measured during the experiment (see Fig. 2b). The inverse of density and the thickness data are expressed as a percentage of the initial values (obtained 15 hours after the snowfall). Note that the thickness estimation of individual layers is inaccurate because of the difficulty of distinguishing the limit between two layers. The uncertainty is particularly significant for the middle layer, for which the errors come from the estimation of both the upper and lower borders. For these reasons, the slight differences of the packing rate that can be observed between the three layers are not significant. The total layer-thickness measurement is much more accurate, as it does not involve any difficult estimation of the border position. Despite evident noise due to measurement errors and small-volume analysis, these five curves are in fairly good agreement and seem to follow the classical dry-snow packing laws. This allows us to conclude that:

Fig. 2. Porosity and snow-layer evolution with time. Porosity evolution (a) and comparisons between thickness of snow layers and inverse of density evolutions (b). There is fairly good agreement between snow-layer measurement and numerical estimation at microstructural scale.

4.2. Specific surface area evolution

The SSA evolution with time t is plotted in Figure 3. By fitting this evolution with a logarithmic function, we obtain the following time dependency with a correlation coefficient of 0.986:

Fig. 3. SSA plotted in logarithmic scale. Time is counted from the begining of the snowfall. The SSA evolution follows a logarithmic law as mentioned in Reference Cabanes, Legagneux and DominéCabanes and others (2003) and Reference Legagneux, Lauzier, Dominé, Kuhs, Heindrichs and TechmerLegagneux and others (2003) (see Equation (8)).

(8)

This law is in fairly good agreement with the results based on CH₄ adsorption at liquid nitrogen temperature (Reference Cabanes, Legagneux and DominéCabanes and others, 2003; Reference Legagneux, Lauzier, Dominé, Kuhs, Heindrichs and TechmerLegagneux and others, 2003).

4.3. Anisotropy

The angular distributions A_D(α_x), A_D(α_y) and A_D(α_z) were estimated on each tomographic image. The early stages of the metamorphism show a strong angular variability of the distributions: this can be explained by the numerous large plates contained in the images and the relatively small region studied. After the first 66 hours of the experiment, the variability of the distributions is decreasing and one can observe a slight but persistent anisotropy along the vertical direction (z axis). This phenomenon is illustrated in Figure 4, in which the distributions were averaged on the eight last stages of the experiment (66–2011hours).

Fig. 4. Anisotropy visualization on polar diagram: a slight but persistent packing is visible along the vertical direction (z axis) for times >66 hours.

In this graph, the distribution of the equivalent sphere was plotted to help in comparing distributions along the different axes. For obtaining the equivalent sphere distribution, Equation (5) was applied with R = 1/C_mean, C_mean being the mean curvature estimated from the eight considered images.

For a coordinate angle of about π/2, A_D(α_z) is significantly smaller than the distribution of the equivalent sphere (and significantly smaller than A_D(α_x) and A_D(α_y)). This indicates that fewer normals are orthogonal to the z axis than are orthogonal to the other axes. This observation is consistent with the values of the distributions near 0 and π and indicates a slight but persistent anisotropy along the z axis.

This slight packing along the vertical direction suggests gravity effects during metamorphism. They should be taken into account to find out and simulate the mechanisms implied in isothermal metamorphism of snow.

4.4. Curvature evolution

Curvature computations were applied on ten microtomographic images. For each processed image, the mean curvature distribution H_D(C) was estimated and is plotted in Figure 5d. In this graph one can see:

Fig. 5. Curvature evolution during isothermal metamorphism of a dry snow sample. Grains are clearly growing and rounding during the metamorphism. (a–c) First stage (a), intermediate stage (b) and last stage (c) of the isothermal experiment. Image edges are 256 voxels (~2.5 mm) wide. (d) The evolution of curvature distribution H_D(C) during the metamorphism.