Vapor flux and recrystallization during dry snow metamorphism under a steady temperature gradient as observed by time-lapse micro-tomography

Dry snow metamorphism under an external temperature gradient is the most common type of recrystallization of snow on the ground. The changes in snow microstructure modify the physical properties of snow, and therefore an understanding of this process is essential for many disciplines, from modeling the effects of snow on climate to assessing avalanche risk. We directly imaged the microstructural changes in snow during temperature gradient metamorphism (TGM) under a constant gradient of 50 K m −1, using in situ time-lapse X-ray micro-tomography. This novel and non-destructive technique directly reveals the amount of ice that sublimates and is deposited during metamorphism, in addition to the exact locations of these phase changes. We calculated the average time that an ice volume stayed in place before it sublimated and found a characteristic residence time of 2–3 days. This means that most of the ice changes its phase from solid to vapor and back many times in a seasonal snowpack where similar temperature conditions can be found. Consistent with such a short timescale, we observed a mass turnover of up to 60 % of the total ice mass per day. The concept of hand-to-hand transport for the water vapor flux describes the observed changes very well. However, we did not find evidence for a macroscopic vapor diffusion enhancement. The picture of temperature gradient metamorphism produced by directly observing the changing microstructure sheds light on the micro-physical processes and could help to improve models that predict the physical properties of snow.


Introduction
Snow is a highly porous material consisting of a complex network of sintered ice crystals. Under environmental conditions, the temperature of alpine snow is usually very close to its melting point. In terms of the homologous temperature T H , i.e. the temperature relative to the melting point, snow under alpine conditions mostly occupies values of 0.9 < T H ≤ 1, corresponding to temperatures between −27 • C and 0 • C. The high homologous temperature, combined with the particularly high vapor pressure of water, can cause rapid sintering and large structural modifications by vapor diffusion. In particular, the naturally occurring temperature gradients in a snowpack cause a water vapor pressure gradient in the pore space with an associated vapor flux. This induces recrystallization and changes in morphology, and thus changes in the physical properties of snow. The structural changes are summarized by the term dry snow metamorphism under an external temperature gradient (Yosida, 1955;Sturm and Benson, 1997). More precisely, large temperature gradients in an alpine snowpack can be sustained for extended periods, especially near the ground. In this regime, fluctuations are often small compared to the magnitude of the gradient, and therefore we focus on temperature gradient metamorphism (TGM) under steady conditions in the following, neglecting complicated daily variations of the gradient near the surface (Pinzer and Schneebeli, 2009b).
It is the close relationship between physical properties of snow and its microstructure that promotes the interest in snow metamorphism. Mechanical properties are determined by the snow structure (Schneebeli, 2004) and depth hoar, a late stage morphological shape caused by temperature gradient metamorphism, is a key factor in the formation of avalanches (Schweizer et al., 2003). Considering optical properties, snow has a feedback on climate (Flanner and Zender, 2006) and the global energy budget (Dozier et al., 2009) by means of its large albedo (reflection coefficient). In addition, a snow layer can provide efficient thermal insulation and thus a decoupling between the ground and the atmosphere (Cook et al., 2007), a property that is related to the low thermal conductivity of snow Calonne et al., 2011;Shertzer and Adams, 2011). Atmospheric chemistry (Grannas et al., 2007) is influenced by the large internal surface area in snow, which provides a catalyst for heterogeneous photochemical reactions. Antarctic ice cores contain trapped air bubbles and show long-term changes in air content that are related to metamorphic processes in the snowpack before pore close-off (Raynaud et al., 2007). This change in air content seems to correlate with changes in insolation and possibly the formation of depth hoar (Hutterli al., 2009).
From an operational perspective, e.g. for avalanche forecasting, numerical snowpack models also have to consider metamorphic processes (Brun et al., 1989;Lehning, 2002). These models use empirical parameterizations based on the work of Marbouty (1980) and Baunach et al. (2001). Lehning (2002) treated grain growth and mass transport across layers as two separate processes, where part of the vapor bypasses the ice matrix and does not take part in the hand-to-hand vapor transport.
The close relation between snow microstructure and the physical properties of snow implies that an understanding of the underlying processes of snow metamorphism is necessary to improve the modeling efforts in the above mentioned fields. Despite significant efforts to understand, model, and parameterize the heat and mass transfer through snow, some details of the physics of snow metamorphism remain unclear. In particular, there has been a long standing controversy about the role of the ice structure in the macroscopic water vapor transport, which will be described in detail in Sect. 2.
The main reason for this controversy might be that controlled measurements in the field and in the laboratory are very challenging, often not repeatable, and mostly destructive. Laboratory experiments are easier to control, but traditional microscopy techniques are also destructive. Contrary to traditional microscopy, X-ray micro-tomography is a non-destructive technique which already has been applied recently for 3-D snow imaging (Coleou et al., 2001;Schneebeli and Sokratov, 2004;Kaempfer et al., 2005;Chen and Baker, 2010).
In the present work, we imaged the microstructure of snow samples in the laboratory during temperature gradient metamorphism. For this purpose, an X-ray desktop micro-CT was supplemented with a sample holder that allowed us to maintain a steady temperature gradient and to observe the structural changes inside an undisturbed snow sample undergoing temperature gradient metamorphism. This direct view of the structural evolution and mass transfer within snow leads to new insight into the ongoing processes. For the first time, we can directly observe where water molecules sublimate and where they are deposited. Moreover, using image analysis techniques and image-based numerical simulations of the temperature distribution, we can calculate the diffusive vapor flux, mass transfer rates, and the time scale for recrystallization in snow. In the following Sect. 2, we briefly review the controversial discussions about the macroscopic vapor flux, the key process of snow metamorphism. Section 3 contains a detailed description of our experiments and image processing procedures that were used to address the questions raised in Sect. 2. Results are presented and discussed in Sect. 4.

