Articles | Volume 15, issue 9
The Cryosphere, 15, 4381–4398, 2021
The Cryosphere, 15, 4381–4398, 2021

Research article 10 Sep 2021

Research article | 10 Sep 2021

Orientation selective grain sublimation–deposition in snow under temperature gradient metamorphism observed with diffraction contrast tomography

Orientation selective grain sublimation–deposition in snow under temperature gradient metamorphism observed with diffraction contrast tomography
Rémi Granger1,2, Frédéric Flin2, Wolfgang Ludwig3, Ismail Hammad3, and Christian Geindreau1 Rémi Granger et al.
  • 1Univ. Grenoble Alpes, CNRS, Grenoble INP (Institute of Engineering Univ. Grenoble Alpes), 3SR, 38000 Grenoble, France
  • 2Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Etudes de la Neige, 38000 Grenoble, France
  • 3ESRF, BP 220 38043 Grenoble CEDEX, France

Correspondence: Rémi Granger ( )


In this study on temperature gradient metamorphism in snow, we investigate the hypothesis that there exists a favourable crystalline orientation relative to the temperature gradient, giving rise to a faster formation of crystallographic facets. We applied in situ time-lapse diffraction contrast tomography on a snow sample with a density of 476 kg m−3 subject to a temperature gradient of 52 Cm-1 at mean temperatures in the range between 4.1 and 2.1 C for 3 d. The orientations of about 900 grains along with their microstructural evolution are followed over time. Faceted crystals appear during the evolution, and from the analysis of the material fluxes, we observe higher sublimation–deposition rates for grains with their c axis in the horizontal plane at the beginning of the metamorphism. This remains the case up to the end of the experiment for what concerns sublimation while the differences vanish for deposition. The latter observation is explained in terms of geometrical interactions between grains.

1 Introduction

The snow cover consists of a porous medium made of ice crystals, air and sometimes liquid water. The microstructure of this material is in constant evolution, and the original fresh snow is transformed under the action of various mechanisms called snow metamorphism. In particular, a temperature gradient (TG) in a snow layer implies a gradient of water vapour density, which creates a vapour flux in the pore space. Under these circumstances, ice crystal interfaces may be significantly out of equilibrium with the vapour phase if the macroscopic temperature gradient is higher than 10–20 Cm-1. Of course, these values should just be considered as orders of magnitude, since the local temperature gradient is strongly depending on the geometrical configuration at the grain scale, as e.g. pointed out by Hammonds et al. (2015) and Hammonds and Baker (2016). Growth and decay under TG conditions depend strongly on the kinetics of the deposition and sublimation processes on different crystallographic facets, and as a consequence, kinetic forms (faceted crystals and depth hoar) of snow crystals appear (Colbeck1983). An important aspect of temperature gradient metamorphism is that the grain size increases while the density remains constant (Marbouty1980; Schneebeli and Sokratov2004; Srivastava et al.2010; Pinzer et al.2012; Calonne2014). Indeed, due to the local temperature field, some parts of the microstructure grow at the expense of others (Akitaya1974; Pinzer et al.2012), leading to an almost null net flux of ice far from the system boundaries. In addition, the kinetic coefficient is not equal on the prismatic and basal faces of ice (see e.g. Yokoyama and Kuroda1990; Libbrecht2005; Furukawa2015). In consequence, growth and decay of an ice crystal in a vapour gradient depend on the orientation of the crystal relative to the thermal (and also water vapour saturation pressure) gradient, with higher deposition and sublimation rates when the more active faces efficiently catch and release the flux.

In that context, Adams and Miller (2003) and later Miller and Adams (2009) conjectured that, in temperature gradient metamorphism, favourably oriented crystals might grow at the expense of others as they grow more quickly towards the common vapour source. Because of that selective growth of crystals, the distribution of crystallographic orientation, called the fabric, is expected to evolve towards an anisotropic fabric. However, few studies report data on the subject. Takahashi and Fujino (1976) and de Quervain (1983) measured crystal orientation, but the data were too scarce to draw conclusions.

More recently, Riche et al. (2013) measured the fabric evolution over natural and sieved samples under a 50 Cm-1 temperature gradient for 12.5 weeks, at a mean temperature of 20 C. They used the Automatic Ice Texture Analyser (AITA) (Wilson et al.2007), which enables measurement of crystallographic orientations of thin sections. They observed an evolution of the fabric from strong cluster type (with c axis of ice crystals mainly in one direction) to weak girdle type (with c axis mainly in a plane) for the natural samples, while no evolution for the denser sieved sample was observed. Calonne et al. (2017) analysed snow–firn profiles from Antarctica in terms of (i) microstructure, using microcomputed tomography, and (ii) fabric, using AITA. They showed that the snow fabric is correlated with its microstructure, with the most isotropic fabric being found in the densest snow and smallest specific surface area (SSA). A similar study applied to EastGRIP data (Greenland) revealed some strong cluster-type textures at shallow depth (Montagnat et al.2020). All of these recent studies suggest that metamorphism affects the fabric of snow and firn.

While AITA is an efficient instrument to obtain statistical data on fabric and on different samples, it requires thin sectioning of the sample. Consequently, it is not possible to follow the time evolution of individual grains and observe the role of crystal orientation at the grain scale. To bring such complementary understanding into the picture, we present here an experiment of temperature gradient metamorphism, monitored over a duration of 3 d by conventional tomography and diffraction contrast tomography (DCT). DCT is a near-field X-ray diffraction imaging modality exploiting the Bragg diffraction signals created by the individual crystals in the illuminated sample volume. By measuring the position and shape of the associated diffraction spots, it is possible to reconstruct the crystallographic orientation and shape of each of the grains (Ludwig et al.2009; Reischig et al.2013). This technique allows for non-destructive in situ observations of extended sample volumes and has already been used for snow to follow deformation under compression (Rolland du Roscoat et al.2011).

The paper is structured as follows: in the first section, we describe the sample preparation, the experimental setup and the principal image processing steps. In the second section, we present a time series of the evolving sample microstructure as well as the characterization and analysis of the relation between grain orientation and mass fluxes. Finally, we discuss the observed trends, the representativeness of the experiment and the main technical difficulties met with the imaging technique.

Figure 1 Scheme of the cryogenic cell CellDyM2, the DCT-compliant version of CellDyM (Calonne et al.2015).


2 Methods

2.1 Sample and sample environment

