Three-dimensional geometric measurements of snow microstructural evolution under isothermal conditions

Abstract

Snow, from its fall until its full melting, undergoes a structural metamorphism that is governed by temperature and humidity fields. Among the many possible mechanisms that contribute to snow metamorphism, those that depend only on curvature are the most accessible to modelling. In this paper, techniques of volume data analysis adapted to the complex geometry of snow are introduced and then applied to experimental tomographic data coming from the isothermal metamorphism of snow near 0°C. In particular, an adaptive algorithm of curvature computation is described. Present results on the evolution of specific surface area and anisotropy already show that such image-analysis methods are relevant tools for the characterization of real snow microstructures. Moreover, the evolution of the curvature distribution with time provides valuable information for the development of sintering models, in the same way as a possible quantitative calibration of snow-grain coarsening laws.

1. Introduction

Modelling the behaviour of snow microstructure is essential for avalanche-risk forecasting, mainly as a means to provide a parameterization of grain-scale physics into existing or future models at the field scale (Jordan, 1991; Brun and others, 1992; Bartlet and Lehning, 2002). Existing models simulate the evolution of snow layers using an experimental parameterization ofmetamorphism. It is presently obtained from a characterization of snow in the field using grain silhouettes of isolated grains (Lesaffre and others, 1998) or two-dimensional (2-D) thin-section analysis (Good, 1987; Edens and Brown, 1991; Brown and others, 1994). Recent three-dimensional (3-D) imaging techniques (Coléou and others, 2001) are not yet implemented as a practical means of snow microstructure characterization. Substantial help could be found in the extensive research in the domain of materials and porous media (Bernache-Assolant, 1993; German, 1996; Bullard, 1997a). However, the various features of snow metamorphism (Colbeck, 1997) remain a particularity that needs to be addressed by the snow and ice community.

Isothermal metamorphism is probably the simplest situation, mainly governed by the distribution of mean curvature along the grain/pore interface, like many sintering processes. Near 0˚C, one can expect a reduced influence of possible surface or bulk diffusion at grain boundaries (Colbeck, 1999, 2001; Adams and others, 2001), owing to the large saturation vapour pressure of water at these temperatures. We recently proposed a simple evaporation scheme and performed preliminary tests using the first results of an experiment of controlled isothermal metamorphism at –2°C (Flin and others, 2003).

These tests do not presently allow any conclusions regarding the validity of the evaporation–condensation approach vs the more common (but much less tractable) diffusive scheme.

The characterization of microstructures is once again necessary in order to go farther in metamorphism modelling. Here we present some new developments of 3-D snow-image analysis and discuss their relevance to some typical behaviour observed in a real metamorphism experiment. After a short presentation of the experimental procedure, the numerical methods used in our work are described. The porosity, specific surface area, anisotropy and curvature distribution are then presented as a function of time and briefly discussed as possible guidelines for further modelling.

2. Experiment of Isothermal Metamorphism

A 3 month long experiment of isothermal metamorphism at –2°C was run at Col de Porte, Chartreuse mountain, French Alps, to provide microtomographic 3-D data for the validation of metamorphism models.

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 (Coléou and others, 2001) using a specially designed refrigerated cell (Brzoska 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 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. Geometrical Analysis

Presented here are the numerical algorithms we used to obtain the geometrical information describing our snow 3-D numerical images: porosity, specific surface area, anisotropy and curvature distribution.

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 (Serra, 1982; Chermant, 1992). In Brzoska 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 (Flin and others, 2001; Coeurjolly 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 Serra (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

The unit outward normal vector is first calculated on each surface voxel p 2 D, as in section 3.2.

The angles α_x(p), α_y(p) and α_z(p) are then computed.

The angular distributions on D, A_D(α_x), A_D(α_y) and A_D(α_z) are evaluated by counting the number of voxels corresponding to each class of angle.

On each histogram obtained, this number is expressed in percentage of the total surface voxels belonging to D and plotted in polar coordinates.

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 (Adamson, 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., Bullard and others, 1995; Rieger and others, 2002). In a preceding paper (Brzoska 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:

The determination, in each voxel p, of two orthogonal planes containing requires systematic rotations of a small discrete neighbourhood around p. This leads to some discretization effects and intensive computation time.

The use of a large fixed neighbourhood (about 20 voxels wide) for 2-D curve fitting implies that small objects cannot be processed properly.

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 (Borgefors, 1986; Verwer, 1991); other methods compute the exact Euclidean distance transform (Hirata, 1996; Meijster 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 (Verwer, 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 (Prewitt, 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 (Flin and others, 2001; Coeurjolly and others, 2003) and can be directly deduced from this computation.

C(p) is then given by the following relation (Sethian, 1996; Bullard, 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. Results and Discussion

Here are presented the evolution of porosity, specific surface area, anisotropy and curvature distribution during the experiment of isothermal dry-snow metamorphism. All computations were applied on 250 voxel (~2.5mm) wide images. Except for the SSA where a logarithmic scale was used, all times were counted from the first sampling.

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.

There is no significant discrepancy between the evolution of the three snow-layer thicknesses.

The precision on porosity obtained by small-volume numerical analysis seems to be equivalent to that obtained by a packing measurement.

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 Cabanes and others (2003) and Legagneux 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 (Cabanes and others, 2003; Legagneux 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.

a strong sharpening of the distribution, which denotes the rounding and smoothing of shapes

a decrease of the average mean curvature, which denotes the increase of the average grain-size

no significant evolution for negative curvatures

a wide and asymmetrical peak at the beginning of the metamorphism. This distribution is typical of fresh snow due to the presence of large planes (zero curvature) and acute shapes (high curvature) (see Fig. 5a)

an acute and symmetrical distribution at the end of the experiment, typical of rounded grains (see Fig. 5c).

5. Conclusion

Geometrical analysis tools have been presented and applied to an experiment monitoring the isothermal metamorphism of natural snow in three dimensions. The evolutions of porosity, SSA, anisotropy and curvature distributions have been obtained and discussed, showing a correct description of the well-known behaviour of an aging snow layer. A short-term application would be a quantitative calibration of the present empirical laws of grain coarsening and packing used in snow-cover models. Additionally, this documented dataset of the isothermal metamorphism constitutes a valuable tool for guiding basic and applied research in snow physics.

Presently, only a small part of the raw data from this experiment is actually exploited. To avoid the many theoretical problems due to sampling below the representative elementary volume, the presently developed techniques will have to be applied to the full-resolution dataset, especially for studies related to packing, heat transfer and micromechanics. Finally, these measurements provide quantitative data for the validation of metamorphism models and offer interesting prospects for improving models and understanding the physical phenomena that occur in snow microstructure.