Current understanding of vapor transport during temperature gradient metamorphism
The underlying mechanism for TGM is vapor diffusion due to temperature-induced concentration gradients above neighboring ice surfaces (De Quervain, 1958, 1973Akitaya, 1974;Gubler, 1985;Colbeck, 1983). Yosida (1955) introduced the metaphor of "hand-to-hand transport" for the process of vapor transfer from one ice grain to another: each ice grain accumulates water molecules where its surface is colder than adjacent ice surfaces, and sublimates water molecules where its surface is warmer. The collective effect of many local "hand-to-hand" processes is a net vapor transport through a snow layer and an associated recrystallization of the ice structure. The transport of water vapor occurs via numerous local sublimation and deposition steps. This is different from the transport of non-reactive gases in porous media, which have to find their way around the ice structures through the pore space (Pinzer et al., 2010;Schwander et al., 1988). In that case, the net flux of the gas is determined by the diffusion coefficient of the gas and the tortuosity of the structure -a concept that was introduced to describe the increased path length through the pore space compared to the Euclidean distance between two points. Yosida (1955) investigated the vapor flux through snow experimentally by measuring the evolution of mass in control volumes separated by fine wire meshes along a vertical snow column. His experiments were inspired by the idea that ". . . since the ice particles present themselves as obstacles to the diffusion of water vapour, and, moreover, since it must be expected that they themselves give out or take in the diffusing water vapour . . . , the law of macroscopic diffusion may have to be expressed by some mathematical expression different in form from equation (1)" (Yosida, 1955).
With Eq. (1) Yosida is referring to the diffusion equation that describes the water vapor diffusion locally on the grain level, expressed here in terms of the particle concentration c(x, y, z, t) (m −3 ), where D 0 denotes the water vapor diffusion coefficient in air. By measuring the changes in mass of the separate control volumes over time while applying an external gradient, Yosida deduced an effective diffusion constant D eff for the macroscopic diffusion of water vapor through snow that was 4-5 times higher than D 0 . This seminal paper is still influencing the modeling of snow metamorphism today. A very similar setup was used by Satyawali (2000) to support the notion of diffusion enhancement, but his experimental data show only a very slight increase. Nevertheless, Satyawali claims that "the effective vapor-diffusion coefficient for snow is much greater than in air", and also derives a density dependence from a model. Sokratov and Maeno (2000) used the same principle, but instead of meshes to separate control volumes, they used the variations in density of the sieved snow to analyze the mass transfer for calculating the flux. Sokratov and Maeno (2000) conclude that the effective diffusivity is the same as in air, and stated that "No systematic dependencies of the effective diffusion coefficient were found on the snow density (porosity) or applied temperature gradient within these ranges." Giddings and LaChapelle (1962) presented a theoretical model based on vapor transport in porous media. Their conclusion is that "in practice, D eff should be a slight amount less than D 0 . . . ." However, the field experiment of Giddings and LaChapelle could not elucidate this question, and in fact they had to adjust the equation used to model grain growth with a free parameter. This is not surprising in hindsight, because instead of measuring the water vapor flux, they measured the growth in grain size. This result seemed to question the theory, and so the paper was rarely cited.
In contrast, Colbeck (1993) developed an analytical model based on a simplified lattice geometry. He concludes " . . . that D eff in snow is about five times greater than its value in air and decreases slowly as the average pore size increases." He supports this work with earlier theoretical studies by Colbeck (1983), Gubler (1985), and Sommerfeld (1983), which predicted an enhanced vapor flux. Experimental results by Voitkovskii et al. (1988), which indicate that the diffusion coefficient in snow is similar to the one in air, were discussed by Colbeck, but were considered to suffer from experimental problems.
It is obvious that the measurements of the water vapor diffusion coefficient in snow are not consistent, and that the issue of diffusion enhancement is discussed controversially in the literature. As Yosida (1955) already has discussed, the wire meshes can have a detrimental effect on the measured flux, because as soon as a gap forms, an artificially high temperature gradient is introduced across the gap, biasing the measured vapor flux.
Based on above models and data, the notion of diffusion enhancement was generally accepted in the snow community. Miller et al. (2003) and Miller and Adams (2009) extended the approach introduced by Colbeck (1993) by using a finite difference approximation and simulating observed grain-growth rates with some parametrization. In their model, continuous growth of an individual crystal is assumed. Christon et al. (1994) conducted a thorough threedimensional numerical simulation of temperature-gradient metamorphism using finite elements. The underlying geometry was inspired by the three grain morphologies introduced by Kry (1975). Christon et al. created an idealized periodic grain-cell geometry containing chains, links, and grains. The simulation of one idealized cell was quasi-steady, sublimation and deposition were not considered. They concluded that the diffusivity was enhanced, but with an upper limit of approximately two. Kaempfer and Plapp (2009) improved this approach by using a phase-field model, which incorporates sublimation and deposition. They showed qualitatively in their 2-D simulations that ice is rapidly sublimated and deposited without a necessarily large change in the shape of the structure and therefore of the crystals. Christon et al. (1994) and Kaempfer and Plapp (2009) both showed that, for conditions of interest here, curvature effects are negligible compared to temperature gradient effects. Sturm and Benson (1997) were not directly interested in the effective vapor diffusion coefficient but more in the metamorphic changes, exemplified in the growth of crystals and their shape, and the macroscopic water vapor transport. They measured in the field the evolution of a subarctic snowpack under large temperature gradients. With their careful measurements of grain growth (deduced from grain size evaluations obtained by sieving), grain shape and mass redistribution across different layers derived from δD and δ 18 O fractionation, they derived an inter-particle flux and a layer-tolayer flux. They concluded that "Calculated layer-to-layer vapor fluxes are ten times higher than inter-particle fluxes, which implies that depth-hoar grain growth is limited by factors other than the vapor supply." They observed that ". . . gain and loss of water molecules due to sublimation from grains takes place at a rate many times higher than the rate at which grains grow", which implies that net vapor flux through the snowpack and grain growth are coupled processes, but grains do not grow at the same rate as the net vapor flux suggests.
However, there are open questions in the context of the macroscopic water vapor transport. The key questions are: (1) What is the effect of the tortuosity of the ice structure on the macroscopic vapor transport? (2) How large is the contribution of the latent heat to the total heat flux? (3) How long does an individual water molecule reside inside a control volume?
A direct observation of the mass redistribution process during temperature gradient metamorphism could address these questions, in particular whether the diffusion coefficient is influenced by the ice geometry. The subject of this study is to present experiments that shed light on some of these questions, based on microscopic observations combined with numerical simulations.