To impose controlled temperature conditions on the sample, we used CellDyM2 (see Fig. 1), a modified version of the cryogenic cell developed by Calonne et al. (2015): as for CellDyM, the snow specimen is inserted in a sample holder placed between two copper columns, at the extremities of which temperature is controlled using Peltier modules. The sample holder and the copper columns are placed in a vacuum chamber for thermal insulation from room temperature. However, CellDyM2 was specifically adapted to DCT experiments. Indeed, the original aluminium sample holder, with typically large crystalline grains, can create parasitic diffraction spots and thus has been replaced by a system of three components:

  • a 1 mm thick polymethyl methacrylate (PMMA) holder of inner diameter 10 mm,

  • a 3 mm thick graphite cylinder of inner diameter 4 mm, and

  • a thin polypropylene (PP) tube inside the graphite cylinder.

The PMMA container, combined with two rubber O-rings, provides a vacuum-tight assembly. The graphite, of thermal conductivity 95 Wm-1K-1 and grain size below ≈5µm, ensures a good thermal conductivity between the two copper columns and, thus, a good control of the temperature gradient without problematic parasitic diffraction spots. Finally, the PP cylinder avoids the diffusion of water vapour inside the porous graphite and preserves the sample from potential contamination by carbon dust. The vacuum in the chamber is ensured by a system of two pumps: a primary scroll pump (Anest Iwata ISP-90) produces a pressure of about 500 Pa, which permits a turbomolecular pump (Pfeiffer Vacuum HiPace 30), directly mounted on the turntable, to reach the final pressure of 0.009 Pa.

Figure 2 Evolution of the mean temperature and downward temperature gradient imposed on the sample.


The cylindrical sample volume with 4 mm diameter and 4 mm height contained snow of rounded grain (RG type, according to the international classification, Fierz et al.2009) with a density of ≈450kg m−3 and was submitted to a downward temperature gradient. The choice of this snow type was guided by our objective to easily measure local phase change fluxes against the orientation of crystals relative to the temperature gradient direction. For example, working on fresh snow would have led to possible image resolution issues as well as to a rearrangement of the grains during the DCT experiment, thus complexifying the following analysis. To ensure constant vapour flux boundary conditions to the snow sample, the top and bottom surfaces of the whole 10 mm high cylinder consisted of 3 mm thick ice lenses, serving as a matter sink or source, respectively. In order to minimize parasitic diffraction spots, these lenses were machined from monocrystals, using a specific core drill. The RG sample was prepared in a cold room at the Centre d'Etudes de la Neige (CEN): the original fresh snow, picked up at the upper station of La Grave cable car, was stored at 20 C for 7 months before it was metamorphosed at 5 C in equilibrium conditions for 3 d. A specimen was then cored in a homogeneous zone and assembled on the bottom copper column with the bottom ice lens. The top ice lens was then stuck to the top copper column, which was hence placed at the upper part of the sample. The assembly was then stored in a specific copper handling tube, preserving the sample from possible temperature gradients during storage. The whole assembly was then transported to the European Synchrotron Radiation Facility (ESRF) within a portable freezer and inserted in the cell at 20 C.

As soon as the cell was closed, the pumping of the air out from the vacuum chamber started, leading to a pressure of 0.046 Pa in less than 27 min. It took ≈7h to reach the stable pressure of 0.009 Pa. Then, the regulation of both of the Peltier modules was progressively set to 2.10 C. This value was reached in about 1 h, and the sample stayed in equilibrium conditions for 21 additional hours. Then, the proper temperature gradient experiment began: this corresponds to the instant t=0 in this paper.

The thermal conditions applied to the snow sample during the experiment are shown in Fig. 2. The experiment consists of three phases that we will now refer to as phase 1, phase 2 and phase 3. Mean temperature and temperature gradients were respectively set to 2.09 C and 17 Cm-1 during phase 1, 2.09 C and 52 Cm-1 during phase 2, and 4.07 C and 52 Cm-1 during phase 3. These conditions were determined by performing numerical simulations of heat transfer in the cell using COMSOL Multiphysics, as presented by Calonne et al. (2015). In all these simulations, the temperatures imposed by the Peltier modules were used as inputs in the numerical model.

2.2 Diffraction contrast tomography: experimental setup

The experiment was performed at beamline ID19 of the ESRF in a period of four-bunch operation of the storage ring. A beam energy of 40 keV was selected using a Silicon double crystal monochromator in Laue–Laue setting. Detection was performed using a detector system consisting of a 200 µm LuAG screen coupled via visible light optics to a sCMOS camera (2048×2048 pixels), resulting in an effective pixel size of 6.5 µm and a field of view of 13 mm. The distance between the rotation axis and the detector was 60 mm.

DCT scans consisted of 7200 projection images, acquired during continuous rotation of the sample over 360 with an angular step size of 0.05 and exposure time of 50 ms, resulting in a total scan time of about 1 h. In order to reduce background, the area of the scintillation screen illuminated by the direct beam was covered by a 5×5mm2 beamstop, made from 0.5 mm lead sheet. Absorption scans, acquired with the same detector, shifted laterally by 5 mm, consisted of 900 images with angular step size of 0.4 and 0.05 exposure time (≈10min scan time). DCT and absorption scans were performed alternatively so that there was systematically an absorption scan performed in the 20 min before a DCT scan.

Reconstruction was performed following the workflow detailed in Reischig et al. (2013). A brief outline of the processing route is summarized below.

  1. Preprocessing. Distortion and background correction are applied to the image stack.

  2. Segmentation of diffraction spots. Double soft thresholding is applied.

  3. Matching Friedel pairs. When a crystalline plane (hkl) meets the Bragg's condition and leads to a diffraction spot for a position ω of the rotation stage, then the Bragg condition is met again at ω+180 for the (h¯k¯l¯) plane and the two associated diffraction spots are called a Friedel pair, which are identified automatically using criteria based on the position, intensity and shape of the segmented diffraction spots. Over 360, a crystallographic plane in one grain can meet the Bragg conditions up to four times and so gives rise to two Friedel pairs. A Friedel pair allows inferring the trajectory of the diffracted beam and hence deriving the direction of the (hkl) plane normal.

  4. Indexing grains. The position and orientation of grains are determined by identifying consistent Friedel pairs of diffraction spots corresponding to the same grain.

  5. Final spot collection. Additional diffraction spots by forward simulation and final selection of diffraction spots are collected.

  6. Reconstruction. The 3D grain shape using the selected diffraction spots for each grain using the SIRT algorithm is reconstructed (Palenstijn et al.2013, 2015).

  7. Assembly of the 3D grain volume. The individually reconstructed grain volumes are assembled into the common sample volume.

  8. Masked dilation. A morphological dilation of the reconstructed grain map is applied. Such a dilation is constrained to the domain occupied by ice, as determined by segmentation of the sample volume obtained from the absorption contrast scans.

