Skip to main content Accessibility help


  • Access
  • Cited by 8


      • 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.

        Explicit iterative computation of diffusive vapour field in the 3-D snow matrix: preliminary results for low flux metamorphism
        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.

        Explicit iterative computation of diffusive vapour field in the 3-D snow matrix: preliminary results for low flux metamorphism
        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.

        Explicit iterative computation of diffusive vapour field in the 3-D snow matrix: preliminary results for low flux metamorphism
        Available formats
Export citation


The metamorphism of seasonal snow is classically considered as limited by vapour diffusion in the pore phase. To account for the lack of knowledge of the ice–vapour reaction coefficient near 0 ° C, the assumption of a reaction-limited metamorphism was first tested in three-dimensional simulations at low and very low temperature gradients; however, the validity of such results is difficult to verify experimentally. By a reasoned use of traditional iterative schemes, vapour diffusion is now simulated in three dimensions on tomographic snow data, mapping the gradient of vapour pressure near the grains. Repeating this process may provide a way to simulate the isothermal metamorphism without grain packing at a reasonable expense of computation time. Preliminary results are compared with existing computations made within the reaction-limited hypothesis.

1. Introduction

The metamorphism of dry seasonal snow on the ground is classically described by water-vapour transfer between grains (Colbeck, 1983a). Depending on the degree of thermodynamic disequilibrium, the microstructure will evolve towards either faceted (Yosida and others, 1955) or rounded grain shapes (Colbeck, 1983a). Among other things, the formation of facets requires a sustained supersaturation with regard to vapour. In a snow cover, this usually happens where a temperature gradient (TG) exists. The formation of rounded grains, which occurs close to equilibrium (no TG, with a temperature T close to 0 ˚ C), is often called isothermal metamorphism of temperate snow. This is the subject of the present study.

The observed vapour transfer between metamorphosing snow grains results from some coupling of reaction and diffusion processes. As in many problems of kinetics, the macroscopic morphology will show the action of the slowest process: it is said that the problem is limited by this process. The main question is whether isothermal snow metamorphism is either reaction or diffusion-limited.

The diffusion-limited hypothesis is clearly the most popular in the snow science community (see, e.g. Colbeck, 1997; Legagneux and Dominé, 2005). In ceramic science, most materials of practical interest have a very low saturating vapour pressure p sat at the temperature they are processed. Conversely, the few vapour-grown industrial materials such as SiC are processed under high temperature, low total pressure and high supersaturation, (see, e.g. Chen and others, 2001). An overview of the SiC bulk growth technique is available at page/research/sic/Chapter3.html#3.1.

Literature on the regime of very low supersaturation for materials at a high partial vapour pressure p vap (typically, for a snow sample at 2 ˚ C, p vap ∼ p sat 517.3 Pa) is scarce. The reaction coefficient α ∈ [0, 1] is formally defined as the ratio between the number of reacting (adsorbing) and incident molecules on a surface; it is the standard way of addressing the question of reaction vs diffusion (Mutaftschiev, 2001). A low value of α implies the growth process is reaction-limited. A high value is related to a diffusion-limited growth, except for vacuum conditions where the mean free path of water vapour compares to the crystal size. Existing measurements of α on ice were often carried out at low temperature (T<−100 ˚ C) and vacuum conditions, forcing a reaction-limited growth, which in turn tells us how to relate the growth rate to α. Regardless, results range within more than two orders of magnitude (Brown and others, 1996).

As an exception, a recent work (Libbrecht, 2003) presents data near 0 ˚ C with low values of α for low supersaturations σ ≡ p vap /p sat 1. Although such conditions of very low growth rates may lead to higher errors in the determination of α for experimental reasons, this finding would advocate a reaction-limited metamorphism, contrasting with the common assumption (Colbeck, 1997).

