Skip to main content Accessibility help


  • Access
  • Cited by 43


      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

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

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

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

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Three-dimensional geometric measurements of snow microstructural evolution under isothermal conditions
        Available formats
Export citation


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 ND be the number of voxels p in D; we have:


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 SD 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 SD 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:


where nx(p), ny(p) and nz(p) are the projections of along the three coordinate directions of the voxel grid.

SD 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:


where ρ is the ice density and l0 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, ADx), ADy) and ADz) 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:


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πR2, we have:

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


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

Note that AD(α) 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:


where the values 1/r1(p) and 1/r2(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 d5-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):


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

The curvature distribution HD(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)).


This law is in fairly good agreement with the results based on CH4 adsorption at liquid nitrogen temperature (Cabanes and others, 2003; Legagneux and others, 2003).

4.3. Anisotropy

The angular distributions ADx), ADy) and ADz) 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/Cmean, Cmean being the mean curvature estimated from the eight considered images.

For a coordinate angle of about π/2, ADz) is significantly smaller than the distribution of the equivalent sphere (and significantly smaller than ADx) and ADy)). 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 HD(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 HD(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.


We thank X. Thibault and the team of line ID19 at the ESRF for the valuable experimental support.


Adams, E. E., Miller, D.A. and Brown, R. L.. 2001. Grain boundary ridge on sintered bonds between ice crystals. J. Appl. Phys., 90(11), 5782–5785.
Adamson, A.W. 1990. Physical chemistry of surfaces. Fifth edition. NewYork, etc, John Wiley and Sons.
Bartelt, P. and Lehning, M.. 2002. A physical SNOWPACK model for the Swiss avalanche warning. Part I. Numerical model. Cold Reg. Sci. Technol., 35(3), 123–145.
Bernache-Assolant, D. 1993. Chimie-physique du frittage. Paris, Edition Hermés.
Borgefors, G. 1986. Distance transformations in digital images. Computer Vision, Graphics and Image Processing, 34(3), 344–371.
Brown, R. L., Edens, M.Q. and Sato, A.. 1994. Metamorphism of fine-grained snow due to surface curvature differences. Ann. Glaciol., 19, 69–76.
Brun, E., David, P., Sudul, M. and Brunot, G.. 1992. A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting. J. Glaciol., 38(128), 13–22.
Brzoska, J.-B., Lesaffre, B., Coléou, C., Xu, K. and Pieritz, R. A.. 1999a. Computation of 3-D curvature on a wet snow sample. Eur. Phys. J. Appl. Phys., 7(1), 45–57.
Brzoska, J.-B. and 7 others. 1999b. 3-D visualization of snow samples by microtomography at low temperature. ESRFNewsletter, 32, 22–23.
Brzoska, J.-B. and 7 others. 2001. Computation of the surface area of natural snow 3-D images from X-ray tomography: two approaches. Image Anal. Stereol., 20(2), Suppl. 1, 306–312.
Bullard, J.W. 1997a. Digital-image-based models of two-dimensional microstructural evolution by surface diffusion and vapour transport. J. Appl. Phys., 81(1), 159–168.
Bullard, J.W. 1997b. Numerical simulations of transient-stage ostwald ripening and coalescence in two dimensions. Mater. Sci. Eng., Ser. A, 238(2), 128–139.
Bullard, J.W., Garboczi, E. J., Carter, W.C. and Fuller, E. R. Jr., 1995. Numerical methods for computing interfacial mean curvature. Comput. Mat. Sci., 4(2), 103–116.
Cabanes, A., Legagneux, L. and Dominé, F.. 2003. Rate of evolution of the specific surface area of surface snow layers. Environ. Sci. Technol., 37(4), 661–666.
Chermant, J.L. 1992. Caractérisation des poudres et des céramiques. Paris, Edition Herme's.
Coeurjolly, D., Flin, F., Teytaud, O. and Tougne, L.. 2003. Multigrid convergence and surface area estimation. In Asano, T., R. Klette and Ch. Ronse, eds. 11th International Workshop on Theoretical foundations of computer vision: Geometry, Morphology and Computational Imaging. Berlin, Springer-Verlag 101–119. (Lecture Notes in Computer Science 2626.)
Colbeck, S. C. 1997. A review of sintering in seasonal snow. CRREL Rep. 97-10.
Colbeck, S. C. 1998. Sintering in a dry snow cover. J. Appl. Phys., 84(8), 4585–4589.
Colbeck, S. C. 2001. Sintering of unequal grains. J. Appl. Phys., 89(8), 4612–4618.
Coléou, C., Lesaffre, B., Brzoska, J.-B., Ludwig, W. and Boller, E.. 2001. Three-dimensional snow images by X-ray microtomography. Ann. Glaciol., 32, 75–81.
Edens, M.Q. and Brown, R. L.. 1991. Changes in microstructure of snow under large deformations. J. Glaciol., 37(126), 193–202.
Flin, F., Brzoska, J.-B., Lesaffre, B., Coléou, C. and Lamboley, P.. 2001. Computation of normal vectors of discrete 3-Dobjects: application to natural snow images from X-ray tomography. ImageAnal. Stereol., 20(3), 187–191.
Flin, F., Brzoska, J. B., Lesaffre, B., Coléou, C. and R. A., Pieritz. 2003. Full three-dimensional modelling of curvature-dependent snow metamorphism: first results and comparison with experimental tomographic data. J. Phys. D, 36(10A), A49–A54.
German, R. M. 1996. Sintering theory and practice. NewYork, etc., John Wiley & Sons, Inc.
Good, W. 1987. Thin sections, serial cuts and 3-D analysis of snow. International Association of Hydrological Sciences Publication 162 (Symposium at Davos 1986−Avalanche Formation, Movement and Effects), 35–48.
Hirata, T. 1996. A unified linear-time algorithm for computing distance maps. Inf. Process. Lett., 58(3), 129–133.
Jordan, R. 1991. A one-dimensional temperature model for a snow cover: technical documentation for SNTHERM.89. CRREL Spec. Rep. 91116.
Legagneux, L., Lauzier, T., Dominé, F., Kuhs, W. F., Heindrichs, T. and Techmer, K.. 2003. Rate of decay of specific surface area of snow during isothermal experiments and morphological changes studied by scanning electron microscope. Can. J. Phys., 81(1/2), 459–468.
Lesaffre, B., Pougatch, E. and Martin, E.. 1998. Objective determination of snow-grain characteristics from images. Ann. Glaciol., 26, 112–118.
Meijster, A., Roerdink, J. B.T.M. and Hesselink, W.H.. 2000. A general algorithm for computing distance transforms in linear time. In Goutsias, J., L. Vincent and D. S. Bloomberg, eds. Mathematical morphology and its applications to image and signal processing. Dordrecht, Kluwer, 331–340.
Prewitt, B. 1970. Object enhancement and extraction. In Lipkin, B. and A. Rosenfeld, eds. Picture processing and psychopictories. New York, Academic Press, 75–149.
Rieger, B., Timmermans, F. J. and vanVliet, L. J.. 2002. Estimation of curvature on surfaces in 3-D grey-value images. In Deprettere, E. F., A. Belloum, J.W. J. Heijnsdijk and F. Van der Stappen, eds. Proceedings 8th Annual Conference of the Advanced School for Computing and Imaging, June 19–21 2001, Lochem, The Netherlands. Delft, Advanced School for Computing and Imaging, 170–177.
Serra, J. 1982. Image analysis and mathematical morphology. London, Academic Press.
Sethian, J.A. 1996. Level set methods, involving interfaces in geometry, fluid mechanics, computer vision, and materials sciences. Cambridge, etc., Cambridge University Press.
Verwer, B. J.H. 1991. Local distances for distance transformations in two and three dimensions. Pattern Recogn. Lett., 12(11), 671–682.