Figure 3 Effect of step 8 to improve the reconstruction quality: before (a) and after (b) correction.


The last step of masked dilation improves the accuracy of shape reconstruction in porous (and/or multiphase) media by enforcing fidelity of the grain shape at external grain interfaces (or phase boundaries). Ice crystals in snow are nearly perfect and in our case of dimensions comparable to or bigger than the Pendellösung lengths: as a consequence, dynamic (multiple) diffraction leads to noticeable deviations of the local diffracted intensity from the idealized, linear projection model (kinematic diffraction conditions) upon which our implementation of the reconstruction algorithm is based. These non-linearities and remaining diffraction spot overlaps, not filtered out in previous consistency checks, give rise to inconsistent grain shape reconstructions, i.e. streaks and modulations of the reconstructed local scattering power inside the grain volumes. Since the final grain map is assembled from a segmentation of these individual greyscale grain volumes, these artefacts can result in (i) inaccurate contours created by overlap between neighbouring grains, (ii) inaccurate contours created by overlap with the pore space of the material, and (iii) internal holes or holes connected to the contour of the segmented grains.

Since the localization of ice phase can be easily obtained from the absorption scans, the masked dilation can be used to correct part of the inaccuracies introduced by dynamical diffraction (for instance internal holes). Figure 3 shows the effect of this correction on a horizontal slice of the sample. While this procedure enables recovery of the correct information for a large fraction of the volume elements (voxels), one has to be aware that it may also lead to the assignment of wrong orientations for voxels in the vicinity of the grain boundaries. This point will be discussed in Sect. 4 in the light of the obtained results.

2.3 Data processing

Reconstruction of the absorption scans was performed using the filtered back projection (FBP) algorithm, and ring artefacts were corrected using an in-house algorithm.

Because some movements of the stages and detector may occur during data acquisition, it is necessary to register the time series of the reconstructed volumes. To simplify that step, reference grooves have been previously made on the exterior side of the PP cylinder following both the axis and the circumference. First, all absorption volumes were registered on a reference volume (taken at t=27h 32 min, at the beginning of phase 2). Registration was performed using the SimpleITK library and consists of

  • smoothing the image with a Gaussian filter with kernel size of two voxels,

  • finding the optimal translation in the horizontal plane,

  • finding the optimal vertical translation, and

  • applying the total transform to the not smoothed reconstructed volumes with linear interpolation.

The optima were defined by maximizing the correlation computed in a hand-selected region of interest ROI focusing on the grooves. Registered volumes were then segmented using an energy-based segmentation algorithm (Hagenmuller et al.2013).

Some grains may be missing during the DCT reconstruction because the diffracted beam may not hit the detector or because the quality and quantity of grouped diffraction spots may not be sufficient for reconstruction. Correcting such data requires manual monitoring and processing throughout the sample volume. For that reason, only a selected subset of the DCT scans was processed in that study. Times at which a processed volume is available are represented by grey dotted lines in Fig. 2. The correction consists in comparing the reconstructed absorption and DCT scans to detect missing grains. This is done by comparing the former and latter volumes in the time series to get the orientation of the grains which remain fixed during the experiment, as no rearrangement occurs.

After that, because DCT scans were always performed shortly after the absorption scans, they were directly registered on the closest absorption volume by finding the optimal translation in space that maximizes the correlation between the two volumes, without ROI selection.

Figure 4 Time evolution of the microstructure in a cube of 2.34 mm edge length cropped out of the centre of the specimen. The interfaces are coloured based on the local mean curvature.


Figure 5 Grain map evolution. Each colour represents a single crystalline orientation. The meaning of the colours is given by the colour wheel, representing orientations in a pole figure centred on a vector X+Y+Z. The five axes under the colour wheel show the colour corresponding to the main directions associated with X, Y and Z.


Figure 6 Parametrization of the different unit vectors and angles: z is the vertical component of the image frame, n is the local normal outgoing from the ice and c is a vector with positive vertical component representing the c-axis direction.


3 Results

3.1 Image series

Figure 4 shows images of the microstructure at t=25.1h and t=78.6h of a cubic subsample of edge length 2.34 mm in the centre of the sample. The colours represent the mean curvature computed with an algorithm developed by Flin et al. (2004, 2005) – see also Calonne et al. (2014), Haffar et al. (2021) and ranges from 46 to 46 mm−1. First, we can see that the initial microstructure of the snow consists of small well-rounded grains. The global evolution is small as full recrystallization of the grains is not observable. Nevertheless, all the typical features of temperature gradient metamorphism can be observed, e.g.

  • a strong faceting of the downward-oriented surfaces, as revealed by the grain geometry and the associated curvature field, and

  • a global migration of the ice microstructure towards the warmer side of the sample.

In Fig. 5, we present the result of the DCT acquisition obtained after the processing route presented in Sect. 2.3. The orientation of each grain is characterized by a unit vector c which, by convention, has a positive vertical component. On the figure, the direction of c axis is represented by a colour given by the colour wheel. From a grain-to-grain close-image comparison, we can see that no grain rearrangement can be observed (no c-axis rotation) but just a modification of the grain shapes and boundaries, which is consistent with the current knowledge of TG snow metamorphism.

Figure 7 Time evolution of the facet orientation. The colours describe α, the angle between the c axis obtained from the DCT information and the normal to the interface: green represents basal faces, and red corresponds to cases where the surface is parallel to the c axis. White arrowheads point out plate-like crystals formed at the end of the experiment, where consistency between orientation and morphology can be checked.


For further analysis, we introduce the parametrization represented in Fig. 6. At a given point of the surface, three vectors have importance in the physics of the metamorphism: (i) the vertical direction, given by z pointing upward, which is the direction of the global macroscopic vapour flux; (ii) the normal at the interface n, taken here as pointing into the pore space; and (iii) the c axis given by c. We can then define three angles: α=(c,n), β=(z,c) and γ=(z,n). Because we are only interested in the direction of c, the range of values is [0, 90] for α and β, and it is [0, 180] for γ.

Figure 7 represents the microstructure with the interface coloured by the angle α so that the basal faces (α=0) are green and the prismatic faces (α=90) are red. By looking particularly at the plates formed during the metamorphism, we can observe that the orientation of the c axis is consistent with their morphology.

Figure 8Density measurements. (a) Representative elementary volume estimation for density, computed from centred cubic subvolumes of variable sizes. (b) Vertical variation of density computed on centred discs of one-voxel height and 3.3 mm in diameter.