For this reason, three-dimensional (3-D) simulations of snow metamorphism were carried out assuming reaction-limited conditions for isothermal (Flin and others, 2003) and TG (Flin and others, 2007) metamorphism. The requirement of dealing with extremely low vapour fluxes makes the design of a conclusive experiment very difficult, for example measuring the equilibrium shape of isolated crystals (Eaglesham and others, 1993). Here, the results of a simple model of diffusion-limited isothermal metamorphism are presented, then applied to the same 3-D data as the formerly developed reaction-limited model (Flin and others, 2003), in order to prepare further intercomparisons between the two approaches.

2. Modelling

2.1. Models for isothermal metamorphism

In low-flux metamorphism at high temperatures (close to 0 ˚ C) and high p sat, the following can be assumed, regardless of the limiting process:

  • 1. The mass transfer in the pore phase only occurs by molecular diffusion. There is no convection in snow, even for moderate thermal gradients (Brun and Touvier, 1987; Kaempfer and others, 2005).

  • 2. Heat conduction across the grain phase is the fastest process (Sekerka, 2004). The latent heat is quickly removed through the highly conductive ice (Colbeck, 1983b) and the temperature T is considered constant and uniform in the snow matrix.

  • 3. Grain metamorphism is the slowest process. For a given state of the microstructure, steady-state solutions are expected for both vapour diffusion (referred to as diffusion) and sublimation/condensation (referred to as reaction) processes.

  • 4. The grain’s interface stays microscopically rough everywhere, i.e. it behaves like a surface of type K (Mutaftschiev, 2001) saturated with kinks. Thus we do not expect faceted shapes in high-temperature isothermal metamorphism of snow.

As mentioned in the introduction, two different mechanisms can be assumed to be controlling the isothermal metamorphism: a reaction mechanism at the ice–pore interface and a diffusion mechanism through the pore space. After recalling the main features of the reaction-limited hypothesis, we will focus on the diffusion-limited assumption.

Reaction-limited process

In case of a reaction-limited process, the vapour field is assumed to be constant (steady state) and uniform (fast diffusion) at a value p amb(T) such that no macroscopic vapour flux is observed at the scale of the sample in respect to assumption (1). This happens when vapour is in equilibrium with any infinite surface whose mean curvature κ (n being the unit local normal, κ ≡ (div n)/2) equals the integral mean curvature 〈κ〉 of snow in the sample. The ice surface is considered out of equilibrium with respect to the surrounding uniform pressure field, and its growth is governed by the Langmuir–Knudsen equation that gives the mass vapour flux j vap:


where R and M are the ideal gas constant and molar mass of vapour, respectively. p sat(κ, T ) is given by Kelvin’s law:


where p sat (T ), γ and V spe are the saturating vapour pressure over a macroscopically flat, microscopically rough surface of ice, the ice–vapour surface tension and the specific volume of ice, respectively. If the curvature field is available on a sample, local growth rates can be computed in a straightforward manner. This made early simulations of metamorphism possible (Flin and others, 2003).

Diffusion-limited process

In the present paper, the metamorphism is assumed to be limited by a vapour diffusion process. In the pore space, p vap should obey the stationary equation of diffusion:


where r is the position vector. The grain surface S is considered to be in equilibrium with the vapour, at a uniform and constant temperature T . This provides a Dirichlet boundary condition to the vapour field, i.e.


After solving Equation (3) by using the boundary conditions of Equation (4), the mass flux of vapour jvap is given by Fick’s law:


Here, D is the diffusion coefficient of water vapour inm2 s 1, c is the concentration in number (mol mol 1) and ρ is density. Since both air and vapour are ideal gases, Fick’s law can be written:


The growth rate R is then directly related to the gradient of vapour pressure along the surface unit normal vector n:


In both situations, the mean curvature κ(r) of the grain’s surface governs the local flux of water vapour either directly in case of a reaction-limited process (Equation (1)) or as a boundary condition (Equation (4)) for a diffusion-limited process.

2.2. Continuous problem

The present process of vapour diffusion occurs in the physical continuous space. It is considered regardless of any digitization process, the first of which is tomographic image acquisition. This is basically a Dirichlet problem for Laplace’s equation whose boundary conditions are constant but not uniform, unlike most common problems of electrostatics (Durand, 1966b). Besides, it should be noted that the above boundary condition (Equation (4)) only concerns the physically relevant boundaries of the system (the grain’s interfaces), but not the practical spatial boundaries of the data file (the ‘edges’ of the cube). This point is addressed later.