Methods
We observed the changes of the snow structure using Xray microtomography (CT) in time-lapse mode. The sample holder was temperature controlled and a temperature gradient was set along the vertical z-axis (Schneebeli and Sokratov, 2004). The time-lapse CT images were used to extract the characteristic residence time of ice and to determine the water vapor flux experimentally, as described in Sects. 3.3.2 and 3.3.3. The following paragraphs provide details about the experimental protocol and the image processing.

Snow samples
We measured three time series with different snow samples subjected to temperature gradients of approximately 50 K m −1 . The snow was collected in the field. Series 1 was sieved (1.4 mm sieve) into the sample holder and left sintering at constant temperature of −15 • C for one week. The sieving procedure is physically very similar to a redeposition of snow by wind. Series 2 and 3 were undisturbed samples extracted from a homogenous layer in a natural snowpack. The three samples covered a volumetric density range from 0.28 (257 kg m −3 ) to 0.34 (312 kg m −3 ). The experimental conditions are given in Table 1. An ice layer of 4 mm thickness was frozen to the lower heat flux plate, such that no gap could form between the plate and the snow during heating by sublimation. The contact was checked before and after each experiment by a scan over the full sample height.

Data acquisition
The sample was heated at the bottom and passively cooled at the top, as described by Schneebeli and Sokratov (2004). The samples had a diameter of 51 mm and a height between 18 and 21.3 mm, depending on the thickness of the ice plate that created a well defined boundary condition. A µCT80 (Scanco Medical AG) computed X-ray micro-tomograph (CT) was programmed to scan the samples every 8 h (Series 1 and 2) or every 4 h (Series 3). The scan time for acquiring 1000 projections per 180 degrees was approximately 2 h in stacks of 104 slices. Because the scan was always from the same position, the slices had equal spacing in time. The voxel size was 25 µm for Series 1, and 18 µm for Series 2 and 3 (one voxel corresponds to a three-dimensional pixel). The reconstructed CT images were filtered with a Gaussian filter (support 2 voxels, standard deviation 1 voxel) and the threshold for segmentation was determined by matching the volumetric density of the binary images to the gravimetric density of the sample. Using the minimum between the peaks of air and ice in a grayscale histogram gives the same threshold. A subvolume of 300 × 300 × 198 voxels was evaluated for Series 1 (corresponding to 278 mm 3 ) and Series 2 (corresponding to 103 mm 3 ), and a subvolume of 504 × 504 × 198 voxels (293 mm 3 ) for Series 3.
Although the volume of snow in the sample holder was comparatively small (about 40 000 mm 3 ), the total sample volume was about 200-times larger than the representative elementary volume -defined as the minimum volume that is representative for bulk properties of different snow typeswhich is less than 6 3 mm 3 (Kaempfer et al., 2005;Brzoska et al., 2008;Zermatten et al., 2011;Calonne et al., 2011).

Quantitative data analysis
The three-dimensional images show the sintered ice structures and the pores. Together with the time-axis they form a four-dimensional dataset. The short time spacing between consecutive images allowed us to follow visually and to quantify numerically the amount of sublimation and deposition on each grain (ice structure). We used this information to measure vapor mass flux using two complementary algorithms. The first algorithm works directly on the images and uses particle image velocimetry (PIV). The second algorithm uses direct numerical finite-element simulation at one timestep to calculate the temperature and vapor pressure field in the ice and pore structure, and by these data the vapor flux. The algorithm is explained in Sect. 3.3.3. The dynamics of recrystallization on single grains was quantified by the ice residence time, which was calculated based on the measured age of each voxel. Residence time was never measured before, because it cannot be observed by optical methods without disturbance. In the following sections we develop the necessary methods.

Morphometric data
Structural parameters of the segmented ice structure were extracted with the software tools of the CT device (Image Processing Language, Scanco Medical). Volumetric density was determined by counting the ice voxels and dividing by the total number of voxels present in the investigated subvolume, while specific surface area (SSA) was measured by triangulation of the ice-air interface (Kerbrat et al., 2008). The structural parameter of most relevance for the present study is the number of ice structure intersections per unit length, S.N, which was calculated from the segmented image as the inverse of the mean distance between the skeletonized ice structures (Hildebrand et al., 1999).

