Tomography-based monitoring of isothermal snow metamorphism under advective conditions

Time-lapse X-ray microtomography was used to investigate the structural dynamics of isothermal snow metamorphism exposed to an advective airflow. The effect of diffusion and advection across the snow pores on the snow microstructure were analysed in controlled laboratory experiments and possible effects on natural snowpacks discussed. The 3-D digital geometry obtained by tomographic scans was used in direct pore-level numerical simulations to determine the effective permeability. The results showed that isothermal advection with saturated air have no influence on the coarsening rate that is typical for isothermal snow metamorphism. Isothermal snow metamorphism is driven by sublimation deposition caused by the Kelvin effect and is the limiting factor independently of the transport regime in the pores.


Introduction
Snow is a bi-continuous material consisting of fully connected ice and pore space (air) (Löwe et al., 2011).Because of the proximity to the melting point, the high vapour pressure causes a continuous recrystallization of the snow microstructure known as snow metamorphism, even under moderate temperature gradients (Pinzer et al., 2012;Domine et al., 2008).The microstructural changes of snow towards equilibrium under conditions of constant temperature are referred to as isothermal snow metamorphism (Colbeck, 1997a;Kaempfer and Schneebeli, 2007).This is a coarsening process whose driving force is the reduction of the surface free energy of the complex ice-air interface.The energy reduction is caused by mass transport processes such as vapour diffusion (Neumann et al., 2009), surface diffu-sion (Kingery, 1960b), volume diffusion (Kuroiwa, 1961), and grain boundary diffusion (Colbeck, 1997a(Colbeck, , 1998(Colbeck, , 2001;;Kaempfer and Schneebeli, 2007).Viscous or plastic flow (Kingery, 1960a), and sublimation condensation with vapour transport (German, 1996;Hobbs and Mason, 1963;Legagneux and Domine, 2005;Maeno and Ebinuma, 1983) are also suggested to play an important role.The Kelvin effect is seen as the driving force for isothermal snow metamorphism (Bader et al., 1939;Colbeck, 1980).Recent studies indicate that sublimation deposition is the dominant contribution for temperatures close to the melting point, whereas surface diffusion dominates at temperatures far below the melting point (Vetter et al., 2010).Snow has a high permeability, which facilitates diffusion of gases and, under appropriate conditions, airflow (Gjessing, 1977;Colbeck, 1989;Sturm and Johnson, 1991;Waddington et al., 1996).Both diffusion and advective airflow affect heat and mass transports in the snowpack (Cunningham and Waddington, 1993;Albert, 1993;McConnell et al., 1998).In the dry snow zone of an ice sheet, Sowers et al. (1992) described a convective zone located just below the surface in which the air is rapidly flushed by convective exchange with the overlying atmosphere.A rapid decrease of the airflow velocity inside a snow layer ( ≤ 0.01 m s −1 ) for high wind speed (≈ 10 m s −1 ) above the snow surface (pore size ≈ 1 mm) are numerically estimated by Neumann (2003).In addition, Colbeck (1997b) confirmed the rapid decrease of airflow velocities inside a snowpack.Advective flow of air may have a direct effect on snow-air exchange processes related to atmospheric chemistry (Clifton et al., 2008;Grannas et al., 2007), and snow metamorphism (Albert and Gilvary, 1992;Albert et al., 2004), and can change the chemical composition of trapped atmospheric gases in ice cores (Legrand and Mayewski, 1997;Neumann and Waddington, 2004;Severinghaus et al., 2010).However, no prior studies have experimentally analyzed the effect of saturated airflow on the vapour transport and the recrystallization of the snow crystals using non-destructive technique in time-lapse experiments.Over-or undersaturated air leads to a rapid growth or shrinkage of snow structures exposed to such conditions, as exemplified in the growth of surface hoar (Stössel et al., 2010).However, saturation vapour density of the air is reached in the pore space within the first 1 cm of the snow sample, regardless of temperature or flow rate (Neumann et al., 2009;Ebner et al., 2014).The change in shape of the snow crystals during metamorphism also affects the permeability, which, in turn, will continue to affect the shape of the snow structure.Although long-term isothermal metamorphism occurs in nature only in the centre of the polar ice caps (Arnaud et al., 1998), it is important to reduce physical complexity of experiments in order to understand the basic mechanisms governing metamorphism.
The objective of this paper is to study the effect of saturated airflow on the vapour transport and the coarsening rate of snow under isothermal conditions.We designed experiments in a controlled refrigerated laboratory and used time-lapse computed tomography (micro-CT) to obtain the discrete-scale geometry of snow (Schneebeli and Sokratov, 2004;Kaempfer and Schneebeli, 2007;Pinzer and Schneebeli, 2009;Chen and Baker, 2010;Pinzer et al., 2012;Wang and Baker, 2014;Ebner et al., 2014).The extracted 3-D digital geometry of the snow was used to calculate the specific surface area and porosity.Direct pore-level simulations (DPLS) were applied to determine the effective permeability by solving the corresponding mass and momentum conservation equations (Zermatten et al., 2011(Zermatten et al., , 2014)).

Methodology
Isothermal experiments with fully saturated airflow across snow samples were performed in a micro-CT at laboratory temperatures of T lab = −8 and −15 • C. Figure 1 shows a schematic of the experimental setup (Ebner et al., 2014).It is to be noticed that the accurateness of the isothermal conditions between the top and base of the sample was less than 0.2 • C and a temperature gradient of about 6.7 K m −1 was possible.However, this was still in the uncertainty of the thermistors ±0.2 K (Ebner et al., 2014) and therefore a quasiisothermal condition was given.Two different snow types with high specific surface area were considered to evaluate the structural change in the earlier stage of isothermal metamorphism of new snow, more in detail.Partly decomposed snow (DFdc) was used for low flow rate ("sa1" and "sa2") whereas large rounded snow (RGlr) was used for higher flow rate ("sa3" and "sa4") to prevent destruction of the fragile snow structure (Fierz et al., 2009).Snow identical to natural snow was used for the snow sample preparation (wa-ter temperature: 30 • C; air temperature: −20 • C) (Schleef et al., 2014).It was sieved with a mesh size of 1.4 mm into two boxes, and sintered for 13 and 27 days at −15 and −5 • C, respectively, for increasing strength and coarsening (Kaempfer and Schneebeli, 2007).A cylinder cut out (diameter: 53 mm; height: 30 mm) from the sintered snow was filled into the sample holder (Ebner et al., 2014).The snow samples were analysed during 96 h with time-lapse micro-CT measurements taken every 8 h, producing a sequence of 13 images.Table 1 summarizes the morphological parameters of the snow.Four different runs were chosen based on the Peclet number (Pe = u D d p /D where u D is the superficial velocity in snow, d p is the pore diameter, and D = 2.036×10 −5 m 2 s −1 is the diffusion coefficient of water vapour in air) to compare the advective and diffusive transport rates inside the pore space.Experimental runs were performed at 1 atm pressure and volume flow rates of 0 (no advection), 0.36, 3.0, and 5.0 L min −1 , corresponding to Pe = 0, 0.05, 0.47, and 0.85.Higher Pe numbers were experimentally not possible, as the shear stress due to airflow could destroy the snow structure, and we restricted the flow rate to the corresponding maximum Pe ≈ 0.8 extracted from the simulations of Neumann (2003) and Colbeck (1997b).Assuming an isothermal snowpack, Pe > 1 is unlikely in nature because of the following: (1) low density snow, which always has a very low strength, will be destroyed due to the high airflow velocity; (2) Pe depends on the temperature due to changing diffusivity.Seasonal temperature fluctuations of −60 to −30 • C are typical for surface snow layer in Antarctic regions, and lead to Pe variations of up to 25 %.Theoretically, Pe ≈ 1.2 could be realistic at −60 • C for a superficial velocity of ≈ 0.06 m s −1 (experiment "sa4").However, simulations by Neumann (2003) showed a rapid decrease of the airflow velocity inside the snow layer (≤ 0.01 m s −1 ) for a high wind speed (≈ 10 m s −1 ) above the snow surface (pore size ≈ 1 mm).This leads to a maximum Pe ≈ 0.8; (3) Pe > 1 would be possible for depth hoar (Alley et al., 1990) or kinetic growth crystals, like sublimation crystals (Gallet et al., 2014, Adams andWalter, 2014) close to the surface but they were only formed under light winds conditions.According to the reported Beaufort number (in Alley et al., 1990), this will be a maximum wind speed of ≈ 2-3 m s −1 (see also Gallet et al., 2014) above the surface.In addition, they were developed in the slopes of older dunes, leading to an additional decrease of the actual wind speed (≈ 1 m s −1 ) above the layers.Based on the simulations of Neumann (2003) an airflow velocity inside the snow layer of ≤ 0.002 m s −1 would be realistic.To reach a Peclet number > 1 under this condition, the mean pore size must be at least 10 mm, which would be a very extreme case for depth hoar or the kinetic growth crystals formed close to the surface.Additionally, Adams and Walters (2014) showed that the top layer of such kinetic growth crystals consists of long slender needle Table 1.Morphological and flow characteristics of the experiments: volume flow ( V ), corresponding Peclet number (Pe), Reynolds number (Re), initial superficial velocity in snow (u D,0 ), initial snow density (ρ 0 ), initial porosity (ε 0 ), specific surface area (SSA 0 ), initial pore diameter (d p ), temperature in the cold laboratory (T lab ), and the sintering time of the snow.crystals connected in a cross-hatch pattern which has a low strength and will be destroyed by such a strong flow.The acceleration voltage in the X-ray tube was 70 kV, with an intensity of 114 µA, and a nominal resolution of 18 µm.The samples were scanned with 2000 projections per 360 degree, with an integration time of 200 ms per projection, taking 1.5 h per scan.The innermost 36.9 mm of the total 53 mm diameter were scanned, and subsamples with a dimension of 7.2 mm × 7.2 mm × 7.2 mm were extracted for further processing.Absolute z position varied up to a maximum of 50 voxels between subsequent scans due to the weight of the sample holder.To correct for the z position, a linear encoder was built into the micro-CT.A 3 × 3 × 3 median filter and Gaussian filter (σ = 1.4,support = 3) was applied to the reconstructed images.Otsu's method (Otsu, 1979) was used to automatically perform clustering-based image thresholding to segment the grey-level images into the ice and air phases.Morphological properties in the two-phase system were determined based on the geometry obtained by the micro-CT.The segmented data were used to calculate a triangulated ice matrix surface and tetrahedrons inscribed into the ice structure.Morphological parameters such as porosity (ε) and specific surface area (SSA) were then calculated.The opening size distribution with spherical structuring elements on the micro-CT scans was used to estimate the mean pore size (d p ) (Haussener et al., 2012).The effective permeability was calculated using the finite volume technique CFD (computational fluid dynamics simulation software from ANSYS; ANSYS, 2010) by solving the continuity and Navier-Stokes equations (Zermatten et al., 2011(Zermatten et al., , 2014) ) for laminar flow where p is the pressure, µ is the dynamic viscosity of the fluid and u D its superficial velocity, ρ is the fluid density, K is the permeability, F is the Dupuit-Forchheimer coefficient, and γ is a dimensionless factor.The first term is the result of viscous effects, predominant at low velocities, whereas the second and third terms describe the inertial effects, which become important at higher fluid velocities.
As the viscous effect was still the dominant case (Re ≈ 1) in the experiment, only permeability K was considered for further discussions.A grid convergence study based on the pressure drop (Zermatten et al., 2014) was carried out to find the optimal representative elementary volume (REV) (6.0 mm×6.0 mm×3.0 mm).An in-house tetrahedron-based mesh generator (Friess et al., 2013) was used to create the computational grid on the segmented data.The computational domain consisted of a square duct containing a sample of snow.The boundary conditions consisted of uniform inlet velocity, temperature and outlet pressure, constant wall temperature at the solid-fluid interface, and symmetry of the sample at the lateral duct walls.The square duct was 5 times the length of the sample to ensure a fully developed velocity profile at the entrance of the snow sample (Fig. 2).The largest mesh element length was 0.153 mm, and the smallest possible mesh element measured 9.56 µm, with average 60 million volume elements for each segmented snow sample.

Results and discussion
The discussions of the observed results are only based on the investigated volume.Influences of the flow on the base, top and lateral boundaries of the overall sample were not considered due to lack of structural observations.A representative temporal temperature profile of the snow sample for both laboratory temperatures of T lab = −8 and −15 • C is shown in Fig. 3. Variations in temperature up to 1.7 and 1.4 • C were due to heat dissipated by the Xray tube and temperature fluctuations inside the cold laboratory (Ebner et al., 2014).A longer sintering duration at −5 • C of the snow for experiment "sa3" and "sa4" was used to increase the mean thickness of the ice matrix.This avoided the destruction of the snow structure due to shear stresses caused by the airflow.The structural analysis of the snow samples was conducted on the complete tomography domain (7.2 mm × 7.2 mm × 7.2 mm).A smaller subset of 110×42×110 voxels (2 mm×0.75 mm×2 mm) was selected to visualize the 3-D evolution (Fig. 4).It showed no significant change in the grain shape, even for different airflow velocities, and only a slight rounding and coarsening was seen for experiments "sa1" and "sa2".A strong translation effect due to settling of sub-layering snow was visible for "sa1" and "sa2".The initial ice matrix did not change with time; only coarsening processes on the ice grain surface were observed (Fig. 5).Sublimation of 4.5 and 4.9 % of the ice matrix and deposition of 4.1 and 5.9 % on the ice matrix were observed for "sa3" and "sa4" (Fig. 6).The data were extracted by superposition of vertical cross-sections at 0 and 96 h with an uncertainty of 6 %.The mass sublimated preferentially at locations of the ice grain with low radii due to Kelvin-effect and was relocated on the grain leading to a smoothing of the ice grain.Our observed results were supported by the vapourpressure map simulated by Brzoska et al. (2008) and the applied airflow velocity did not affect the relocation process.
The well-sintered snow showed very little settling under its own weight (Kaempfer and Schneebeli, 2007) and, consequently, no significant change in porosity was observed.This supports the hypothesis that further densification is limited by coarsening kinetics (Kaempfer and Schneebeli, 2007; A typical temperature profile for experiment "sa1", "sa2" and "sa3", "sa4".The temperature rise was caused by the X-ray tube and fluctuations inside the cold laboratory (Ebner et al., 2014).The accurateness of the isothermal conditions between the top and base of the sample throughout the experiment is less than 0.2 • C which is still in the uncertainty of the thermistors ±0.2 K (Ebner et al., 2014).Schleef and Löwe, 2013).A spatially constant porosity distribution at t = 0 and t = 4 days is seen in regime.The average deviation between t = 0 and t = 4 days was 0.5, 1.8, 0.5, and 0.5 % for "sa1", "sa2", "sa3" and "sa4".
Our segmented 3-D data accurately reproduced the original snow sample and the temporal porosity distribution confirmed that no settling and densification occurred in the investigated volume (Fig. 8).The gravimetric porosity ε grav at the beginning and at the end of each experiment was measured by weighing.The measured density values were converted to porosity (ε grav = 1 − ρ s /ρ ice ), and compared to the value of porosity computed by DPLS on the micro-CT geometry.The computed values differed from the measured ones by 1.4 and 0.1 % at the beginning and 4.1 and 2.3 % at the end for experiments "sa3" and "sa4".The qualitative progression of the spatial SSA of the scanned snow height for four discs of 7.2 mm × 7.2 mm × 1.8 mm (Fig. 9) did not change significantly with height.This suggested that the snow properties were homogeneous throughout the sample and duration of the experiments.The slight decrease of the spatial SSA for experiment "sa4" is explained by the distribution not initially being completely homogeneous.The coarsening process led to a decrease of the SSA over time (Fig. 10), which was higher for group "sa1" and "sa2" compared to "sa3" and "sa4".The difference was caused by the 34 % lower initial SSA of group "sa3" and "sa4".Applying the theories developed by Legagneux et al. (2004) and Legagneux and Domine (2005), the evolution of SSA of the ice matrix could be modelled well.The model proposed is  given by Legagneux and Domine (2005 where SSA 0 is the initial SSA at time t = 0, n is the growth exponent, and τ a parameter related to grain growth and a form factor. Table 2 shows the fitted parameters and the cor-  (Ratke and Voorhees, 2002).Ostwald ripening describes the coarsening of solid particles with a given size distribution, considering disconnected grains that do not undergo settling.The driving force in the model is the reduction of the SSA and the model hypothesis is based on the concept that mass transfer occurs by sublimation due to curvature effects, transport through the gas phase and deposition.Theoretically, the growth exponent n is approximately 2 when surface processes are rate limiting and 3 when diffusion is rate limiting.Experiment "sa1" and "sa2" had a higher value of n, indicating a strong coarsening process due to sintering and that surface processes were rate limiting (Legagneux et al., 2004;Legagneux and Domine, 2005).Experiment "sa1" and "sa2", and "sa3" and "sa4" had similar fitting parameters and a low value of n, suggesting that surface effects were rate limiting (Legagneux et al., 2004;Legagneux and Domine, 2005).The lower value of n for experiment "sa3" and "s4" was due to the longer sintering time of 27 days at −5 • C before the experiments were started leading to a very little change in the microstructure of the snow.When the sintering times of 13 and 27 days were included in the model, the fitting parameters indicated a consistent growth exponent n for each experiment (Table 3) and a good agreement with the theory.They expressed strong coarsening and surface processes for each experiment.Notice, Eq. ( 2) extremely depends on the initial state, which is well illustrated by the large difference obtained for n values of "sa3" and "sa4" between Tables 2 and 3. Concluding, the calculated values indicated that surface processes caused the limiting rate rather than the diffusion step and no significant influence of advective transport could be observed.
The effect of decreasing SSA on the permeability was not elucidated in our experiments.A SSA decrease of at least 5 % in the experiments could not be reproduced in the permeability.However, the computational uncertainty up to 16 % (Zermatten et al., 2014) in the permeability is still in the range  Zermatten et al., 2014).
to cover the correlation between SSA and permeability.The effect of increasing airflow velocity had no influence on the flow characteristics (Fig. 11).The temporal evolution of permeability for experiment "sa2" showed a decrease of 8 % for the first 40 h and remained constant afterwards.Experiments "sa1", "sa3" and "sa4" showed no significant change in the permeability, which is consistent with the negligible change in density.The average fluctuations of the permeability K between each time step and the slight decrease at the beginning in "sa2" showed small differences that were below the precision of the numerical method with an uncertainty up to 16 % (Zermatten et al., 2014).Only the first time step of "sa3" showed a particularly high difference of 17.3 %, but neither the porosity nor SSA showed significant differences reflecting this value.This difference could therefore be due to an error during the measurement or during the meshing procedure.
Four isothermal metamorphism experiments of snow under saturated advective airflow were performed, each with duration of 4 days.The effects of the main transport processes, diffusion and advection, were analysed inside the pore space.The airflow velocities were chosen based on the Peclet number.Pe > 0.85 for natural surface conditions were not possible due to the destruction of the snow structure and is not frequent in natural snowpacks due to the low airflow velocities in snow (Neumann, 2003;Colbeck, 1997b).Every 8 h the snow microstructure was observed by X-ray micro-tomography.The micro-CT scans were segmented, and porosity and specific surface area were calculated.Effective permeability was calculated in direct pore-level simulations (DPLS) to analyse the flow characteristic.
The experimental observations supported the hypothesis that further densification was limited by coarsening kinetics and further confirmed a constant porosity evolution (Kaempfer and Schneebeli, 2007).Curvature caused sublimation of small ice grains and ice structures with small curvature radii leading to a slight decrease in SSA.Compared to rates typical for isothermal snow metamorphism, no enhancement of mass transfer inside the pores of isothermal advection with saturated air was observed.Sublimation deposition caused by the Kelvin effect was the limiting factor independently of the transport regime in the pores.

Figure 1 .
Figure 1.Schematic of the experimental setup and the sample holder.A thermocouple (TC) and a humidifier sensor (HS) inside the humidifier measured the airflow conditions.Two thermistors (NTC) close to the snow surface measured the inlet and outlet temperature of the airflow (Ebner et al., 2014).

Figure 2 .
Figure2.Schematic of the computational domain with an enlarged subsample of snow.In the snow sample, the dark gray part represents the ice, whereas the mesh is built in the pore space.

Figure 4 .
Figure 4. Evolution of the 3-D structure of the ice matrix during isothermal metamorphism under advective conditions.Experimental conditions (from left to right) at different measurement times from beginning to the end (top to bottom) of the experiment.The shown cubes are 110 × 42 × 110 voxels (2 mm × 0.75 mm × 2 mm) large.

Figure 5 .
Figure 5. Residence time of ice particles within in a slice (5.7 mm× 5.7 mm) parallel to the flow direction for (a) "sa3" and (b)) "sa4" by overlapping time-lapse tomography pictures.The period of 8 h was sufficiently short to calculate the residence time of each ice voxel with an uncertainty of 6 %.

Figure 6 .
Figure 6.Superposition of vertical cross-section parallel to the flow direction at time 0 and 96 h for (a) "sa3" and (b) "sa4".Sublimation and deposition of water vapour on the ice grain were visible with an uncertainty of 6 %.

Figure 7 .Figure 8 .
Figure 7. Spatial porosity profile of the scanned area at the beginning and at the end of each experiment.The spatial variability within the reconstructed volume was measured in four discs of 7.2 mm × 7.2 mm × 1.8 mm.

Figure 9 .
Figure 9. Spatial SSA profile of the scanned area at the beginning and at the end of each experiment.The spatial variability within the reconstructed volume was measured in four discs of 7.2 mm × 7.2 mm × 1.8 mm.

Figure 10 .
Figure 10.Temporal evolution of the specific surface area, SSA, of the ice matrix obtained by triangulated structure surface method.The computed fit is of the form SSA(t) = SSA 0 τ τ +t 1/n .

Table 2 .
Values of the fitted growth rate τ and growth exponent n for the evolution of the SSA and the corresponding normalized root mean square error (NRMSE).

Table 3 .
Values of the fitted growth rate τ and growth exponent n for the evolution of the SSA including the sintering time of 13 and 27 days, and the corresponding normalized root mean square error (NRMSE).