We consider a problem where all boundary values are fixed (inner Dirichlet problem). Physically, this may be thought of as coating the sample edges with a continuous layer of bulk ice. The existence and uniqueness of a harmonic solution (which solves Laplace’s equation) to the Dirichlet problem (prescribed value of the parameter at the boundaries) is a classical result. However, the standard proof assumes a uniform value of the parameter over any connected surface (in electrostatics, the potential). This result was long ago extended (De la Vallée Poussin, 1932) to non-uniform values at the boundaries of a finite system, provided the following conditions are fulfilled.

  • 1. The domain is of finite size.

  • 2. The pore space D is connected (there may be more than one connected grain cluster in the domain).

  • 3. p vap is continuous at the surface S of any connected grain cluster.

  • 4. Any point denoted by its position vector r ∈ S fulfils the so-called Poincaré condition: there exists a neighbourhood N(r) and a right cone C whose apex is r and such that C ∩N(r) ∩ D = . Practically, this means that even dihedral angles (not accounted for in the present study) are allowed; only inward cusps are not.

It is then mathematically meaningful to solve the inner Dirichlet problem with boundary conditions that are non-uniform but continuous.

2.3. Numerical methods

Most standard methods for solving Laplace’s equation rely on digitizing the Laplacian operator on a regular grid, then solving the linear system resulting from various efficient implicit or semi-implicit schemes. Unfortunately, the present-day size of 3-D image data files, which often amounts to more than 10003 volume elements (voxels), prevents the direct use of techniques involving matrix inversion or Gauss pivoting (Press and others, 1992). Besides, the algorithmic complications of adaptive non-structured meshes, which allow a substantial reduction in size of the linear system to be solved (Frey and George, 2000), seem to be little adapted to the complex topology of 3-D data from sintered materials. For these reasons we used traditional explicit iterative (non-matrix) schemes based on the Gauss–Seidel method (Durand, 1966a; Patankar, 1980; Press and others, 1992). The present work uses the early Jacobi method, where the result at the nth iteration uses only neighbouring values from the (n−1)th iteration as a way to understand the assumptions made at the data file boundaries (see section 2.4).

The Laplacian operator Δp was digitized on all first neighbours of a voxel, following a method described by Durand (1966b). The 26 first neighbours of a given voxel can be sorted by their centre-to-centre distance d to the central voxel (denoted by index 0) comprising: 6 side voxels of d = 1, 12 edge voxels of and 8 corner voxels of By considering Equation (3), the third-order Taylor expansion of p yields


This is the common digitization used to account for Δ = xx + yy + zz = 0. The fifth-order expansion includes all first neighbours and leads to


Digitization Equation (9) was used everywhere in the pore space D except when it was not defined, that is on voxels close to the surface and whose neighbourhood contains at least one voxel in the ice phase. In the latter case, digitization Equation (8) was used. This approximation is mathematically consistent since both digitizations are Taylor expansions of the pressure field. To quantify the effect of the change of digitization scheme close to the ice–pore interface, Laplacian maps have been computed and provide reasonable errors.

2.4. Real boundary conditions

The aim of modelling metamorphism is to provide snapshots of the snow microstructure that are physically relevant to a snow layer at the macroscopic scale. The first requirement is that the studied sample should be larger than the representative elementary volume (REV) with respect to the studied parameter (Bear, 1988). This important point is addressed in section 3.

To be representative of a macroscopic layer, numerical simulations also need to use a physically relevant parameterization of the data file boundaries. For instance, such a parameterization could be in the form of assuming the borders of the data file are surrounded by ice everywhere; this would provide a very simple Dirichlet condition to the problem. However, such an assumption would result in unrealistic curvature values close to the file borders. More realistic edge conditions can be obtained by combining the Dirichlet and Neumann conditions, for example. Periodization of the data is a more standard approach for similar problems (Yao and others, 1993), especially where a spectral resolution is considered. However, the continuity of second derivatives of the vapour field is broken at the edges even when using mirror symmetries. This generates many unphysical sources and sinks at the edges of the data file, so that some judicious smoothing of the edges still occurs.