Figure 9 Density evolution: grey vertical lines represent the times of change of temperature conditions. Red dots correspond to times for which DCT information is available. The horizontal black dashed line represents mean density, which stays around 476kg m−3 (±5 %).


3.2 Density

Representative elementary volume (REV) analysis for density was performed on absorption volumes at different times of the experiment, using cubic centred subvolumes of increasing sizes. The results are depicted in Fig. 8a and show that the density reaches a constant value for cubes of edge about 2.34 mm. In order to complement the characterization of the homogeneity, we computed the density of horizontal discs of the sample, centred on the vertical axis of the sample with a diameter of 3.3 mm (so that it corresponds to the cubic REV of 2.34 mm edge) and one voxel (6.5 µm) high. This enables one to obtain a vertical density profile. The results are shown in Fig. 8b. The high values on top and bottom correspond to ice lenses (ρice=917kg m−3). Height and time fluctuations are visible, but the sample is overall homogeneous in the vertical direction. We also computed the density for each absorption volume on a cubic REV of size 2.34 mm. The result is shown in Fig. 9. We observe a constant density throughout the experiment at the initial value of (476±5 %)kg m−3.

Figure 10 Specific surface area: REV estimation for different times. The vertical grey line represents the size of the REV volume taken for the remainder of the study.


Figure 11 Evolution of the specific surface area during the experiment. The grey vertical dashed lines represent the times of change of temperature conditions. Red dots correspond to times for which DCT information is available.


3.3 Specific surface area

First, REV for SSA was computed on centred cubic subvolumes of increasing sizes, at various times during the experiment. SSA computations were performed using an improved version (see Dumont et al.2021) of the VP projection algorithm (Flin et al.2005, 2011). Figure 10 summarizes the results of this analysis. For all times, values are stable for subvolumes of size larger than 2 mm. So the REV of density (2.34 mm) can be chosen.

Figure 11 shows the evolution of SSA during the experiment. Starting at (18±5 %)m2 kg−1, the SSA decreases slightly during the first 30 h and then increases more and more rapidly to reach about (20±5 %)m2 kg−1.

Figure 12 Evolution of the eigenvectors of the orientation tensors represented in the pole figure with vertical axes as reference (a) and eigenvalues (b). Colours of dots for eigenvector on the left correspond to eigenvalue colours on the right.


3.4 Orientation tensor

Given the crystallographic orientation of each ice voxel, we can compute the second-order orientation tensor (Fisher et al.1993; Riche et al.2013):

(1) a ( 2 ) = 1 N i i = 1 N i c i c i ,

with ci being the c-axis orientation of the ith voxel in the volume, the tensor product and Ni the number of ice voxels. The orientation distribution is inscribed in an ellipsoid whose principal axes are given by the eigenvectors (v1, v2, v3) of a(2) and the length of the principal axis by its eigenvalues. Denoting the eigenvalues as a1, a2, a3 with 1>a1a2a3>0, we can identify three limit cases, defining three references of fabric.

  • a1>a2a3: cluster-type fabric, majority of c axis in one direction (defined by v1).

  • a1a2>a3: girdle-type fabric, majority of c axis in a plane (defined by v1 and v2).

  • a1=a2=a3=1/3: isotropic fabric, same distribution in the three directions.

Figure 13 Cluster index versus girdle index for the different times indicated. Grey dashed line represents the y=x line.


In order to evaluate the proximity of the measured distribution to those types, cluster index, girdle index and strength can be defined by ln(a1/a2), ln(a2/a3) and ln(a1/a3) respectively. All indexes equal 0 for a perfectly isotropic fabric. Figure 12 represents eigenvectors on the left panel, in pole figure representation, while the associated eigenvalues for the DCT images against time are shown on the right-hand side. Eigenvalues lie in the range 0.279 to 0.405 without significant visible evolution over time. Eigenvectors are stable with time, with one being parallel to temperature gradient and the two others being contained in the horizontal plane.

Figure 13 is a scatter plot of the cluster index against the girdle index for the treated DCT volumes. The cluster index, girdle index and strength index are respectively 0.188, 0.156 and 0.343 at the initial stage and 0.266, 0.083 and 0.349 at the final stage, with mean values of 0.217, 0.117 and 0.335. Thus, the fabric does not really change during the experiment, and all orientations are similarly represented through the sample.

3.5 Pore size distribution

Pore size distributions were computed on absorption images using a sphere fitting algorithm. The results are shown for several times in Fig. 14, with vertical lines denoting the mean of the corresponding distribution. Overall, the mean stays at about 190 µm and the distribution does not show any drastic evolution in shape. In particular, no trend is observable for the biggest pores, while a small trend is observable for the smaller ones. Indeed, the quantity of pores of size in the range 30 to 80 µm increases from about 3 % to 4 % of the total pore volume, which represents about one-third of its original value.

Figure 14 Pore size distribution, with two voxels = 13 mm bin size computed using a fitting sphere algorithm. Vertical dotted lines represent the mean of the associated distribution.


3.6 Relation between flux and grain orientation

The main objective of this study is to look at possible relations between local fluxes and crystalline orientation. In particular, we would like to test at the local scale the validity of orientation selective growth of grains under temperature gradient metamorphism.

Figure 15 Interpretation of the diagram of mean flux against α=(c,n) and γ=(z,n). The small red arrows represent the facets taken into account.


Figure 16 Interpretation of the diagram of mean flux against β=(z,c) and γ=(z,n). The small red arrows represent the facets taken into account.


From the 30 registered and segmented absorption images (Ak)k1,30, we computed local interface speeds vk, where k is the index of the absorption image. At a given point of the surface of the volume k,30 the normal nk was computed using the ADGF algorithm (Flin et al.2005). Denoting by dkl the displacement of the interface along nk from Ak to Al, we computed both dkk+1 and dkk-1 using a specific algorithm (Flin et al.2014, 2018) whose principle is relatively close to the one presented by Krol and Löwe (2016). In that case, the error on the measured distance is mainly due to errors on the interface position induced by the discretization and is estimated to one voxel = 6.5 µm. Then local normal velocity vk at the instant tk of acquisition of Ak was estimated using

(2) v k = d k k + 1 t k + 1 - t k  if  k = 1 t k - t k - 1 t k + 1 - t k - 1 d k k + 1 t k + 1 - t k - t k + 1 - t k t k + 1 - t k - 1 d k k - 1 t k - t k - 1  if  k 2 , 29 - d k k - 1 t k - t k - 1  if  k = 30 .

