Introduction
In the terrestrial environment, ice exists near its phase-change temperature. Consequently, snow on the ground is a thermodynamically active material with a granular microstructure that is continuously changing. The morphologic structure of the snowpack influences virtually all of its thermomechanical and optical properties including strength, density, thermal conductivity, heat capacity, and reflectivity or albedo. These properties dictate many critical snow responses and are foundational to snow science. The critical dependence of snow properties on microstructure makes metamorphism an extremely important area of study. Typically, dry snow will metamorphose either toward a faceted granular structure when exposed to a temperature gradient or to rounded sintered ice grains under macroscopically isothermal conditions. Both processes may occur simultaneously at different levels within the same snow-pack. Which of these processes dominates is largely dependent on the magnitude of the local temperature gradient. Larger temperature gradients, nominally taken to be >10 K m−1 (Reference ColbeckColbeck, 1982), lead to the development of the structurally fragile snow with a faceted crystal morphology, termed the kinetic growth form (Reference ColbeckColbeck, 1982). The process is particularly prevalent near the relatively warmer snow base, in which case the resulting structure is termed depth hoar (Reference SeligmanSeligman, 1936). Lower temperature gradients yield a structure composed of rounded sintered grains, commonly referred to as the equilibrium growth form (Reference ColbeckColbeck, 1987).
Several theories have been published to explain metamorphism mechanisms for rounded grains and faceted grains individually, but limited research spans both regimes (Reference Arons and ColbeckArons and Colbeck, 1995). Reference Lehning, Bartelt, Brown, Fierz and SatyawaliLehning and others (2002) present approaches using mixture theory to describe ‘equi-temperature’ metamorphism and vapor transport within a lattice network of cubic grains to model kinetic growth. Reference Christon, Burns and SomerfeldChriston and others (1994) explored kinetic growth metamorphism by solving heat- and mass-conservation as well as phase-change differential equations using three-dimensional finite-element techniques. Their approach did not include the effects of ice surface geometry on excess vapor pressure, which are important for equilibrium growth but not a significant limiting factor during kinetic growth. As such, their model focused on temperature gradient metamorphism only. Reference Flin and BrzoskaFlin and Brzoska (2008) present a simple vapor diffusion model for kinetic growth that includes crystal modeling to predict rounded and faceted morphologies based on temperatures and surface geometry. They use a Kossel crystal model for sublimation and condensation, but neglect latent heat in the ice network. Reference Flanner and ZenderFlanner and Zender (2006) address curvature growth and temperature gradient growth with simplified spherical ice geometries to study the evolution of snow albedo. They present a coupled system with smooth transitions between the two regimes, but are only concerned with rounded ice geometries. Reference Kaempfer, Schneebeli and PinzerKaempfer and others (2008) and Reference Kaempfer and PlappKaempfer and Plapp (2009) utilized a microcomputed tomography imager to non-destructively image snow undergoing temperature gradient metamorphism. They used microstructure images to develop a model for a numerical tool that examined the heat and mass transfer with phase change through a diffuse boundary. Their results suggest heat- and mass-flow ‘inhomogeneities’ on the size scale of grains. Their approach is unique with the coupling of experimental methods and numerical simulation. Reference Miller, Adams and BrownMiller and others (2003) developed a theory for snow metamorphism independent of temperature conditions. They presented simplified ice-grain and pore-space geometries where vapor mass in the pore and energy in the ice were conserved. Latent-heat transfer into or out of the ice network was included. They coupled the conservation equations with heat and mass sources from phase change at the ice/pore interface that included ice surface energy. The resulting coupled nonlinear equations were solved by finite-difference methods for discretized ice and pore networks. Reference Miller, Adams and BrownMiller and others (2003) showed good correlation with rounded morphologies and identified transition temperature gradients for kinetic growth. They examined important parameters defining the transition to kinetic growth, but did not present an approach where individual grains were allowed to grow kinetically after transition.
The goal of this research is to utilize Reference Miller, Adams and BrownMiller and others’ (2003) approach as the basis for heat and mass transfer and expand it to include a kinetic growth capability that considers hopper-type crystal habit. In keeping with a unified approach, the method should be independent of temperature boundary conditions and should always be applicable as conditions dictate. The resulting model could simulate a broad spectrum of dry snow microstructure metamorphism.
Model Development
Existing model background
The current project uses the heat, mass and phase-change interaction developed by Reference Miller, Adams and BrownMiller and others (2003). Before the development of a new kinetic growth model is introduced, their approach is briefly summarized. The microstructural geometry and governing equations are described by Reference Miller, Adams and BrownMiller and others (2003). Details with derivations are found in Reference MillerMiller (2002). They modeled a network of rounded ice grains connected to each other by concave necks with accompanying pore spaces. Their network was the basis of this study and was used to extend the approach to kinetic growth. Reference Miller, Adams and BrownMiller and others (2003) present a quasi-two-dimensional approach, as heat and vapor are free to move in one dimension with phase-change interaction between vapor and solid phases in a second dimension (e.g. Fig. 1). Metamorphism was simulated through solutions of the heat and mass conservation equations, which were coupled through phase change. In the pore spaces, the steady-state mass conservation equation with mass sources/sinks and Fick’s law for vapor diffusion (including the differential form of the Clausius–Clapeyron equation) are
J is the mass flux from diffusion through the pore (kg m−2 s−1), J ec is the mass flux to or from the pore due to phase change (kg m−2 s−1), A ec is the ice surface area undergoing phase change (m2), V pore is the volume of the pore space (m3), D is the water-vapor diffusion coefficient in air (2.02 × 10−5 m2 s−1),R is the water-vapor gas constant (462 J kg K−1), T is the pore temperature (K), L is the water latent heat of sublimation (2.838 × 106 J kg K−1), MS is the mass supply or sink to the pore space from phase change (kg m −3 s−1)and P is the pore pressure (Pa) (Reference Miller, Adams and BrownMiller and others, 2003). By assuming the vapor is an ideal gas, treating the pore space as saturated and assuming a temperature gradient in only one dimension, the relations in Equation (1) can be combined and expanded giving a non-linear ordinary differential equation in terms of pore temperature (Reference ColbeckColbeck 1982, Reference Colbeck1983) but including mass sources or sinks from phase change (Reference MillerMiller, 2002). The mass vapor flux through the pore space and mass flux due to phase change are represented in Figure 1.
Fig. 1. Relationship between crystal habit (β), crystal axes relative to temperature gradient (γ) and crystal orientation relative to temperature gradient (α). Mass flux through the pore space (J) and mass flux due to phase change ( J ec) are shown.
To solve the mass conservation differential equation, Reference Miller, Adams and BrownMiller and others (2003) chose to examine constant temperature conditions. While temperatures within the microstructural model were allowed to change during metamorphism, the boundary nodes were maintained at constant temperature. This implies no limit on availability of vapor entering or leaving the model boundary, except as prescribed by vapor diffusion. Vapor crossing the model boundaries comes from or goes to the surrounding snowpack (e.g. mass supply by sublimation of ice grains in what is often referred to as the hand-to-hand transfer of mass (Reference YosidaYosida and others, 1955)). Temperature gradients in the pore are used to calculate mass flux crossing boundaries. Given the small scales of ice and pore structures considered, and their intent to individually model small areas at various locations within a snowpack, this assumption was deemed reasonable.
Next, Reference Miller, Adams and BrownMiller and others (2003) considered conservation of energy through the ice network. The steady one-dimensional heat conduction equation through the ice network with variable area is given as
k is the ice thermal conductivity (2.2 W m−1 K−1), A is the cross-sectional area for heat conduction through the ice grain and neck network (m2), θ is the ice temperature (K), y is the dimension parallel to the temperature gradient (m), and HS is a heat source term that is coupled to the pore through the phase change at the boundary (W m−1) (Reference Miller, Adams and BrownMiller and others, 2003). Equation (2) is very similar to heat conduction of a variable geometry cooling fin except convective heat exchange is replaced by heat from phase change (Reference Incropera and DeWittIncropera and DeWitt, 1985). Since the ice geometry changes in time and the phase-change boundary is moving, Reference LockLock (1990) defines this as a ‘quasi-steady’ problem even though heat capacitance terms are neglected. Notice that the heat (HS) and mass (MS) source terms are coupled through the phase change mass flux ( J ec) at the ice/vapor interface.
Reference Miller, Adams and BrownMiller and others (2003) then present the energy balance at the phase change surface accounting for ice surface geometry, ice surface temperature, and pore temperature. The result is the third coupled differential equation,
where x is the direction of latent-heat flow during phase change (m).
Reference Miller, Adams and BrownMiller and others (2003) finally offer a finite-difference iterative solution to Equations (1–3) that ensures mass conservation in the pore, energy conservation in the ice (including latent heat) and accounts for phase change. The geometry is then updated, based on integrating mass flux from phase change, as the simulation continues. Reference Miller, Adams and BrownMiller and others (2003) used this model to examine density, grain size, bond size and temperature for macroscopically isothermal conditions with excellent correlation to established trends and experimental data. They replicated observed phenomena such as large grains growing at the expense of small grains and intergranular bond development. In the presence of a temperature gradient, a transition to kinetic growth was defined when the vapor pressure of the pore becomes sufficient to overcome convex surface curvature effects of the ice. This model is used here as the basic tool for the movement of vapor in the pore space and heat in the ice network with the phase-change interaction between them. With this approach established, a kinetic growth extension to the model is now offered.
Facet geometry and crystal habit
The hexagonal crystal structure of ice Ih gives it unique crystal habit configurations when grown from vapor. Depending on temperature, either the basal or prism plane will generally exhibit dominant growth. Their relative growth rates determine the primary geometric habit of faceted crystals (Reference LockLock, 1990; Reference LibbrechtLibbrecht, 2005). If c-axis growth dominates, prism-like structures appear. Conversely, plate-like structures form when a-axis growth is dominant. Early experiments by Reference Nakaya, Sato and SekidoNakaya and others (1938) showed that crystal habit varied with temperature and supersaturation. While specific environmental and crystal habit relationships were not developed, their experiments indicated that basic crystal habit was primarily determined by temperature. Several authors have determined the effects of supersaturation and temperature on the growth of ice from the vapor, but Reference KobayashiKobayashi’s (1957, Reference Kobayashi and Oura1967) experiments rigorously confirmed primary crystal habit dependence on temperature, with secondary features such as dendrite growth or hopper crystals (e.g. Fig. 2c) determined by vapor supersaturation.
Fig. 2. (a) Facet crystal geometry showing basal plane. (b) Two crystal facet faces with relative a- and c-axis growth determining primary crystal habit. The aspect ratio of the crystal is determined by a- and c-axes growth velocities. β is the angle from the dominant growth direction to the crystal face and is a function of the relative growth rates of each crystallographic axis. All faces for a particular crystal have the same angle due to symmetry assumptions. The size of the crystal is taken as the length of the crystal edge l/sin(β). (c) Examples of modeled single-facet faces superimposed on hopper-type crystals. The crystal on the left is c-axis dominant growth and the crystal on the right is a-axis dominant growth. The model will accommodate one to six faces. (Hopper crystal representations taken from Reference LockLock, 1990.)
While the very complex features of many crystal habits with complicated secondary features are difficult to model (e.g. Reference LibbrechtLibbrecht, 2005), primary crystal habit in this model is incorporated through simple geometric shapes. Here the faceted crystals nucleate and grow on Reference Miller, Adams and BrownMiller and others’ (2003) round ‘parent’ grains. The rounded form provides the base to grow the kinetic crystal morphologies. If the relative crystallographic growth rates (a vs c axes) are known, the primary crystal habit is established for kinetic growth. Reference Lamb and HobbsLamb and Hobbs (1971) presented relative crystallographic growth rates as a function of temperature. Using their data at specific temperatures, crystal habit is determined here by selecting the relative growth velocities of the a and c axes. Only the crystal habit is utilized from Reference Lamb and HobbsLamb and Hobbs (1971); actual growth velocities are addressed later. Figure 2 shows the proposed simple facet geometry to include primary crystal habit. The simple geometric structures considered here approximate hopper-type structures that result from high supersaturation growth (Reference LockLock, 1990). This geometry is based on the hopper crystal habit since the cup (or partial cup) morphology is often observed in temperature-gradient-driven snow metamorphism (Fig. 2c). The full symmetry of the well-developed atmospherically grown hopper crystal will not be expected to develop when confined to the snowpack constraints. The approach only considers combinations of prismatic shapes with the aspect ratio of the triangular surface determined by relative crystallographic growth velocities. Each triangle is referred to as a facet face whose shape is determined by temperature. The approach also allows for one to six faceted crystal faces to develop. In Figure 2a, the basal plane showing the characteristic prism length, l, is illustrated. The c axis is perpendicular to the basal plane. Figure 2b shows two sides of a potential crystal with combined a- and c-axis growth. In Figure 2c, a single facet face from the model is superimposed on hopper crystals for plate and needle-type morphologies. At any particular temperature, one axis generally has a greater growth velocity than the other (unless they are exactly equal, which is possible at three temperatures (Reference Lamb and HobbsLamb and Hobbs, 1971)). In order to quantify and model the growth rate of each axis, an angle from the dominant growth direction (a or c depending upon temperature) to the crystal facet face is defined by the angle β. β is determined by the relative growth velocity of the two crystallographic directions. If single axis growth is prevalent, β will be small, but the approach can accommodate all combinations of relative growth rates. If a- and c-axes growth rates are equal, β will be 45°. Since β is measured from the dominant growth axis, it is never greater than 45°.
Once crystal habit has been established, the actual growth velocity of the dominant axis can be calculated. Spiral defect propagation theory offers one explanation of the kinetic growth of ice crystals. Reference Burton, Cabrera and FrankBurton and others (1951) developed an extensive model for step propagation enabling the calculation of crystallographic face growth velocity, which is well supported in the literature. Reference HobbsHobbs (1974) and Reference LockLock (1990) interpret Reference Burton, Cabrera and FrankBurton and others (1951) with the following relationship for the growth velocity of the a or c axis as
where Vl is the facial velocity (m s−1),l is either the c or a axis (depending on which is dominant), α is the molecular absorption coefficient (0.01; Reference YokoyamaYokoyama, 1993), m is ice molecular mass (2.989 × 10−26 kg), k is Boltzmann’s constant (1.38 × 10−23 J K−1), ρ ice is ice density (917 kg m−3) and ΔP is the excess vapor pressure (Pa). ΔP 1 is a critical excess vapor pressure, given by
where η is the ledge energy per unit step length (2 × 10−11 J m−1 from Reference YokoyamaYokoyama, 1993), a 0 is the lattice parameter (4.519 × 10−10 m for a-axis growth and 7.357 × 10−10 m for c-axis growth; Reference Petrenko and WhitworthPetrenko and Whitworth, 1999), X s is the mean absorption migration distance (400a 0 m; Reference YokoyamaYokoyama, 1993) and P 0 is the equilibrium vapor pressure (Pa). For ΔP 1 ≫ ΔP, VI ∝ΔP 2, and for ΔP1 ≪ ΔP, VI ∝ΔP. Such behavior has been noted in laboratory studies. Reference LibbrechtLibbrecht (2005) questions the adequacy of critical parameters such as the molecular coefficient which is likely a function of supersaturation and temperature. For this study, the parameters α, η, a 0 and X s are all assumed constant. If, in the future, more accurate representations of these values become available, the model can be easily updated. Kinetic growth is always possible if there is excess vapor pressure relative to the ice surface, but if the vapor is in equilibrium or undersaturated relative to the ice surface, kinetic growth is not possible. For understaturated conditions, sublimation is still allowed under the original Reference Miller, Adams and BrownMiller and others (2003) model. Equations (4) and (5) describe crystal growth when enough excess vapor pressure exists in the surroundings to produce kinetic growth. With the dominant axis growth velocity described in Equation (4), the overall crystal growth velocity with both axes contributing to growth is
This crystal growth velocity represents the growth of a crystal edge which results from the simultaneous growth of the a and c axes.
Heat and mass transfer during kinetic growth
With the crystal habit and the crystal growth velocity established, the heat and mass source terms used in the conservation Equations (1) and (2) need to be developed for kinetic growth. To establish source terms, the mass flux from the vapor to the ice during kinetic growth is required. Even though detailed geometry of faceted crystals has been studied little, some geometric parameters have been measured. Reference Sturm and BensonSturm and Benson (1997) conducted extensive experimental research in Alaska on depth hoar. They collected numerous depth-hoar samples for mass and size characterization. Using mass fraction techniques in various-sized sieves, Reference Sturm and BensonSturm and Benson (1997) established a polynomial fit of sieve size and grain mass as
where mj
is the grain mass average for the sieve,ψis a fit constant (126.6 kg m−3) and 
is the average size of the j and j + 1th sieves (m). This mass to grain-size approximation is useful in relating mass addition to crystal growth. Reference Sturm and BensonSturm and Benson (1997) empirically relate depth-hoar grain mass with a grain dimension that is similar to grain size used here. For the current problem, the length of a single crystal edge is used for the average sieve size (Fig. 2b), 
. While a narrow yet long crystal may fit through a sieve smaller than this dimension, the relatively planar characteristics of these crystals makes this dimensional approximation reasonable. Equation (7) becomes
With the crystal mass related to a crystal dimension by Equation (8), the thickness of the crystal can be established. Assuming a constant thickness for the entire crystal and using the triangular faces shown in Figure 2, a crystal thickness can simply be calculated. Kinetic growth is assumed to take place at the crystal’s exposed edges or ‘tips’ through rectangular surface areas (Fig. 2c). With the growth velocity and crystal geometry defined, the mass flow from phase change and the resulting latent-heat contributions may be calculated. The resulting mass source (MS) to be used in Equation (1) is
where V pore is the volume of the pore element containing the crystal tip and supplying vapor for growth (Reference MillerMiller, 2002). The negative sign for the facet mass source term represents condensation from the vapor to the ice phase. The growth velocity of the crystal is geometrically related to the dominant growth direction velocity, VI , through the cosine term appearing in Equation (9) but originating in Equation (6). The heat source term for facet growth is derived in a similar manner. Using the mass flow from phase change, the facet heat source in the ice becomes (Reference MillerMiller, 2002)
Equations (9) and (10) can now be used in Equations (1) and (2) respectively.
Crystal Orientation in the Presence of a Temperature Gradient
With crystal habit defined as a function of relative crystallographic growth rates (β) and crystal geometry developed, the orientation of a crystal within the snowpack can be considered. Crystal habit has been defined relative to crystallographic coordinates; now orientations of the crystals relative to a temperature gradient are considered. It is assumed that faceted crystals nucleate on existing or parent grains and grow with the same crystallographic orientation as the parent grain. Reference Adams and MillerAdams and Miller (2003) demonstrated that ice crystals grown from the vapor phase on existing ice substrates adopt the same crystallographic orientation as the parent. They also proposed a grain growth theory for depth hoar where dominant configurations grow more rapidly than less preferentially aligned neighbors. While the theory is supported by observation, they state further quantitative work is required to confirm the conjecture. Atmospherically formed snow crystals will be deposited onto the ground in a generally random crystallographic orientation. When subjected to a temperature gradient, the ‘new’ crystals will pick up the orientation of the parent (primary/nucleate) snow grains onto which they accrete. A particular crystal may have a growth advantage due to parent crystal orientation, since it may be better positioned to quickly grow deeper into warmer vapor (if it is not restricted by the proximity of other grains). If a particular grain’s crystallographic orientation relative to the temperature gradient requires the tip to grow toward the warmer areas, the growth rate will be enhanced as the crystal develops. For these orientations, the relatively high thermal conductivity of ice (as compared to the pore spaces) results in a cooler tip reaching into a progressively warmer (as the crystal grows) saturated pore space (i.e. supersaturated with respect to the ice surface). The result is increased excess vapor pressure at the ice tip, enhancing the growth rate. New crystals that are favorably oriented could crowd out other less optimally oriented crystals (although this is not currently implemented). Since crystals grow with the same crystallographic orientation as the parent grain, and observed depth-hoar crystal morphologies are similar amongst grains in the same snowpack location (Reference AkitayaAkitaya, 1974), this dominant grain theory seems plausible. An analogous process takes place in lake ice grown under high growth-rate conditions (Reference Gow and LangstonGow and Langston, 1977; Reference LockLock, 1990; Reference Petrenko and WhitworthPetrenko and Whitworth, 1999). Initial basal plane growth rates dominate for liquid water at approximately 273 K, allowing vertical basal planes to grow deeper faster than other orientations. As the crystals reach deeper, they grow faster and eventually crowd out c-axis vertical crystals.
To examine the dominant grain growth scenario, the variation of crystal orientation relative to the temperature gradient must be established. In Figures 1 and 2, the crystal habit was defined by the angle (β) measured from the dominant crystallographic growth axis to the crystal edge and was determined using the relative axial growth velocities as a function of temperature (Reference Lamb and HobbsLamb and Hobbs, 1971). As shown in Figure 1, the angle from the temperature gradient to the dominant growth axis is defined by the angle y. This feature is added to simulate various crystallographic orientations relative to the temperature gradient so distributions of crystal orientations in snow can be represented. Figure 1 shows the relationships between the crystal, crystallographic axes and temperature gradient.
Updating microstructural geometry
With provisions for heat and mass flow through the ice/pore network, a definition of a simple microstructural geometry with habit and orientation relative to the temperature gradient established, the numerical simulation may commence. The heat-conduction, vapor mass-continuity and phase-change equations are solved at every time-step using Reference Miller, Adams and BrownMiller and others’ (2003) model. During kinetic growth, the mass flux due to phase change, J ec, is calculated at each time-step from terms in Equations (4–8). The latent-heat contribution from phase change is added to the ice network directly at the parent grains (HS); conduction from the grain tip to root is not considered. The facet grain temperature is assumed to be constant due to high (relative to pore) thermal conductivity of the ice. The inclusion of heat conduction from the crystal tip to root at the base grain is being considered in future updates. Over a small time-step, the mass flux is integrated, giving the total addition of mass to the crystal. The crystal geometry is then updated and the simulation advances to another time-step.
Results and Analysis
A simple set of microstructural conditions was used first to examine kinetic growth trends. Then specific comparisons are made with previously published experiments of kinetic growth. In examining kinetic growth trends, following the work of Reference Miller, Adams and BrownMiller and others (2003), an initially rounded set of parent grains with 0.5 mm radius and snow density of 150 kg m−3 was used as the base network to simulate kinetic growth. The model’s sensitivity to these parameters is discussed by Reference Miller, Adams and BrownMiller and others (2003) as well as their roles in modeling metamorphism. The current effort is focused on a kinetic growth approach, so these microstructural parameters are not varied here but do impact the potential rates of heat and mass diffusion. The temperature gradient was allowed to increase until the kinetic growth on the parent grains developed. In Reference Miller, Adams and BrownMiller and others (2003), the transition to kinetic growth occurred when the excess vapor pressure over the ice surface was sufficient to require rounded parent grains to grow kinetically.
Facet growth vs temperature
The change in facet growth rate as a function of temperature is not a simple result. In addition to the temperature dependence of metamorphic processes due to molecular mobility, the dominant crystal habit varies with temperature. The growth rates of the a and c axes are not the same under identical conditions. The critical excess vapor pressure in Equation (5) is greater for c-axis growth due to a larger lattice parameter, a 0 (Reference Petrenko and WhitworthPetrenko and Whitworth, 1999). As a result, a-axis dominant growth rates are greater than c-axis growth rates under identical thermal conditions. Additionally, to isolate the temperature dependence, only one crystal orientation relative to the temperature gradient is presented. Varying crystal orientation relative to the temperature gradient is examined in the next subsection. The angle from the temperature gradient to the crystal was fixed at 20° (the combination of crystal habit and crystal orientation is shown in Figure 1). These microstructural conditions were simulated with a 45 K m−1 temperature gradient over a range of temperatures. The results in Figure 3 show the kinetic growth rate under identical temperature gradients but at different temperatures. As expected, the overall growth rate does decrease with decreasing temperature. In addition, the transitions between a- and c-axis dominant growth are evident, and take place at temperatures indicated by Reference Lamb and HobbsLamb and Hobbs (1971). The enhanced a-axis growth rates are apparent when compared to c-axis growth rates.
Fig. 3. Facet growth vs temperature. α = 20°, ρ = 150 kg m−3, temperature gradient is 45 Km−1 and parent grain radius is 0.5 mm. The crystal growth rate changes with dominant growth direction and temperature.
Facet growth vs orientation
The faceted crystal growth rates may be integrated in time with various crystal orientations relative to the temperature gradient. In Figure 1, the total angle from the crystal edge to the temperature gradient, α, is given by
α combines the parent crystal orientation relative to temperature gradient, γ, and crystal habit, β, by locating the crystal edge relative to the temperature gradient. In Figure 1, β is always positive and less than 45° since it is measured relative to whichever crystallographic axis is dominant. Negative values of γ can better align the crystal with the temperature gradient. Simulation results for 225 hours of parent grains with 0.5 mm radius and density of 150 kg m−3 at an average temperature of 263 K and temperature gradient of 25 K m−1 are shown in Figure 4. β varied from 42.3° to 44.7° depending on the actual temperature of each crystal. Crystallographic orientations in the snowpack, γ, were varied simulating different orientations relative to the direction of the temperature gradient, α. Note that γ is dependent on the crystallographic orientation of the nucleate parent grain. Figure 4 suggests that the direction of the facet face relative to the temperature gradient is a primary factor in determining which crystals will have highest growth rates. Crystals more closely aligned with the temperature gradient grow quickly, supporting Reference Adams and MillerAdams and Miller’s (2003) proposed dominant grain growth theory discussed previously. Figure 4 demonstrates that certain crystals are preferentially oriented for higher growth rates when compared to less optimally aligned neighbors. These preferentially oriented crystals grow rapidly due to the larger excess vapor pressure created as the tip reaches warmer vapor. In addition to the facet size, the facet growth rate was examined during the same 225 hour time period. By observing the growth curve slopes, Figure 4 shows increasing growth rate with increasing size. Again, this is attributed to the increasing excess vapor pressure encountered by a crystal tip as it grows into warmer saturated vapor. It is evident that crystal orientation is an important factor, critical to determining which crystals may outpace others. While only a few select orientations were evaluated for this study, any grain orientation distribution could be implemented in the model.
Fig. 4. Facet length (l/sin(β)) vs time for various crystal orientations, with angle from the temperature gradient displayed. Smaller values of α represent closer alignment to the temperature gradient. Temperature gradient is 25 Km−1, average temperature is 263 K, parent grain radius is 0.5 mm and p = 150 kg m−3.
Facet grain-size dispersion
Data in Figure 4 demonstrate not only enhanced growth rates for preferentially orientated grains, but suggest increased grain-size dispersion with time. Initially equally sized parent grains grow at different rates due to crystal orientation relative to the temperature gradient. As a consequence, grain-size variance increases during kinetic growth.
Experimental comparison
Reference Fukuzawa and AkitayaFukuzawa and Akitaya (1993) developed kinetic crystals in the laboratory by exposing snow to very high and prolonged temperature gradients. Their experiment simulated the growth of faceted crystals near the surface of low-density snow on clear, cold nights. To further understand the relationship between temperature gradients and crystal growth, they exposed snow to gradients ranging from 100 to 300 K m−1 for up to 50 hours. Previous studies (Reference MarboutyMarbouty, 1980) had examined kinetic growth features, but at much smaller temperature gradients. Fukuzawa and Akitaya’s study was singled out for comparison because the authors were very meticulous in their experimental set-up and data collection. Additionally, snow experiencing these large gradients for relatively long periods is considered stressing for the current numerical approach.
There are no widely accepted standardized methods to quantitatively describe the geometry of kinetic growth crystals. With rounded crystals, grain radii are generally used as the primary size measure, but complex kinetic growth crystals do not lend themselves to such simple parameters. As a result, attempting to compare grain sizes between several studies is difficult. Reference Fukuzawa and AkitayaFukuzawa and Akitaya (1993) took photographs of their samples and used an image analyzer to determine grain size. They designated grain size as the diameter of a circle with the same projected surface area as the grain. For grains with predominantly symmetric planar features, this measure is similar (but not identical) to the crystal size defined in Figure 2. For simple structures such as hexagonal plates, the two size measures are similar, but as grains become highly ornate or cupped, the measurement methods diverge. The current study does not assume a circular planar area, but rather uses a crystal ‘length’ (l/sin(β) in Fig. 2b) as a measure of size which could be larger than the diameter of an equivalent circular area, particularly if crystals develop in cupped or elongated configurations. If grains become cupped or elongated, the current theory will predict a larger grain size than Reference Fukuzawa and AkitayaFukuzawa and Akitaya (1993). In addition to grain size, Fukuzawa and Akitaya calculated one standard deviation for grain-size distributions. This measurement is rarely reported, yet is highly informative. As the grains grew in their experiments, the size dispersion increased.
Even though crystal size definition is an issue for comparing results, the current model was used to simulate Fukuzawa and Akitaya’s experiments. For all experiments, Fukuzawa and Akitaya reported initial samples of very fine ice particles with an initial equivalent diameter of 0.12 mm, forming the ‘base’ of round grains for the model. The model initial kinetic crystal size is then given by l/sin(β ) where l is the radius of the parent grain. The result is a larger initial grain size than reported by Fukuzawa and Akitaya, but is necessary to account for differences in grain size definitions and the non-circular nature of Fukuzawa and Akitaya’s initial snow samples. Additionally, the current numeric approach has the nucleation site at the center node of a parent grain only, requiring a slightly larger initial facet size to extend to the vapor pore adjacent to the node. In reality, the facets will nucleate on the parent grain surface, not from the center of an existing grain. This approach simplified the computations and could be eliminated in the future by allowing latent-heat addition at locations other than grain center nodes (Reference Miller, Adams and BrownMiller and others, 2003). Facets with five crystal orientations (relative to the temperature gradient) were then allowed to grow from the base grains. For the test average temperature of 257 K, the model yields a-axis dominant growth, with β∼20°.
Figure 5 presents the model and experimental data comparisons. Fukuzawa and Akitaya’s observed average sizes and standard deviations (represented by the vertical bars) are shown. For both experiments, the average crystal size increased linearly and growth rates increased with increased temperature gradients. Additionally, the grain-size standard deviation doubled during both experiments, although the relative standard deviation remained largely unchanged. The model predictions are consistent with the experimental trends of increasing grain size with time, increased growth rates with temperature gradient, and increased size dispersion. The simulation’s increase in grain-size distribution results from varying crystal orientations, providing some physical explanation for size variation increases. The model’s crystal growth rates are constant when the crystals are relatively small, but experience increasing rates as time progresses. In the 200 Km−1 experiment, the model reasonably predicts the observed grain size and distribution. As the temperature gradient increased to 250 K m −1, the model does not track the experiment accurately. In the initial stages (first 24 hours), the model predicts the experimental results reasonably well, but starts to diverge from the experiment as the grain size increases when subjected to the larger temperature gradient. In the 200 K m−1 case, the model results would likely eventually diverge too, as the trend shows increasing growth rates with time. There are several potential explanations for the divergence in size when the model results are compared to the experiments. As the crystals get large, the experimental growth rate remains constant while the model diverges from the experiment with larger growth rates. One potential reason for the divergence in size is the size definition. Fukuzawa and Akitaya calculate an equivalent circular diameter, even if the crystals lack axial symmetry. The current method defines size as the longest crystal dimension (Fig. 2). As the crystals grow, the model allows for naturally occurring dimensional aspect ratios through the crystal axes’ relative growth defined by crystal habit; no equivalent diameter is necessary. To help illustrate the size definition differences, consider a single triangular facet face of a crystal from Figure 2b. If such a crystal was observed, Fukuzawa and Akitaya’s equivalent circular diameter would be <50% of the size predicted by the model. The model allows for one to six facet faces, further complicating the comparison. Additionally, the current approach assumed saturated vapor in all pore spaces, which in certain cases may not be true after long periods of exposure to large temperature gradients, depending upon vapor supply sources. It was also assumed that the average temperature of the kinetic crystal was the same as the parent grain from which it is growing. Due to the high thermal conductivity of ice (as compared to the pore space), this is considered reasonable until the grains become ‘large’. Since kinetic growth is taking place toward the warmer direction of the gradient, the excess vapor pressure at the crystal tip continues to increase due to its connection in the cooler region at the parent grain. The current approach only limits excess vapor pressure at the facet tip by the pore’s ability to diffuse to that location. Heat from phase change is numerically added at the center of the parent grain (Reference Miller, Adams and BrownMiller and others, 2003), not diffused through the crystal. The addition of heat conduction from the crystal tip to its root is envisioned for a future update. It also seems logical to expect some physical interference between grains growing and contacting neighboring grains, especially when they become large. The current model does not account for growth limitations created by finite pores resulting in crystal interference, contributing to the overestimation of growth. In summary, the model did predict the trends observed in Reference Fukuzawa and AkitayaFukuzawa and Akitaya (1993), but did not consistently track the average size and growth rate when grain sizes and temperature gradients became large. These differences may be attributed to size definition issues and limitations to the model based on current assumptions.
Fig. 5. Kinetic growth experiment comparison for temperature gradients of (a) 200 K m−1 and (b) 250 K m−1 at average temperature of 257 K. Vertical bars represent one standard deviation in experimental grain size and are shown around average grain size. The model shows increasing growth rates, especially for crystals more closely aligned with the temperature gradient. The model predicts increasing grain-size dispersion based on alignment with temperature gradient also. For the smaller temperature gradient (a), the model and experimental results are similar, but diverge as crystals become large (b).
Conclusions
The model agrees with fundamental established trends of kinetic growth in a dry snowpack. It reasonably predicts the crystal growth rate and size dispersion increase of a laboratory experiment for a 200 K m−1 temperature gradient, but does not accurately predict the crystal size at 250 K m−1 for the entire experiment. The results support a previously postulated dominant crystal growth theory whereby preferentially orientated grains (based on local conditions) have growth rates that outpace other orientations.
Extended applications of the model are currently being considered. For example, recent experimental research at the Montana State University Subzero Science and Engineering Research Facility has focused on the development of near-surface facets (Reference Morstad, Adams and McKittrickMorstad and others, 2007; Reference SlaughterSlaughter and others, 2009). Strong temperature gradients can develop near the snow surface, resulting in kinetic growth on and just below the surface. The current model could be used to more thoroughly examine fundamental conditions that result in these growth mechanisms.
The approach could be refined and improved. The assumed ice and pore network is very simple. The addition of increased coordination number and grain interference would bring further realism to the microstructure. Additionally, the approach could be adapted to different thermodynamic conditions such as crystal growth under ice layers or addition of sub-saturated pore spaces around faceted grains resulting in sublimation. The addition of heat conduction during kinetic growth from the crystal tip to the root is also envisioned.
The addition of kinetic growth to Reference Miller, Adams and BrownMiller and others’ (2003) unified approach produced a numerical model that uses generalized thermal conditions, allows the transfer of heat and mass to develop naturally, and simulates snow metamorphism at a microstructural level. This tool can evaluate snow metamorphism from macroscopically isothermal conditions, transition to kinetic growth and kinetic crystal growth. Using simple facet geometries, the kinetic growth model accounts for crystal habit as a function of temperature and allows for differing crystallographic orientations. The model conserves heat and mass, including phase change interactions. Ice and pore geometries are updated in time as interactions progress, resulting in a time-transient metamorphism model.
Acknowledgements
The research reported here was supported by the US National Science Foundation under grant No. EAR-0635977. The authors extend a special thanks to two anonymous reviewers who provided comprehensive and insightful evaluations of this paper. Their review significantly enhanced the paper.