As a first step in this direction, we chose to simply discard the outer voxels in the summation (i.e. we averaged over inner voxels only), leading to a truncated expansion of Δp close to the edge. When removing voxels from one side of the neighbourhood (1 side + 4 edges + 4 corners) in a given direction of the grid (say, the xi axis), the corresponding derivatives k/∂xk i of all orders k ≥ 1 are introduced in the Taylor expansion of p side, although they are not accounted for in the digitization of Δp0. Applying the digitization in this way (with ‘missing voxels’) implicitly assumes that these derivatives are set to zero.

The first-order condition is ∂p/∂xi = 0; this is a flux condition that replaces the full 3-D problem with a combination of the Neumann and inner Dirichlet conditions. The second- order condition, 2 p/∂x 2 i = 0, allows mirror symmetries at the edges (Durand, 1966b). This provides a method of turning this composite problem into a sequence of nested one-, two- and three-dimensional (1-D, 2-D and 3-D) inner Dirichlet problems, ensuring the existence and uniqueness of the solution. In using a third-order digitization such as Equation (8), the interface should be smooth (differentiable) when crossing the edge of the data file. This can be thought of as an implicit extension of the data at constant mean curvature. An additional justification for this choice was the fact that κ governs the thermodynamics of both reaction and diffusion problems.

2.5. Algorithm

Following assumption (3) of section 2.1, grain interfaces are considered motionless throughout the computation. The vapour field is computed in the diffusion limit using the following procedure.

  • 1. The field of unit normal vectors n at the surface S of snow grains is adaptively computed (Flin and others, 2005).

  • 2. The corresponding field of mean curvature is computed using κ = (divn)/2 (Flin and others, 2004).

  • 3. The boundary conditions are fixed on the surface S using Equation (2).

  • 4. Similarly, the vapour field is initialized in the pore space D by Kelvin’s equation at p amb = p 0 sat 〈κ〉, 〈κ〉 being the value of the mean curvature averaged over the whole sample.

  • 5. The vapour field is then computed using the iterative Jacobi method (Press and others, 1992). The pressure of the voxel v(r) at iteration (i + 1) is dependent upon the pressures at iteration (i) of the first neighbours of v(r). The method in which is estimated depends on the position of the first neighbours of v(r).

case 1 If all the first neighbours are in the pore space D or on the surface S, digitization Equation (9) is used. We then have:


case 2 If at least one neighbour belongs to the ice, digitization Equation (8) is used. Hence:


In both cases, when at least one neighbouring voxel of v(r) does not belong to the image file, it is not taken into account in the estimation. In other words, the digitization is restricted to the available data by computing the weighted mean of on the inner voxels only (see section 2.4).

6. Procedure (5) is iterated until the field of p vap(r) converges.

The gradient p (r) is then evaluated on the closest non-adjacent layer to S using Prewitt’s digitization of the gradient operator (Prewitt, 1970). The resulting growth rate R is then computed using Equations (6) and (7). If time-lapse evolution of snow microstructure was required, the next step of the simulation would consist of modifying the grain surface S according to the field of surface fluxes. The computation of p(r) and R(r ∈ S+εn) from the newly obtained surface would result, if repeated, in a morphological iterative process that would provide the microstructure evolution of snow with time.

3. Preliminary Results and Outlook

Vapour field and growth rate computations were first checked qualitatively on two spheres of different radii. The test showed matter accretion (positive flux) on the large sphere, close to the small one, unlike the result with a reaction model (Fig. 1). This is the sort of difference we can look for between further simulations of metamorphism using both the reaction and diffusion processes. For both models, p amb was initialized as the saturation pressure, corresponding to