From the normal velocity of the interface, the mass flux ϕk of vapour changing into ice can be computed using

(3) ϕ k = ρ i v k .

Voxels were first classified depending on the values of α and γ in 2D bins (of size 0.25 and 1, respectively). Then, a mean flux value was computed for each bin. This enables us to represent a colour map of the mean fluxes against α=(c,n) and γ=(z,n) at the different times of the experiment. An example of this representation is given in Fig. 15. The pixels on the right of the diagram represent prismatic faces (α=90), while the pixels on the left represent basal faces (α=0). The top of the diagram corresponds to surfaces pointing upward (γ=0), while the bottom of the diagram represents surfaces pointing downward (γ=180); pixels at mid-height of the diagram correspond to surfaces whose normals are horizontal (γ=90). A similar representation can be done by representing the mean flux against β=(z,c) and γ=(z,n) (with respective bin size of 2 and 1; see Fig. 16). In that case, the right of the diagram represents grains with a horizontal c axis (β=90), while the left of the diagram corresponds to grains with a vertical c axis (β=0).

Figure 17 Mean fluxes against α=(c,n) and γ=(z,n).


Figure 18 Mean fluxes against β=(z,c) and γ=(z,n).


Figure 19 Mean fluxes against α=(z,c), for γ[0,20] (faces oriented upward), γ[80,100] (face oriented vertically) and γ[160,180] (face oriented downward). Shaded areas correspond to associated uncertainties, coming from the measurement of the interface displacement between two images and estimated to be one voxel, i.e. 6.5 µm.


Representations for all times of mean flux against (α, γ) are given in Fig. 17 and against (β, γ) in Fig. 18. To help in the interpretation of these diagrams, Fig. 19 presents the values shown on Fig. 17 averaged over three ranges of γ: [0,20], [80,100] and [160,180], against α. The shaded areas represent error arising from the measurement of the distance between two successive interface positions, which is estimated to be one voxel, i.e. 6.5 µm.

At t=25.1h, from all figures, the bin-mean fluxes have modulus less than 2×10-7kgm-2s-1 in all cases. At t=33.3h, mean fluxes for γ<90 are negative (denoting sublimation) and positive (deposition) otherwise. This is due to the fact that the mean vapour fluxes through the pore space are oriented upward: the interfaces whose normals point downward (γ>90) catch vapour which deposits on ice into facets, while the interfaces oriented upward sublimate with rounded shapes (see e.g. Fig. 4 and Calonne et al.2014; Krol and Löwe2018).

In addition, we can observe a second variation depending on the value of α. Higher rates of sublimation–deposition with absolute bin-mean fluxes of 5×10-7kg m2 s−1 are indeed found for α>45. For example, from Fig. 19, for α varying from 5 to 90, one can see a trend for the fluxes varying from 4×10-7 to 5×10-7kg m2 s−1 for γ in [160, 180] and from -4×10-7kg m2 s−1 to -5×10-7kg m2 s−1 for γ in [0, 20]. While there is no dependence on α for γ=90, the closer γ gets to 0 or 180, the stronger is that dependence (see Fig. 17). Until t=78.8h, the same dependence on γ is observed, and the dependence on α has the same shape for γ<90, while it is not visible for γ>90. Overall, at all points of the diagram, the values are lower than at t=33.3h. The evolution between t=33.3h and t=78.8h is regular so that times t=49.9h, t=59.9h and t=71.8h are intermediate steps. The same observations can be made in Fig. 18. The dependence on γ differentiates between faces subject to deposition and those subject to sublimation. For t≥33.3h, voxels belonging to crystals with a c axis in the horizontal plane or close to this plane (β[45,90]) are more active.

Figure 20 Evolution inside a small subvolume of 520×520×520µm3. From left to right, columns represent the c-axis orientation, α, β, γ and the computed fluxes at the interface. The rows correspond to the different times of the experiment as indicated on the left of the figure.


3.7 Example of evolution inside a small subvolume

Figure 20 illustrates the various points observed from the data analysis. It shows the evolution of the microstructure inside a small cubic subsample of edge size 80 voxels = 520 µm, along with different local quantities: c-axis orientation, α, β, γ and interfacial fluxes. For α, the colour map is so that prismatic and basal faces are represented in red and green, respectively. We can observe the formation of three faceted crystals, each with a different orientation. The one pointed out by the red ellipse is almost horizontal, while the one pointed out by the white ellipse is almost vertical. The blue localized facet has intermediate orientation. For all of the three plates one can see that the crystals preferentially grow by their prismatic faces (red area in second column).

The points of the prismatic surface in the area contoured by the red ellipse have coordinates (90, 90) in the (α, γ) space (Fig. 20) and weaker flux than those of the surface pointed out by the blue ellipse, whose coordinates are (90, 135) in that same space. Also, we can see some interactions between the plates pointed out by red and white ellipses. In particular, at t=78.8h, we observe that the vertical plate (white ellipse) grows at the expense of the horizontal one. Thus, while it is not a representative volume of the microstructure, this example illustrates well the observed trends in Figs. 17, 18 and 19.

4 Discussion

4.1 Fluxes and dependence on orientation

At t=25.1h, the fluxes and variations on fluxes are small. The temperature gradient of 17 Cm-1 is the smallest used. We conclude that in phase 1, metamorphism did not act significantly due to small temperature gradient. This can be explained by the relatively high density of the sample of 476 kg m−3. Indeed, it is known that high-density snow requires higher temperature gradients (Akitaya1974; Marbouty1980; Pfeffer and Mrugala2002) and longer times of evolution to get visible modifications of the microstructure. Under a 52 Cm-1 temperature gradient, the maximum values of the bin-mean fluxes are of about 4×10-7kg m2 s−1. Pinzer et al. (2012) computed a macroscopic flux of about 3.5×10-7kg m2 s−1 for snow evolving under 49 Cm-1 temperature gradient at 3.4 C. Flin and Brzoska (2008) computed the highest fluxes of about 1.6×10-7kg m2 s−1 under 16 Cm-1. Thus, the present computed flux values are in very good agreement with this literature.

