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

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⁻³ subject to a temperature gradient of 52 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{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 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{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 (Colbeck, 1983). An important aspect of temperature gradient metamorphism is that the grain size increases while the density remains constant (Marbouty, 1980; Schneebeli and Sokratov, 2004; Srivastava et al., 2010; Pinzer et al., 2012; Calonne, 2014). Indeed, due to the local temperature field, some parts of the microstructure grow at the expense of others (Akitaya, 1974; 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 Kuroda, 1990; Libbrecht, 2005; Furukawa, 2015). 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 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{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 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{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 ≈450 kg m⁻³ 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 ≈7 h 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 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{1}}$ during phase 1, −2.09 ^∘C and 52 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{1}}$ during phase 2, and −4.07 ^∘C and 52 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{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×5 mm² 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 (≈10 min 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 $\mathit{\omega }+\mathrm{180}{}^{\circ }$ for the $\left(\overline{h}\overline{k}\overline{l}\right)$ 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=27 h 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 $\mathit{X}+\mathit{Y}+\mathit{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.1 h and t=78.6 h 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⁻¹. 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: $\mathit{\alpha }=(\mathit{c},\mathit{n})$, $\mathit{\beta }=(\mathit{z},\mathit{c})$ and $\mathit{\gamma }=(\mathit{z},\mathit{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 8 Density 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 476 kg m⁻³ (±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=917 kg m⁻³). 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⁻³.

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 %) m² kg⁻¹, the SSA decreases slightly during the first 30 h and then increases more and more rapidly to reach about (20±5 %) m² kg⁻¹.

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

$$\begin{array}{}\text{(1)}& {\mathit{a}}^{\left(\mathrm{2}\right)}={\displaystyle \frac{\mathrm{1}}{N_i}}\sum _{i=\mathrm{1}}^{N_i}{\mathit{c}}^i\otimes {\mathit{c}}^i,\end{array}$$

with cⁱ being the c-axis orientation of the ith voxel in the volume, ⊗ the tensor product and N_i the number of ice voxels. The orientation distribution is inscribed in an ellipsoid whose principal axes are given by the eigenvectors (v₁, v₂, v₃) of a⁽²⁾ and the length of the principal axis by its eigenvalues. Denoting the eigenvalues as a₁, a₂, a₃ with $\mathrm{1}>a_{\mathrm{1}}\ge a_{\mathrm{2}}\ge a_{\mathrm{3}}>\mathrm{0}$, we can identify three limit cases, defining three references of fabric.

• $a_{\mathrm{1}}>a_{\mathrm{2}}\simeq a_{\mathrm{3}}$: cluster-type fabric, majority of c axis in one direction (defined by v₁).

• $a_{\mathrm{1}}\simeq a_{\mathrm{2}}>a_{\mathrm{3}}$: girdle-type fabric, majority of c axis in a plane (defined by v₁ and v₂).

• $a_{\mathrm{1}}=a_{\mathrm{2}}=a_{\mathrm{3}}=\mathrm{1}/\mathrm{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 (a_{\mathrm{1}}/a_{\mathrm{2}})$, $\ln (a_{\mathrm{2}}/a_{\mathrm{3}})$ and $\ln (a_{\mathrm{1}}/a_{\mathrm{3}})$ 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 $\mathit{\alpha }=(\mathit{c},\mathit{n})$ and $\mathit{\gamma }=(\mathit{z},\mathit{n})$. The small red arrows represent the facets taken into account.

Figure 16 Interpretation of the diagram of mean flux against $\mathit{\beta }=(\mathit{z},\mathit{c})$ and $\mathit{\gamma }=(\mathit{z},\mathit{n})$. The small red arrows represent the facets taken into account.

From the 30 registered and segmented absorption images $(A_k{)}_{k\in \left[\left[\mathrm{1},\mathrm{30}\right]\right]}$, we computed local interface speeds v_k, where k is the index of the absorption image. At a given point of the surface of the volume $k\in \left[\left[,\mathrm{30}\right]\right]$ the normal n_k was computed using the ADGF algorithm (Flin et al., 2005). Denoting by $d_k^l$ the displacement of the interface along n_k from A_k to A_l, we computed both $d_k^{k+\mathrm{1}}$ and $d_k^{k-\mathrm{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 v_k at the instant t_k of acquisition of A_k was estimated using

$$\begin{array}{}\text{(2)}& v_k=\left\{\begin{array}{ll}\frac{d_k^{k+\mathrm{1}}}{t_{k+\mathrm{1}}-t_k}& \text{ if }k=\mathrm{1}\\ \frac{t_k-t_{k-\mathrm{1}}}{t_{k+\mathrm{1}}-t_{k-\mathrm{1}}}\frac{d_k^{k+\mathrm{1}}}{t_{k+\mathrm{1}}-t_k}-\frac{t_{k+\mathrm{1}}-t_k}{t_{k+\mathrm{1}}-t_{k-\mathrm{1}}}\frac{d_k^{k-\mathrm{1}}}{t_k-t_{k-\mathrm{1}}}& \text{ if }k\in \left[\left[\mathrm{2},\mathrm{29}\right]\right]\\ -\frac{d_k^{k-\mathrm{1}}}{t_k-t_{k-\mathrm{1}}}& \text{ if }k=\mathrm{30}.\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(3)}& {\mathit{\varphi }}_k={\mathit{\rho }}_i\cdot v_k.\end{array}$$

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 $\mathit{\alpha }=(\mathit{c},\mathit{n})$ and $\mathit{\gamma }=(\mathit{z},\mathit{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 ($\mathit{\alpha }=\mathrm{90}{}^{\circ }$), 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 ($\mathit{\gamma }=\mathrm{180}{}^{\circ }$); pixels at mid-height of the diagram correspond to surfaces whose normals are horizontal ($\mathit{\gamma }=\mathrm{90}{}^{\circ }$). A similar representation can be done by representing the mean flux against $\mathit{\beta }=(\mathit{z},\mathit{c})$ and $\mathit{\gamma }=(\mathit{z},\mathit{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 ($\mathit{\beta }=\mathrm{90}{}^{\circ }$), while the left of the diagram corresponds to grains with a vertical c axis (β=0).

Figure 17 Mean fluxes against $\mathit{\alpha }=(\mathit{c},\mathit{n})$ and $\mathit{\gamma }=(\mathit{z},\mathit{n})$.

Figure 18 Mean fluxes against $\mathit{\beta }=(\mathit{z},\mathit{c})$ and $\mathit{\gamma }=(\mathit{z},\mathit{n})$.

Figure 19 Mean fluxes against $\mathit{\alpha }=(\mathit{z},\mathit{c})$, for $\mathit{\gamma }\in [\mathrm{0},\mathrm{20}{}^{\circ }]$ (faces oriented upward), $\mathit{\gamma }\in [\mathrm{80},\mathrm{100}{}^{\circ }]$ (face oriented vertically) and $\mathit{\gamma }\in [\mathrm{160},\mathrm{180}{}^{\circ }]$ (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 γ: $[\mathrm{0},\mathrm{20}{}^{\circ }]$, $[\mathrm{80},\mathrm{100}{}^{\circ }]$ and $[\mathrm{160},\mathrm{180}{}^{\circ }]$, 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.1 h, from all figures, the bin-mean fluxes have modulus less than $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{7}}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ in all cases. At t=33.3 h, mean fluxes for $\mathit{\gamma }<{\mathrm{90}}^{\circ }$ 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 ($\mathit{\gamma }>\mathrm{90}{}^{\circ }$) 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öwe, 2018).

In addition, we can observe a second variation depending on the value of α. Higher rates of sublimation–deposition with absolute bin-mean fluxes of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ kg m² s⁻¹ are indeed found for $\mathit{\alpha }>{\mathrm{45}}^{\circ }$. For example, from Fig. 19, for α varying from 5 to 90^∘, one can see a trend for the fluxes varying from $\sim \mathrm{4}\times {\mathrm{10}}^{-\mathrm{7}}$ to $\sim \mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ kg m² s⁻¹ for γ in [160, 180^∘] and from $\sim -\mathrm{4}\times {\mathrm{10}}^{-\mathrm{7}}$ kg m² s⁻¹ to $\sim -\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ kg m² s⁻¹ 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.8 h, the same dependence on γ is observed, and the dependence on α has the same shape for $\mathit{\gamma }<{\mathrm{90}}^{\circ }$, while it is not visible for $\mathit{\gamma }>{\mathrm{90}}^{\circ }$. Overall, at all points of the diagram, the values are lower than at t=33.3 h. The evolution between t=33.3 h and t=78.8 h is regular so that times t=49.9 h, t=59.9 h and t=71.8 h 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.3 h, voxels belonging to crystals with a c axis in the horizontal plane or close to this plane ($\mathit{\beta }\in [\mathrm{45},\mathrm{90}{}^{\circ }]$) are more active.

Figure 20 Evolution inside a small subvolume of $\mathrm{520}\times \mathrm{520}\times \mathrm{520}$ µm³. 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.8 h, 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.1 h, the fluxes and variations on fluxes are small. The temperature gradient of 17 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{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⁻³. Indeed, it is known that high-density snow requires higher temperature gradients (Akitaya, 1974; Marbouty, 1980; Pfeffer and Mrugala, 2002) and longer times of evolution to get visible modifications of the microstructure. Under a 52 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{1}}$ temperature gradient, the maximum values of the bin-mean fluxes are of about $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{7}}$ kg m² s⁻¹. Pinzer et al. (2012) computed a macroscopic flux of about $\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{7}}$ kg m² s⁻¹ for snow evolving under 49 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{1}}$ temperature gradient at −3.4 ^∘C. Flin and Brzoska (2008) computed the highest fluxes of about $\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{7}}$ kg m² s⁻¹ under 16 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{1}}$. Thus, the present computed flux values are in very good agreement with this literature.

As described in Sect. 3.6, for t=33.3 h, we can see a trend for the most active points of the interface to be with $\mathit{\alpha }\ge \mathrm{45}{}^{\circ }$. 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 $\mathit{\alpha }=\mathrm{90}{}^{\circ }$. 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 $\mathit{\alpha }=\mathrm{45}{}^{\circ }$, 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 $\mathit{\alpha }=\mathrm{90}{}^{\circ }$. 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 Lacmann, 1982; Furukawa, 2015). Correspondingly, from Fig. 18, the most active faces correspond to crystals for which $\mathit{\beta }=\mathrm{90}{}^{\circ }$, 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 $\mathit{\alpha }\simeq \mathrm{45}{}^{\circ }$) disappear as the facets grow. Secondly, after the important growth of plates between t=25 h and t=33 h, 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.3 h, 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⁻³ is high for dry snow, whose densities lie in the range 50 to 500 kg m⁻³. Akitaya (1974) noticed the occurrence of hard depth hoar for a snow density >300 kg m⁻³ exposed to a temperature gradient of 39 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{1}}$. Pfeffer and Mrugala (2002) observed formation of hard depth hoar for snow of 400 kg m⁻³ for a temperature gradient >20 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{1}}$ at a mean temperature of −7 ^∘C for 3 d. With a snow of density 230 kg m⁻³, 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 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{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⁻³ has been subject to a temperature gradient of 17 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{1}}$ for 1 d and then to a temperature gradient of 52 ${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{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.