Fig. 1. (a) Qualitative map of the mean curvature κ of a system of two neighbouring spheres of different radii; (b) growth rate R in the reaction case (spheres do not ‘see’ one another); and (c) R in the diffusion-limited case (vapour transfer is favoured between close grains).

The computations were then applied to tomographic 3-D data (Brzoska and others, 1999) acquired from the European Synchrotron Radiation Facility (ESRF, Grenoble, France). The snow images presented in this paper refer to a snow sample obtained at 2 ˚ C after a metamorphosing time of 12 days (Flin and others, 2004). Figure 2 depicts a 2-D slice of the computed vapour field.

Fig. 2. A 2-D slice of the computed vapour-pressure map in the diffusion-limited case. The snow image presented here was obtained after metamorphosing during 12 days under isothermal conditions. Image edges are 300 voxels (3 mm) wide.

In order to check the representativeness of the computed vapour field, the mean values of vapour pressure were computed in cubic neighbourhoods growing from the centre of the 3-D image. The vapour pressure averages vs the size of the computing volumes are plotted in Figure 3. The same process has been applied to two other samples: a recent snow sample obtained after 2 days of isothermal metamorphism (Flin and others, 2004) and a slightly faceted sample taken after 3 weeks of TG metamorphism (TG = 3 Km 1, average temperature –3 ˚ C) (Flin and others, 2007). From Figure 3, it appears that the REV for vapour pressure significantly depends on the snow microstructure. The REV for the 12 day old sample (A) can be approximated as the volume of a cube of 2.5–3.0mm in dimension. However, the REV for the recent snow sample (B) (TG sample (C)) seems much larger (smaller).

Fig. 3. Vapour-pressure REV estimations for different samples.

Figure 4 shows the growth rate results of both reaction and diffusion-limited models for the same 12 day old sample. (The result for the reaction model is scaled by the effective reaction coefficient α that has not yet been precisely assessed from the experimental data.) One can see a sharper transition between positive and negative growth rates in the diffusion process, which is consistent with results on spheres. Local transport is favoured, then enhances the contrast of the crossover between regions of positive and negative R. It will be interesting to check regions where two (convex) grains of different average curvature are close to each other, where the above situation of two spheres is expected.

Fig. 4. (a)Map of the growth rate Rin the reaction case scaled by the effective reaction coefficient α; and (b) map of R in the diffusion-limited case. Image edges are 300 voxels (3 mm) wide. Growth rates are expressed in 10 10 ms 1.

The most interesting comparisons remain to be made between numerically metamorphosed results of both models from the same input snow data. One can reasonably expect visible changes in curvature histograms between these two simulations, and therefore a comparison with the experimental histogram may answer the difficult question of whether the diffusion or the reaction mechanism is dominant for the given experimental conditions. The only requirement is that models should be run on data files that are representative of samples larger than the REV of vapour transfer. Further developments will have to deal with the intrinsic slowness of iterative schemes. It will first be necessary to replace the Jacobi method with a more efficient scheme. For example, over-relaxation methods (Patankar, 1980) could be used to implement the model on nested cubic grids, adopting the low-resolution result as a guess for the final high-resolution run.


We thank J.B. Johnson and an anonymous reviewer for constructive comments. F. Flin is grateful to Y. Furukawa for his welcome at the Institute of Low Temperature Science as a guest researcher.