As described in Sect. 3.6, for t=33.3h, we can see a trend for the most active points of the interface to be with α45. We think that this corresponds to plates growing from grains, with prismatic faces being more active than basal faces. Indeed, from the definition of α, the prismatic faces correspond to α=90. Additionally, at the beginning of the experiment, the grains are typically rounded so that all values of α are equally represented in terms of surface area. Since surfaces with α=45, which have high index orientations, are rougher at the molecular level, they have higher deposition coefficient and are more active. Moreover, in this experiment, grown plates were really thin: about two to three voxels thick; as a consequence, normals computed on the side of these structures were inaccurate, leading to several tens of degrees of uncertainty. These two effects explain why the highest fluxes are found in the range α 45 to 90 and not just at α=90. This observation is consistent with the mean temperature during the experiment, which was always higher than or equal to 4 C. Indeed, under that temperature condition, the kinetic coefficient is expected to be higher on prismatic faces than on basal faces (Kuroda and Lacmann1982; Furukawa2015). Correspondingly, from Fig. 18, the most active faces correspond to crystals for which β=90, i.e. with the c axis in the horizontal plane. Following Adams and Miller (2003), this seems to confirm that kinetically favourably oriented crystals, in our case with basal faces oriented vertically and, thus, showing prismatic faces downward and upward, are indeed more active. At the end of the experiment, the fluxes at the interface overall decrease. We see three reasons for that: firstly, the proportion of high index orientations (with α45) disappear as the facets grow. Secondly, after the important growth of plates between t=25h and t=33h, the growth, and the whole metamorphism, slows down because of geometric interactions between grains (see e.g. Fig. 14, which shows the average pore size is slightly decreasing with time). Finally, at the end of the experiment, the mean temperature was reduced, thus reducing the sublimation–deposition rates.

4.2 End of experiment: independence of orientation for deposition

While a strong phase change is observed for grains with c axis close to the horizontal plane at t=33.3h, this progressively disappears on deposition surfaces. This suggest that at a certain point, growth intensity becomes less dictated by kinetics. We think that this may be explained by grain-to-grain interactions. At the beginning of phase 2, the microstructure has not evolved much, and a strong gradient triggers the growth of plates through the available pores. When the facets grow in the pore space between several grains, growth locally occurs if the temperature of the grain is lower than the average temperature of the grain in interaction with that pore. But, as the facet grows through the pores, it becomes closer and closer to another grain, and growth becomes more dependent on the local configuration of the microstructure and less in relation to the macroscopic temperature gradient. So, faceted grains develop well in pores and less efficiently after they have grown through the pores (see e.g. Fig. 14). Similar conclusions on the necessity of pore space for growth have been drawn by Akitaya (1974). However, on the sublimating surfaces, as the interface moves away from the pore, such a phenomenon does not occur, and fluxes stay well dependent on orientation.

4.3 Uncertainties due to processing

As explained in Sect. 2.3, a dilation step is required at the end of the image processing route, in order to recover the orientation information missing due to dynamic diffraction. One can wonder if that could bias the Figs. 17 and 18. When dilating the orientation information in one partially reconstructed grain, the correction is correct until the dilation front crosses a grain boundary. Beyond this, the correction is wrong. However, we assume the orientation difference between adjacent grains to be random. In consequence, that step affects the signal-to-noise ratio, which increases if the number of corrected voxels is higher than the number of errors. Errors are located in the vicinity of grain boundaries and represent a very small fraction of the analysed surface voxels, and this implies no systematic error is able to explain the observed trends. However, the improvement of the reconstruction assessing the problem of dynamic diffraction is of great interest, since it would increase the signal quality and permit other analyses. In particular, it seems crucial to properly study the grain boundaries. The problem of dynamic diffraction could be tackled by (i) improvement of diffraction spot thresholding and (ii) modelling of the dynamic diffraction to reconstruct the grains from inhomogeneous intensities. The latter is a current subject of development of the technique.

4.4 Representativeness

REV analysis for density and specific surface area shows that a 2.34 mm cubic subsample was representative of the volume sample as those quantities reach a plateau at that size. The dependence on γ is consistent with a global vertical temperature gradient, suggesting that possible edge effects due to the small size of the sample were negligible as compared to the effects of temperature gradient.

The density of 476 kg m−3 is high for dry snow, whose densities lie in the range 50 to 500 kg m−3. Akitaya (1974) noticed the occurrence of hard depth hoar for a snow density >300kg m−3 exposed to a temperature gradient of 39 Cm-1. Pfeffer and Mrugala (2002) observed formation of hard depth hoar for snow of 400 kg m−3 for a temperature gradient >20Cm-1 at a mean temperature of 7 C for 3 d. With a snow of density 230 kg m−3, they did not observe the formation of hard depth hoar, and concluded a critical density exists between those values. In the present experiment, the transition occurred in phases 2 and 3 where the temperature gradient was 52 Cm-1 with a mean temperature of about 2 and 4 C respectively. So, results must be representative of the process of formation of hard depth hoar. In addition, both Riche et al. (2013) and Calonne et al. (2017) report no evolution of the fabric for high-density snow, which would be consistent with the expectations that geometric-based competition between grains increases with density as the pore size decreases. However, the present results show that differences between grain orientations actually exist for growth and sublimation rates, even at high density. So this suggests that while the effects exist, they are less important in the selection of growing grains at that density and are more likely to occur at lower density.

5 Conclusions

For the first time, DCT monitoring of a temperature gradient metamorphism experiment was carried out. A snow sample of rounded grains (RG) with initial density of (476±5 %)kg m−3 has been subject to a temperature gradient of 17 Cm-1 for 1 d and then to a temperature gradient of 52 Cm-1 for 2 d at a mean temperature of about 2 C for 2 d and 4 C for 1 d. Several points may be concluded from this experiment.

  • The crystalline orientations have successfully been measured using DCT on about 900 growing and faceting grains.

  • The local flux has been computed at each point of the interface, and the obtained value is well consistent with previous studies.

  • Analyses show that the first mass transfers were overall higher on grains with the c axis close to the horizontal plane, for both sublimation and deposition, thus confirming the hypothesis of orientation selective metamorphism, even at high density.

  • Fluxes become independent of orientation at longer times for deposition while differences remain for sublimation. We attribute this effect to grain-to-grain interactions, particularly present at that density, preventing further growth of crystals.

Additionally, this study also opens new perspectives. Better accuracy on orientation, particularly at grains boundaries, can be obtained by a proper management of dynamic diffraction that occurs in ice, in the reconstruction process. Additionally, the observation of orientation selective grain growth and sublimation at high density gives confidence in the fact that this phenomenon actually occurs in more general conditions. The present experiment could be repeated in various conditions to characterize this effect in more detail. Finally, such data are necessary for models that describe the snow metamorphism using the crystalline nature of ice.

Code availability

Codes come from very different sources and are not all available publicly. Please contact the authors for further details.

Data availability

Research data are available on request. Please contact the authors.