Ice residence time and ice turnover rate
The continuous sublimation and recrystallization of the ice structures within the snow allow for the possibility of a complete turnover of the ice, meaning that all molecules of a grain (an ice crystal) are replaced by freshly deposited Table 1. Experimental conditions and structural evolution of the snow. T : mean temperature of sample ( • C), ∇T : temperature gradient (K m −1 ), t: period between CT-scans (h), duration: total time of experiments (h), n: solid volume fraction (%) (Fig. S1), SSA i , SSA f : specific surface area at the beginning and end of the experiments (mm −1 ) (Fig. S2 molecules. In two consecutive 3-D-images, freshly deposited ice mass is detected by the appearance of an ice voxel at a site where there was none in the image before and, accordingly, sublimating ice mass manifests itself in the disappearance of a previously existing ice voxel. In principle, such an effect could also be caused by settling of the sample, as observed by Theile et al. (2011). Settling, however, would increase the density of the snow and would therefore be detected in the temporal evolution of the snow density (see also the discussion in Sect. 4). A basic argument to rule out settling effects can be obtained by calculating the expected rate of settling. Using the mean value of strain rate between equitemperature and temperature-gradient snow (Armstrong, 1980), 0.7 × 10 −9 s −1 , the rate of settling is less than 3.2 µm d −1 for a pressure of 70 Pa, the maximum at the bottom of the snow sample.
We define the residence time T R of ice as the current age of an ice voxel, i.e. T R is the time that has elapsed since the first appearance of the voxel (Fig. 1). To calculate the residence time, we stepped through the time series of images and incremented a time label for each unmodified ice voxel by t. A new voxel was assigned the time t = t, and was therefore affected with an uncertainty of less than t (because it might be younger). This is the only source of error in this method as long as the subsequent images are correctly aligned. Any misalignment of one image in the series would lead to a peak of young voxels and distort the distribution of residence times towards the lower end, which was not observed in our data. After an initial transient regime where the maximum age of the voxels increases with each step, the distribution of residence times is expected to converge towards a steady-state distribution when most of the ice voxels have been turned over at least once. An illustration of the algorithm can be found in Fig. 1.
The volumetric turnover rate R(kg m −3 s −1 ) is defined as where m is the relocated ice mass between two images separated by t, and V is the volume under consideration. By counting the new and disappearing voxels, respectively, the sublimating and depositing mass can be calculated. A difference in sublimation and deposition rates would imply a change in snow density over time. We will return to this point in Sect. 4. The zoomed boxes show the same ice structure (black) at two subsequent time-steps. Subtracting the images yields the sublimated (red) and freshly deposited (green) ice mass. New voxels are assigned an "age" of t = 8 h; the time label for unmodified voxels (gray) is incremented by t. The color-coded image on the right shows a steady-state distribution of residence times T R . Note that the bottom of that image corresponds to the warm side, as indicated by the temperature gradient ∇T .

Vapor mass flux
Yosida's term "hand-to-hand" vapor transport vividly describes the underlying process that is responsible for the morphological changes in the snowpack. The net effect of the hand-to-hand transport from individual snow structures to neighboring structures is a water vapor flux in the pore space in the opposite direction of the temperature gradient. On the microscopic level, vapor flux is governed by the diffusion equation (Eq. 1), together with Fick's law. On a macroscopic level, we still used Fick's law, but the diffusion coefficient was replaced with an effective diffusion constant D eff (m 2 s −1 ): Here, M is the mass of one molecule (kg), and the mass flux j is given in (kg m −2 s −1 ). Estimating the size of this net vapor flux -and thus the influence of the ice structure on the flux -was the goal of the following data analysis. We developed two independent www.the-cryosphere.net/6/1141/2012/ The Cryosphere, 6, 1141-1155, 2012 methods for deducing the vapor flux from measured CT data: the first one is based on the structural changes of the ice matrix between two different time points. The second method numerically solves for the water vapor concentration in the three-dimensional volume of air and ice with a prescribed temperature boundary condition.

Vapor mass flux obtained by particle image velocimetry (PIV)
Because ice structures gain mass at the (colder) bottom and lose mass at the (warmer) top, vapor transport results in an apparent displacement of the structures. The basic idea of this technique is to estimate the vapor flux through the pore space by locally evaluating the apparent displacement of the visible ice structures, using spatial correlations between subsequent CT images. Spatial correlations of similar but displaced images can be measured with particle image velocimetry (PIV), developed in the field of fluid dynamics. It is used to track the motion of tracer particles in fluids (Westerweel, 1997). Here, for every CT image, a displacement vector = ( x, y, z) with respect to a subsequent image was assigned to each ice voxel by shifting one image until the local correlation within a three-dimensional window was maximized. It is important to note that the correlation was maximized by displacing one image in three dimensions, thus obtaining a three-dimensional vector field. The window size and the temporal separation between the scans influence the signal-to-noise ratio of the calculated displacement field, and deserve closer attention. Choosing the window size is a trade-off between spatial resolution and smoothness of the displacement vector field. If the 3-D window is much smaller than a typical structural element, there is no distinct maximum of the correlation, and no unique displacement vector is obtained. On the other hand, a window too large smoothes local variations and finally underestimates the displacement. Based on convergence tests with different window sizes, the optimal window size we chose was 32 voxels, ensuring a well defined maximum in the correlation. The choice of temporal separation t of scans takes into account the changing shapes of ice structures and is a trade-off between strong signal and high spatial correlations. The signal, i.e. the displacement field, increases when images separated by two or three times t are used instead of subsequent CT scans. However, the shape of the structures changes over time, so the correlation decreases. We chose 16 h difference for Series 1 and 2, and 8 h for Series 3, corresponding to every second CT-scan (2 t separation) for each series. As a test of the implementation we chose a typical CT image and translated it manually by several pixels in each direction, which could be well recovered by the PIV algorithm, with an error below 1 %.
The displacement vector field encodes information about the apparent movement of the ice structure. This apparent movement is caused by the invisible vapor flux in the opposite direction, yet both mass fluxes have the same magnitude. Fig. 2. Simplified geometry to illustrate the estimation of vapor flux from the apparent movement of the structure. The wiggly arrows show the schematic vapor flux, orange caps indicate the locations of sublimation. The reference area A proj is used for the derivation using the Wigner-Seitz radius. The displacement vector d, the evaluated slice P and the structural length l are also schematically indicated. The temperature gradient is pointing along the z-axis, corresponding to higher temperatures at the bottom. This definition of the coordinate system, with the x-axis pointing into the paper plane, is used throughout this study.
In other words, we use the apparent mass flux of ice, measured by PIV, to infer the water vapor mass flux in the opposite direction. Assuming an isotropic and homogeneous distribution of ice structures, we will give a short derivation of the equality of the two mass fluxes. In a simple screening approach, which is a first order approximation, each ice structure exchanges vapor only with its nearest neighbors (Fig. 2). It follows that every plane P perpendicular to the temperature gradient "sees" water molecules crossing from the lower side of the plane to the upper side which originate from at most a depth l, the typical structural separation (or pore diameter). For simplicity, we assume spherical particles and discuss the effect of dropping this condition later. The aim is to obtain an expression for the volume (and thus mass) of ice that has to pass the plane P when the structure moves down by a distance d. The projected area that contributes to the flux through the plane can be calculated by taking the mean separation of the structures to be the Wigner-Seitz radius, l = (3/(4π n)) 1/3 , where n is the particle density. Then the projected area contributing to the flux through plane P is given by The above approximation, i.e. the right hand side of Eq. (4), is motivated by simplicity; it results in a 12 % lower area for a typical volume density of = 0.3 and for a volume density of = 0.42, the exact expression and the approximation yield the same value.
Finally, we approximate the volume of the sublimating spherical caps as d A proj , where d is the vertical displacement of the structures. The error for this final approximation is opposite to the error from Eq. (4), and is small if the displacement is small compared to the radius of the sphere; for d = r/5, the relative difference between exact cap volume and the approximation is 0.3 %. With this model, the flux through a plane j z is obtained as which is exactly the mass flux of a downward moving ice structure. The underlying model that led to this expression is probably oversimplified, neglecting vapor exchange between non-neighboring particles and non-spherical shapes, and the derivation is accurate by about 10 % for a typical density of = 0.3. However, it illustrates the concept of how to deduce the mass flux from the apparent structural movement. Abandoning the sphere model, which has very little to do with real snow structures as shown by CT-images (Fig. 3), increases the projected area but does not affect the flux. The average overall slices is defined as the flux j PIV .
A stronger argument for the validity of Eq. (5) is obtained a-posteriori: as will be shown later in the results section, the complete turnover of almost all the ice mass takes place in a rather short time (about two days). At that point, the structure has apparently "flown down" by such a large amount that no original ice structure remains. In that case, Eq. (5) can be interpreted as the apparent visible mass flux of the ice particles moving downward through a plane perpendicular to the temperature gradient. The invisible and oppositely oriented vapor flux must have the same magnitude.
A more severe concern about the reliability of this method is the strong sensitivity to positioning errors in the CT device. If there is a positioning error on the order of a voxel size, then the displacement fields suffer from random fluctuations. In the current experiments, no feedback mechanism for correction of such errors was available, and therefore fluctuations from positioning errors were possibly present in the results.

Vapor mass flux obtained by finite element simulation (FE method)
The vapor flux is driven by concentration gradients between neighboring structures, which are induced by different temperatures at the surfaces of those structures. Given the spatial distribution of ice in 3-D, the Laplace equation for the steadystate temperature distribution can be solved numerically, and thus the concentration boundary conditions at the ice-air interfaces can be obtained supposing local equilibrium. For the calculation of the equilibrium vapor pressure, the interfaces are considered flat, because -in the presence of a temperature gradient of the order of the one in the present study - Fig. 3. A slice of a displacement vector field calculated in three dimensions. Light grey represents sublimated voxels and dark grey freshly deposited voxels; mid-gray indicates no change with respect to the previous image. The reference system and the scale for the displacement vectors are shown on top, the temperature gradient was pointing downward. The physical dimensions of the presented slice are 2.5 × 1.9 mm 2 . curvature effects (Kelvin effect) only play a role for radii of curvature below 10 −7 m (Christon et al., 1994;Flanner and Zender, 2006). Solving for the concentration field subject to these Dirichlet boundary conditions yields the mass flux.
Two assumptions are the basis for the described procedure. First, the calculation assumes that at every moment the diffusion field is determined by the temperature field. However, the vapor flux has a feedback on the temperature field, both because the structures (and thus the pathways through the snow) are changed and latent heat is transported by the sublimation-deposition process, thereby heating deposition sites. The structural changes can be neglected because the time scales for temperature-equilibration (order of minutes) and structural modifications (order of many hours) are very different. Note that growth rates during snow metamorphism are relatively slow (typically on the order of 100 µm d −1 in our experiments). Also, the local influence of the released latent heat can be neglected, since the ice volume where the vapor is deposited is connected to the highly conductive ice network and the heat is immediately conducted away. This is in contrast to the growth of isolated ice crystals in the atmosphere, where the latent heat generated by solidification has to be diffused through the air. In that case, the heat increase has a direct effect on crystal growth.
The second assumption is that the vapor flux is completely diffusion limited, similar to the simulations by Flin and Brzoska (2008). The attachment kinetics of water molecules during ice crystal growth are not yet well understood and in general there is a large gap between crystal growth theory and the phenomenological observations of real crystal www.the-cryosphere.net/6/1141/2012/ The Cryosphere, 6, 1141-1155, 2012 growth (Nelson, 2001). The difficulty in determining the deposition coefficient α from measurements of the growth velocity in crystal growth experiments lies in the interplay between diffusion (of heat and molecules) and the attachment kinetics. Although the experiments on crystal growth are not in reasonable agreement with each other, there is a general consensus that for molecularly rough surfaces α ≈ 1, while on facets α < 1. The impact of the deposition coefficient on the diffusion field depends probably on how many facets there are present in the snow. By inspection of the CT images, the fraction of faceted surfaces with respect to the total surface area was small, typically lower than 20 %. To get an estimate for the vapor flux, we assumed α = 1 on all interfaces.
Releasing either assumption has the effect of slowing down the diffusion. If the latent heat generated by solidification cannot be diffused away, it will decrease the local temperature gradient across the pore. If the surface kinetics lead to a buildup of a supersaturation over the faceted surface, the concentration gradient will be flattened and the flux will be locally decreased. In this sense, our calculation gives an upper limit for the flux, which will be useful for discussing general concepts like diffusion enhancement or contribution of vapor flux to heat transport through snow.
We implemented the finite-element procedure by converting all voxels (ice and air) to eight-node brick finite elements and solving numerically for the temperature distribution in the ice and pore network, taking into account the respective material properties of ice and air (Kaempfer et al., 2005). The thermal conductivities for ice and air were set to 2.34 W K −1 m −1 and 0.024 W K −1 m −1 , corresponding to values at −15 • C. Between −15 • C and 0 • C, the thermal conductivity for air rises by 4 %, while the one for ice drops by 5 % (Slack, 1980;, resulting in an uncertainty of a few percent when compared to the experimental situation. We set the boundary conditions as they were in the experiments with insulation at the sidewalls (homogeneous von Neumann) and a vertical temperature gradient defined by fixed temperatures at the bottom and top nodes (Dirichlet).
To obtain the vapor concentration field c(x, y, z) in the pore space from the temperature distribution, a second Laplace equation with Dirichlet boundary conditions dictated by the temperature field has to be solved. With diffusion limited growth, the boundary condition consists in setting the vapor pressure above the ice surfaces to the equilibrium vapor pressure . Curvature effects on the equilibrium vapor pressure can be neglected since they are small compared to effects due to temperature differences as considered here (Kaempfer and Plapp, 2009). A further observation allowed us to linearize these equations: Since the temperature difference across a pore is very small (assuming a large pore of 1.4 mm diameter and an average temperature gradient of 50 K m −1 , the temperature difference is 0.07 • C), the relative deviation from linearity of Fig. 4. Visualization of (a) the structure for one slice of Series 2, at t = 320 h, (b) temperature distribution (shown in mK w.r.t. the boundary conditions; real temperature at the bottom was −7.05 • C), and (c) temperature gradient (K m −1 ). Note that the temperature boundary conditions were equivalent to a gradient of 50 K m −1 . The coordinate system coincides with the one in Fig. 2, as well as the direction of the macroscopic temperature gradient. The physical dimensions of the slices are 3.6 × 5.4 mm 2 .
the equilibrium vapor pressure is of the order of 5 × 10 −6 . Therefore, using c(x, y, z) = c e (T (x, y, z)), where c e is the equilibrium water vapor concentration at saturation over ice (Koop et al., 2000), taken at the local ice surface temperature, fulfills the Laplace equation to a very good approximation. Using the equilibrium vapor pressure p e and the ideal gas law, the concentration gradient in z-direction is then obtained by Fig. 5. Structural evolution of Series 2, observed in a subvolume (0.9 × 3.6 × 3.6 mm 3 ) of the measured field of view. The temperature gradient is pointing downward along the z-axis, as defined in Fig. 2 The upper row shows the structural changes during the first 72 h (every 24 h), the lower row shows the same subvolume during the last 144 h of the experiment (every 48 h). Due to the non-destructive nature of the experimental technique, we can observe a cup crystal growing in an undisturbed environment. The non-destructive character of the measurements allows for a direct observation of deposition sites of the water molecules. A video showing the complete time series consisting of 83 tomographic measurements is provided in the Supplement.
where k B is Boltzmann's constant, T is the temperature, and the subscript z for the gradients denotes their projection onto the z-axis. The mass flux follows from Fick's law, j z = −D∇c z , where ∇c z is the mean concentration gradient along z in a slice perpendicular to the macroscopic temperature gradient. In Fig. 4, an exemplary gradient field is shown for the structure at time point t = 320 h. Finally, j FE is the average vapor flux over all slices.

Results and discussion
Three time series of metamorphosing snow have been obtained with snow parameters as described in Table 1. The number of individual CT scans for Series 1, 2, and 3 was 48, 83, and 66, respectively.

Visualization of the structural evolution
The images obtained by time-lapse tomography give a general impression of how the structure evolves over time. As was already observed by Schneebeli and Sokratov (2004), the density did not change significantly over time; numerical values of the ice fraction are given in Table 1, together with the standard deviation of the temporal fluctuations. For the interested reader, a graph with the temporal evolution is shown in the Supplement (Fig. S1), along with the evolution of specific surface area (Fig. S2). While the density stays constant over time, a significant coarsening of the structures could be observed. The three series that we measured compare to the "low density snow" of the work by Schneebeli and Sokratov (2004), and we observed the same general trend towards coarser structures. A visualization of the ice matrix at different time points reveals the dramatic structural changes that occur during temperature gradient metamorphism (Fig 5). The great potential of the CT method lies in the in situ capabilities: the structure is imaged without disturbing the sample. This allows for much more than just a general comparison of geometrical parameters. As Fig. 5 vividly illustrates, we can follow the process of sublimation and deposition at the single grain level. For the first time, the growth of a cup crystal in an undisturbed environment can be directly observed. The lower row of Fig. 5 shows a close-up (3.6 × 0.9 × 3.6 mm 3 ) of the field of view that was scanned, with one image every 48 h. The five intermediate images (since t = 8 h between two CT scans) are not shown to illustrate the development of the cup crystal. Movie S1 in the Supplement shows the complete evolution of Series 2. Only the movie shows the reader the full dynamics and a complete visualization of the recrystallization. Interesting timesteps are after 48 h, where the predominantly rounded grains become clearly faceted. A further interesting feature occurs around 500 h. On a larger facet an instability forms in the form of a small ridge, which acts as one seed for the development of a cup crystal (the cut through the cup crystal is created by the visualization). Interestingly, the inside of the cup crystal is rather rough. In the following we will not try to explain why facets and cup crystals appear, beyond the recognition that a temperature gradient alone is not sufficient, especially in the first few days. This implies that the size of the pores may be important for formation of cup crystals as well.

Ice residence time and mass turnover
With the algorithm described in Sect. 3.3.2, the residence time T R of ice voxels was determined. Because all ice voxels start initially with an age of zero, it takes one complete recrystallization cycle before all voxels have been assigned a meaningful residence time. After an initial transient phase of three to five days, the distributions of residence times reached a steady state. The time it took for 98 % of the voxels to be relocated was 136 h for Series 1, 88 h for Series 2, and 103 h for Series 3. Figure 6 shows structures of Series 2 with color coded residence times for three different time points. Since the structures grew at the lower side of the grains, the youngest voxels were found there. It was striking that even for larger structures like the cup crystal, the oldest voxels at the top were mostly younger than 200 h. Indeed, according to the visual impression, most of the voxels of this crystal have an age less than 100 h.
The age distribution after the transient phase was fitted by the exponential function resulting in e-folding times τ on the order of 2 to 3 days (Fig. 7). For example, we observed an increase of τ from 50 to about 70 h during the final 441 h (18 days) of Series 2. To state this in a different way, consider that after a time τ , 63 % of the ice voxels have changed phase at least once. With τ ≈ 70 h, there are ice voxels within the snow that have changed phase about 10 times after the 665 h of the experiment. The slight increase of residence time can be explained by the slow increase of the structural size during temperature gradient metamorphism. A larger structural size at equal vapor flux will cause a longer residence time, because a larger structure takes longer to sublimate under otherwise identical conditions.
Under which circumstances a particular structural element is the seed for a larger structure is an interesting question.  6. Visualization of residence time with data from Series 2, of which the structure is shown in Fig. 5. The bottom of the structures are deposition sites for vapor, therefore the youngest ice voxels can be found there. The topmost parts, where sublimation takes place, host the oldest voxels. Note that even for the largest structures, the residence time of voxels rarely reaches 200 h. The direction of the temperature gradient (which coincides with z-axis) is shown in Fig. 2.
Writing about growing grains as a conserved entity, as is often done in the literature, seems confusing to us since after 2 to 3 days almost nothing of the original grain remains. Only the "memory" of the grain, encoded in the temporal correlation of the structure, survives. A possible description of the coarsening mechanism could be based on the slowly  Fig. 7. Probability distribution of ice residence times at the same three time points shown in Fig. 6. Symbols represent the measured data, while the curves are fitted exponential functions (Eq. 7, correlation coefficient r 2 > 0.95 for all series). The e-folding time τ , shown as a vertical line for each curve, slowly increases over time, but is always only a small fraction of the total experiment duration. fluctuating temperature field: whether a structural element grows or shrinks depends only on the difference in growth rate and sublimation rate. The continuous re-ordering of the ice structure produces fluctuations of the temperature gradients in the pore space on the time scale of many hours. The probability that larger structures survive such fluctuations is higher, while smaller structures are more likely to disappear. Time-lapse tomographic data open the possibility for the analysis of this hypothesis, but is beyond the scope of the current paper.
Counting the sublimated and freshly deposited voxels gives an estimate of the rate at which ice mass (in kg) is relocated per unit volume. We calculated the turnover rates R (Eq. 2) by averaging the sublimated and deposited mass for each series (Fig. 8). Both deposited and sublimated mass were equal on average, as was to be expected by the absence of significant changes in density (Fig. S1). A systematic difference of both quantities would necessarily lead to a net increase or decrease of density. While the turnover rate decreased for Series 1 and 2, it slightly increased for Series 3. Referring to the time scale of one day, the amount of relocated mass was very high. The average daily turnover amounted to 121, 130, and 190 kg m −3 d −1 for Series 1, 2, and 3, corresponding to 47, 43, and 61 % of the total ice mass. We think that R is important for chemical processes, as it influences the release and burial of chemicals (Domine and Rauzy, 2004).
Recalling the hand-to-hand mechanism of vapor transport, a connection between structural coarseness and the turnover rate was expected. In particular, if more small structures are involved in passing a certain flux through the snow volume, then the volumetric turnover rate is expected to be larger. On the other hand, fewer and larger structures would result in less total ice mass that is moved in a certain volume. Therefore, we calculated the structure number S.N (see Sect. 3.3.1) for each of the three series (Fig. 9). Again, Series 3 exhibits a qualitatively different behavior: while S.N decreases for Series 1 and 2, it increases for Series 3, which had a higher density and was exposed to a higher mean temperature. We suspect that this sample was closer to a snow type usually described as "hard depth hoar" (Perla, 1985). In that case it was described empirically that the crystals cannot grow freely and are hindered in the growth by their neighbors. The microscopical explanation is not yet clear.
Plotting the turnover rate versus S.N shows indeed a correlation ( Fig. 10) for Series 1 and 2. Series 3 covered a range of S.N that was too small to reliably fit a linear function to it. A correlation between a structural parameter like S.N (based on structure) and the turnover rate (requires time-lapse data) could be very useful for predicting the latter given just the structure and the thermal boundary conditions. However, extensive data are needed to establish the functional form of such a correlation.

Vapor flux and effective diffusion constant
We calculated the vapor flux using both the PIV method and the FE method, as detailed in Sect. 3.3.3. An example vector field for the PIV method is shown in Fig. 3. For an example of a resulting temperature and gradient distribution obtained by the FE method see Fig. 4.
In addition, we compared these fluxes to a simple continuum model where all the snow is replaced by air and the flux is calculated between two ice plates. In such a continuum model, the flux is given by  Fig. 10. Correlation between the turnover rate and the geometrical parameter S.N. Symbols represent the measured data, correlations for Series 1 and 2 are shown by trend lines (linear regression); no trend line was drawn for Series 3. The arrow indicates the evolution from "newer" to "older" snow. coefficient of water vapor in air, D 0 is evaluated at the mean temperature of the experiment (Massman, 1998), and m H 2 O is the mass of a water molecule. The simple model vapor flux j , the estimate by PIV j PIV , and the numerical solution j FE are compared in Fig. 11. Obviously, the PIV calculations suffer from large fluctuations, which we attributed to uncertainties in sample positioning combined with a strong local shape change that resulted in a mis-calculated displacement vector. Assuming that the fluctuations are of statistical nature, the mean value as well as the trend still give useful information. However, the j PIV curves (also if they were smoothed) are lower than the j FE curves. One likely source of error is that the PIV method operates on voxels, which means that any apparent displacement smaller flux j calculated from continuum theory, ignoring all microstructure and considering snow as a continuum with the diffusion constant of air. Jagged lines: flux j PIV measured by particle image velocimetry. Symbols: flux j FE calculated by numerically solving for the vapor pressure gradients in the microstructure for each time-step independently. Data are shown for all three series, each with different initial porosity and grain size. The vapor flux depends on both absolute temperature and temperature gradient, but is independent of porosity and microstructure. than half a voxel escaped unnoticed. With a typical displacement of two voxels in our case, the PIV method must be considered less reliable than the FE method. A second reason is the smoothing effect due to the correlation window and therefore a systematic underestimation of the flux. While the first issue could be addressed by increasing the resolution of the CT measurement, the smoothing due to a finite window will always play a role.
Despite the possible error in j PIV , Fig. 11 displays two striking features. First, the three methods of calculation coincide reasonably well to conclude that a diffusion enhancement of a factor of 4-5 is not likely for our samples. Second, the flux stays constant over time, despite the dramatic changes that occurred in the structure (see Fig. 5). These observations mean that, for the vapor flux through any plane in the snow, the surrounding structure does not play a role. In other words, the macroscopic vapor flux in snow can be calculated once the temperature gradient and the mean temperature of the snow are known, independently of the microstructure.
It is important to note that there is no fitting parameter involved in these calculations. The ability to derive such an important parameter as the vapor flux without the necessity of approximating the structural details in any way is not only important for the general understanding of vapor transport in snow, but also simplifies future modeling.
As already mentioned, the concept of diffusion enhancement is controversial. While Yosida (1955) and Colbeck (1993) postulated this concept, Giddings and LaChapelle (1962) argued on theoretical grounds that the influence of the snow structure is negligible. With the measurements of the vapor flux in this study, we provided experimental support that the arguments of Giddings and LaChapelle (1962) are correct. Although we covered only a small range of densities and temperature gradients with our experiments, we expect the same to hold true for other snow types as well, based on the theoretical arguments. The lattice model developed by Gubler (1985) and refined by Colbeck (1993) is based on disconnected ice spheres. This type of mass distribution concentrates the mass flux in the gap zones between the spheres. Consequently, this leads to an enhanced macroscopic mass flux in these gap zones. When the flux is averaged over the possible z-positions, the values for the effective diffusion constant go down and approach the one in air.
The fact that we did not observe a significant drop in j PIV also allows a general conclusion about the deposition coefficient α. As discussed earlier, it is expected that α < 1 over faceted surfaces, but no reliable measurements exist to determine the functional form of α versus temperature. With the development of facets over time, and therefore a decrease of α over a small fraction of the ice surface, one would expect a slowing down of the diffusion. We did not observe such an effect in our samples, which means that either the fraction of facets or the decrease of α (or both) were too small to be measured. Therefore, the systematic error in j FE introduced by assuming α = 1 is negligible.
In a natural snow cover -in particular close to the surfacea constant temperature gradient over three weeks is not very common. The fluxes and recrystallization of crystals become more complex when the temperature gradients change over time, or even change sign. The vapor flux will change proportionally to the magnitude of the gradient, and the sign of the gradient determines the direction of the vapor flux. Pinzer and Schneebeli (2009b) made first experiments with a sinusoidally changing temperature gradient. They show that the shape of the evolving snow is very different from TGM under a steady gradient, although the magnitude of the vapor flux is comparable.

Conclusions
We conducted in situ time-lapse tomographic experiments on metamorphosing snow samples under a static temperature gradient that was sufficiently strong to cause the development of depth hoar. The non-destructive nature of X-ray microtomography revealed details of the sublimation-deposition process that could not be observed before: the precise location of sublimation-deposition sites and the amount of ice that was transferred through the pore space over time.
Observing the ice structures as they grow over time helps to show how temperature gradient metamorphism works at the relevant scales. The residence time (lifetime) of ice voxels was directly derived from the data and turned out to be on the order of 2-3 days. Although hints towards a complete turnover appear already in the work by Sturm and Benson (1997), the present study provided the first experimental observation of lifetime of snow grains, and the lifetime was surprisingly low. This was consistent with the high mass turnover of up to 60 % of the total ice mass per day. In view of these highly dynamic changes, it would be confusing to write about growing ice grains, since after a few days the entire ice mass was relocated and nothing of the old core of the grain remained. The process we observed would be better described as "growth by replacement", as larger -and more vertically oriented -structures survive longer than smaller and horizontally oriented ones. Dadic et al. (2010) observed air bubble migration in bubbly ice under a temperature gradient, which is conceptually very similar to our experiments in snow. In fact, the bubbles migrate due to sublimation and deposition at the warmer and colder sides, respectively.
Using the 4-dimensional data, we addressed a long standing problem in snow science: the macroscopic vapor transport through snow. Estimating the vapor flux with both particle image velocimetry and finite element simulations, we showed that, for our experiments, the influence of the ice structure was negligible. Indeed, the 10 % difference between the three different methods indicated that within the density range between 250-310 kg m −3 , which is typical for seasonal snow that is at least a few days old, the temperature gradient and its sign alone were sufficient to calculate the vapor flux.
This lack of dependence on ice structure is in accordance with theoretical considerations, but contradicts the concept of diffusion enhancement (Colbeck, 1993). Our data provide evidence to support the argument that there is no diffusion enhancement in snow (Sokratov and Maeno, 2000).
The terms "grain-to-grain" flux and "layer-to-layer" flux have been used in the literature to describe the net fluxes that are responsible for grain growth and mass changes within a layer, respectively. They should not be confused with the underlying continuous sublimation-deposition process, which is responsible for a much higher flux.
Additional experiments over an extended density range, with the possibility to apply more complex conditions (Pinzer and Schneebeli, 2009a), will help to complete the picture. Our previous work shows the importance of alternating temperature gradients on snow metamorphism (Pinzer and Schneebeli, 2009b). We think that in situ time-lapse tomography combined with FE-simulations is an important tool to understand snow metamorphism.