Bear, J. 1988. Dynamics of fluids in porous media. Mineola, NY, Dover.
Brown, D.E. and 6 others. 1996. H2O condensation coefficient and refractive index for vapour-deposited ice from molecular beam and optical interference measurements. J. Phys. Chem., 100(12), 4988–4995.
Brun, E. and Touvier, F.. 1987. Etude expérimentale de la convection thermique dans la neige. J. Phys. I (Paris), 48(C1), 257–262.
Brzoska, J.-B. and 7 others. 1999. 3-D visualization of snow samples by microtomography at low temperature. ESRF Newsl. 32, 22–23.
Chen, Q.-S., Zhang, H., Ma, R.-H., Prasad, V., Balkas, C.M. and Yushin, N.K.. 2001. Modeling of transport processes and kinetics of silicon carbide bulk growth. J. Cryst. Growth, 225(2–4), 299–306.
Colbeck, S.C. 1983a. Ice crystal morphology and growth rates at low supersaturations and high temperatures. J. Appl. Phys., 54(5), 2677–2682.
Colbeck, S.C. 1983b. Theory of metamorphism of dry snow. J. Geophys. Res., 88(C9), 5475–5482.
Colbeck, S.C. 1997. A review of sintering in seasonal snow. CRREL Rep. 97–10.
De la Vallée Poussin, C. 1932. Extension de la méthode du balayage de Poincaré et problème de Dirichlet. Ann. Inst. H. Poincaré, 2(3), 169–232.
Durand, E. 1966a. Electrostatique, Tome II: problèmes généraux, conducteurs. Paris, Masson.
Durand, E. 1966b. Electrostatique, Tome III: méthodes de calculs, diélectriques. Paris, Masson.
Eaglesham, D.J., White, A.E., Feldman, L.C., Moriya, N. and Jacobson, D.C.. 1993. Equilibrium shape of Si. Phys. Rev. Lett., 70(11), 1643–1646.
Flin, F., Brzoska, J.-B., Lesaffre, B., Coléou, C. and Pieritz, R.A.. 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.
Flin, F., Brzoska, J.-B., Lesaffre, B., Coléou, C. and Pieritz, R.A.. 2004. Three-dimensional geometric measurements of snow microstructural evolution under isothermal conditions. Ann. Glaciol., 38, 39–44.
Flin, F. and 9 others. 2005. Adaptive estimation of normals and surface area for discrete 3-D objects: application to snow binary data from X-ray tomography. IEEE Trans. Image Process., 14(5), 585–596.
Flin, F., Brzoska, J.-B., Pieritz, R.A., Lesaffre, B., Coléou, C. and Furukawa, Y.. 2007. The temperature gradient metamorphism of snow: model and first validations using X-ray microtomographic images. In Kuhs, W.F., ed. Physics and chemistry of ice. Cambridge, RSC Publishing, 181–189.
Frey, P.J. and George, P.L.. 2000. Mesh generation: application to finite elements. Oxford, etc., Hermes Science.
Kaempfer, T.U., Schneebeli, M. and Sokratov, S.A.. 2005. A micro-structural approach to model heat transfer in snow. Geophys. Res. Lett., 32(21), L21503. (10.1029/2005GL023873.)
Legagneux, L. and Dominé, F. 2005. A mean field model of the decrease of the specific surface area of dry snow during isothermal metamorphism. J. Geophys. Res., 110(F4), F04011. (10.1029/2004JF000181.)
Libbrecht, K.G. 2003. Growth rates of the principal facets of ice between 10 ˚ C and 40 ˚ C. J. Cryst. Growth, 247(3–4), 530–540.
Mutaftschiev, B. 2001. The atomistic nature of crystal growth. Berlin, etc., Springer-Verlag.
Patankar, S.V. 1980. Numerical heat transfer and fluid flow. New York, Hemisphere Publishing.
Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P.. 1992. Numerical recipes in C: the art of scientific computing. Second edition. Cambridge, Cambridge University Press.
Prewitt, J.M.S. 1970. Object enhancement and extraction. In Lipkin, B.S. and Rosenfeld, A., eds. Picture processing and psychopictorics. New York, Academic Press, 75–149.
Sekerka, R.F. 2004. Theory of crystal growth morphology. In Müller, G., Métois, J.-J. and Rudolph, P., eds. Crystal growth: from fundamentals to technology. Amsterdam, Elsevier, 55–93.
Yao, J., Frykman, P., Kalaydjian, F., Thovert, J.F. and Adler, P.M.. 1993. High-order moments of the phase function for real and reconstructed model porous media: a comparison. J. Colloid Interface Sci., 156(2), 478–490.
Yosida, Z. and 6 others. 1955. Physical studies on deposited snow. I. Thermal properties. Contrib. Inst. Low Temp. Sci., Ser. A 7, 19–74.