Author contributions

RG and FF developed the cryogenic cell. RG, FF, WL and CG carried out the experiment. WL and IH post-processed the DCT volumes. RG, FF and CG analysed the data and wrote the article.

Competing interests

The authors declare that they have no conflict of interest.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank the reviewers, Ian Baker and Henning Löwe, for their insightful comments and suggestions. We are specifically grateful to Sabine Rolland du Roscoat, Fabien Lahoucine, Laurent Pézard, Jacques Roulle, Elodie Boller, Jean-Paul Valade and Paul Tafforeau for their contributions to the setup of the experiment. Thanks are also due to Armelle Philip, Philippe Lapalus, Alexis Burr, Anne Dufour, Isabel Peinke and Pascal Hagenmuller for their participation in the experiment monitoring. Finally, we would like to thank the ID19 beamline of the ESRF, where the images have been acquired.

CNRM/CEN is part of the LabEx OSUG@2020 (Investissements d’Avenir – grant agreement ANR10 LABX56). The 3SR Laboratory is part of the LabEx Tec 21 (Investissements d’Avenir – grant agreement ANR11 269 LABX0030).

Financial support

This research has been supported by the MiMESis-3D ANR project (grant no. ANR-19-CE01-0009).

Review statement

This paper was edited by Jürg Schweizer and reviewed by Ian Baker and Henning Loewe.


Adams, E. E. and Miller, D. A.: Ice crystals grown from vapor onto an orientated substrate: application to snow depth-hoar development and gas inclusions in lake ice, J. Glaciol., 49, 8–12,, 2003. a, b

Akitaya, E.: Studies on Depth Hoar, Contributions from the Institute of Low Temperature Science, 26, 1–67, 1974. a, b, c, d

Calonne, N.: Physique Des Métamorphoses de La Neige Sèche: De La Microstructure Aux Propriétés Macroscopiques, PhD Thesis, Grenoble, 2014. a

Calonne, N., Flin, F., Geindreau, C., Lesaffre, B., and Rolland du Roscoat, S.: Study of a temperature gradient metamorphism of snow from 3-D images: time evolution of microstructures, physical properties and their associated anisotropy, The Cryosphere, 8, 2255–2274,, 2014. a, b

Calonne, N., Flin, F., Lesaffre, B., Dufour, A., Roulle, J., Puglièse, P., Philip, A., Lahoucine, F., Geindreau, C., and Panel, J.-M.: CellDyM: A Room Temperature Operating Cryogenic Cell for the Dynamic Monitoring of Snow Metamorphism by Time-Lapse X-Ray Microtomography, Geophys. Res. Lett., 42, 3911–3918,, 2015. a, b, c

Calonne, N., Montagnat, M., Matzl, M., and Schneebeli, M.: The Layered Evolution of Fabric and Microstructure of Snow at Point Barnola, Central East Antarctica, Earth Planet. Sc. Lett., 460, 293–301,, 2017. a, b

Colbeck, S.: Theory of Metamorphism of Dry Snow, J. Geophys. Res.-Oceans, 88, 5475–5482,, 1983. a

de Quervain, M. R.: The Institute for Snow and Avalanche Research at Weissfluhjoch/Davos: The First Five Years (1943 to 1948), Ann. Glaciol., 4, 307–314,, 1983. a

Dumont, M., Flin, F., Malinka, A., Brissaud, O., Hagenmuller, P., Lapalus, P., Lesaffre, B., Dufour, A., Calonne, N., Rolland du Roscoat, S., and Ando, E.: Experimental and model-based investigation of the links between snow bidirectional reflectance and snow microstructure, The Cryosphere, 15, 3921–3948,, 2021. a

Fierz, C., Armstrong, R. L., Durand, Y., Etchevers, P., Greene, E., McClung, D. M., Nishimura, K., Satyawali, P. K., and Sokratov, S. A.: The international classification for seasonal snow on the ground, IHP-VII Technical Documents in Hydrology no. 83, IACS Contribution no.1, available at: images/0018/001864/186462e.pdf (last access: 16 August 2021), 2009. a

Fisher, N. I., Lewis, T., and Embleton, B. J. J.: Statistical Analysis of Spherical Data, Cambridge University Press, Cambridge, 1993. a

Flin, F. and Brzoska, J.-B.: The Temperature-Gradient Metamorphism of Snow: Vapour Diffusion Model and Application to Tomographic Images, Ann. Glaciol., 49, 17–21,, 2008. a

Flin, F., Brzoska, J.-B., Lesaffre, B., Coléou, C., and Pieritz, R. A.: Three-Dimensional Geometric Measurements of Snow Microstructural Evolution under Isothermal Conditions, Ann. Glaciol., 38, 39–44,, 2004. a

Flin, F., Brzoska, J.-B., Coeurjolly, D., Pieritz, R., Lesaffre, B., Coléou, C., Lamboley, P., Teytaud, O., Vignoles, G., and Delesse, J.-F.: An Adaptative Filtering Method to Evalute Normal Vectors and Surface Areas of 3D Objects. Application to Snow Images from X-Ray Tomography, IEEE T. Image Process., 14, 585–596,, 2005. a, b, c

Flin, F., Lesaffre, B., Dufour, A., Gillibert, L., Hasan, A., Rolland du Roscoat, S., Cabanes, S., and Pugliese, P.: On the Computations of Specific Surface Area and Specific Grain Contact Area from Snow 3D Images, edited by: Furukawa, Y., Hokkaido University Press, Sapporo, JP, Proceedings of the 12th International Conference on the Physics and Chemistry of Ice held at Sapporo, Japan, 5–10 September 2010, 321–328, 2011. a

Flin, F., Calonne, N., Lesaffre, B., Dufour, A., Philip, A., Rolland du Roscoat, S., and Geindreau, C.: The TG Metamorphism of Snow: toward an Evaluation of 3D Numerical Models by Time-Lapse Images Acquired under X-Ray Tomography?, Physics and Chemistry of Ice (PCI), Hanover, USA, 17—20 March 2014. a

Flin, F., Denis, R., Mehu, C., Calonne, N., Lesaffre, B., Dufour, A., Granger, R., Lapalus, P., Hagenmuller, P., Roulle, J., Rolland du Roscoat, S., Bretin, E., and Geindreau, C.: Isothermal Metamorphism of Snow: measurement of Interface Velocities and Phase-Field Modeling for a Better Understanding of the Involved Mechanisms, Physics and Chemistry of Ice (PCI), Zuerich, Switzerland, 8–12 January 2018. a

Furukawa, Y.: 25 – Snow and Ice Crystal Growth, Editor(s): Tatau Nishinaga, Handbook of Crystal Growth (Second Edition), Elsevier, United Kingdom, 1061–1112,, 2015. a, b

Haffar, I., Flin, F., Geindreau, C., Petillon, N., Gervais, P.-C., and Edery, V.: X-ray tomography for 3D analysis of ice particles in jet A-1 fuel, Powder Technol., 384, 200–210,, 2021. a

Hagenmuller, P., Chambon, G., Lesaffre, B., Flin, F., and Naaim, M.: Energy-Based Binary Segmentation of Snow Microtomographic Images, J. Glaciol., 59, 859–873,, 2013. a

Hammonds, K. and Baker, I.: Investigating the Thermophysical Properties of the Ice–Snow Interface under a Controlled Temperature Gradient Part II: Analysis, Cold Reg. Sci. Technol., 125, 12–20,, 2016. a

Hammonds, K., Lieb-Lappen, R., Baker, I., and Wang, X.: Investigating the Thermophysical Properties of the Ice–Snow Interface under a Controlled Temperature Gradient: Part I: Experiments & Observations, Cold Reg. Sci. Technol., 120, 157–167,, 2015. a

Krol, Q. and Löwe, H.: Analysis of Local Ice Crystal Growth in Snow, J. Glaciol., 62, 378–390,, 2016. a

Krol, Q. and Löwe, H.: Upscaling ice crystal growth dynamics in snow: Rigorous modeling and comparison to 4D X-ray tomography data, Acta Mater., 151, 478–487,, 2018. a

Kuroda, T. and Lacmann, R.: Growth Kinetics of Ice from the Vapour Phase and Its Growth Forms, J. Cryst. Growth, 56, 189–205,, 1982. a

Libbrecht, K. G.: The Physics of Snow Crystals, Rep. Prog. Phys., 68, 855–895,, 2005. a

Ludwig, W., King, A., Reischig, P., Herbig, M., Lauridsen, E., Schmidt, S., Proudhon, H., Forest, S., Cloetens, P., du Roscoat, S. R., Buffière, J., Marrow, T., and Poulsen, H.: New Opportunities for 3D Materials Science of Polycrystalline Materials at the Micrometre Lengthscale by Combined Use of X-Ray Diffraction and X-Ray Imaging, Mater. Sci. Eng. A-Struct., 524, 69–76,, 2009. a

Marbouty, D.: An Experimental Study of Temperature-Gradient Metamorphism, J. Glaciol., 26, 303–312,, 1980. a, b

Miller, D. and Adams, E.: A Microstructural Dry-Snow Metamorphism Model for Kinetic Crystal Growth, J. Glaciol., 55, 1003–1011,, 2009. a

Montagnat, M., Löwe, H., Calonne, N., Schneebeli, M., Matzl, M., and Jaggi, M.: On the birth of structural and crystallographic fabric signals in polar snow: A case study from the EastGRIP snowpack, Front. Earth Sci., 8, 365, 2020. a

Palenstijn, W. J., Batenburg, K. J., and Sijbers, J.: The ASTRA Tomography Toolbox, in: 13th International Conference on Computational and Mathematical Methods in Science and Engineering, CMMSE 2013, 4, 1139–1145, available at: (last access: 8 September 2021), 2013. a

Palenstijn, W. J., Bédorf, J., and Batenburg, K. J.: A distributed SIRT implementation for the ASTRA Toolbox, in: Proceedings of The 13th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine 2015, edited by: King, M., Glick, S., and Mueller, K., 166–169, available at: (last access: 8 September 2021), 2015. a

Pfeffer, W. T. and Mrugala, R.: Temperature Gradient and Initial Snow Density as Controlling Factors in the Formation and Structure of Hard Depth Hoar, J. Glaciol., 48, 485–494,, 2002. a, b

Pinzer, B. R., Schneebeli, M., and Kaempfer, T. U.: Vapor flux and recrystallization during dry snow metamorphism under a steady temperature gradient as observed by time-lapse micro-tomography, The Cryosphere, 6, 1141–1155,, 2012. a, b, c

Reischig, P., King, A., Nervo, L., Viganó, N., Guilhem, Y., Palenstijn, W. J., Batenburg, K. J., Preuss, M., and Ludwig, W.: Advances in X-Ray Diffraction Contrast Tomography: Flexibility in the Setup Geometry and Application to Multiphase Materials, J. Appl. Crystallogr., 46, 297–311,, 2013. a, b

Riche, F., Montagnat, M., and Schneebeli, M.: Evolution of Crystal Orientation in Snow during Temperature Gradient Metamorphism, J. Glaciol., 59, 47–55,, 2013. a, b, c

Rolland du Roscoat, S., King, A., Philip, A., Reischig, P., Ludwig, W., Flin, F., and Meyssonnier, J.: Analysis of Snow Microstructure by Means of X-Ray Diffraction Contrast Tomography, Adv. Eng. Mater., 13, 128–135,, 2011. a

Schneebeli, M. and Sokratov, S. A.: Tomography of Temperature Gradient Metamorphism of Snow and Associated Changes in Heat Conductivity, Hydrol. Process., 18, 3655–3665,, 2004.  a

Srivastava, P. K., Mahajan, P., Satyawali, P. K., and Kumar, V.: Observation of Temperature Gradient Metamorphism in Snow by X-Ray Computed Microtomography: Measurement of Microstructure Parameters and Simulation of Linear Elastic Properties, Ann. Glaciol., 51, 73–82,, 2010. a

Takahashi, Y. and Fujino, K.: Crystal Orientation of Fabrics in a Snow Pack, Low Temp. Sci., Ser. A, 34, 71–78, 1976. a

Wilson, C. J. L., Russell-Head, D. S., Kunze, K., and Viola, G.: The Analysis of Quartz C-Axis Fabrics Using a Modified Optical Microscope, J. Microsc., 227, 30–41,, 2007. a

Yokoyama, E. and Kuroda, T.: Pattern Formation in Growth of Snow Crystals Occurring in the Surface Kinetic Process and the Diffusion Process, Phys. Rev. A, 41, 2038,, 1990. a

Short summary
In this study on temperature gradient metamorphism in snow, we investigate the hypothesis that there exists a favourable crystal orientation relative to the temperature gradient. We measured crystallographic orientations of the grains and their microstructural evolution during metamorphism using in situ time-lapse diffraction contrast tomography. Faceted crystals appear during the evolution, and we observe higher sublimation–deposition rates for grains with their c axis in the horizontal plane.