Introduction

Problems in geophysics and engineering that are affected by snow deformation have proven difficult to solve because of snow’s complex deformation mechanics and its strong dependence on scale caused by inherent spatial heterogeneity and the influence of layering. The range of applied problems influenced by snow deformation is large and includes vehicle mobility, hydrology, avalanche mechanics, ground/air energy balance and electromagnetic sensor performance (Reference MellorMellor, 1977; Reference Gray and MaleGray and Male, 1981; Reference Shapiro, Johnson, Sturm and BlaisdellShapiro and others, 1997). The scales of interest range from the microscale (10⁻³m) through the local scale (1 m) to the landscape scale (10⁴m). Landscape-scale snow properties can be derived from local-scale snow data using distributed or statistical models. However, local-scale snow properties depend strongly on heterogeneities and layering at the microscale, for which physical representations are lacking. Continuum models that treat snow at the local scale are based on the concept of a representative volume element (RVE) (Reference HillHill, 1963; Reference HashinHashin, 1964; Reference KrönerKröner, 1977; Nemat- Nasser, 1986; Reference Amieur, Hazanov, Huet, Parker and EnglandAmieur and others, 1994) in which microscale processes exhibit averaged behavior. In principle, layering and heterogeneity at the local scale can be modeled using a series of RVEs, which depend strongly on microscale (10⁻³m) snow structure and processes. In some cases, the dimensions of layers or heterogeneous features at the local scale are too small to be described by RVEs and must be analyzed at the microscale. In addition, when micromechanical processes are complex, it is difficult to develop accurate RVE descriptions of snow deformation. Consequently, the ability to accurately determine the effective parameters requires that micromechanical processes and microstructures at the micro- to RVE scales be accurately represented (Reference Savage, Herrmann, Hovi and LudingSavage, 1998; Reference Nemat-Nasser and HoriNemat-Nasser and Hori, 1999).

Discrete-particle methods that explicitly describe the dynamics of assemblies of particles and the micromechanical interaction processes between grains offer a means to improve our ability to represent microstructure and complex micromechanical processes. Reference JohnsonJohnson (1998) used a finite- element approach that included surface contacts to simulate uniaxial strain compression of an aggregate of spherical particles. The simulation qualitatively replicated the form of macroscale deformation of the aggregate while tracking the evolution in microstructure responsible for the material behavior. While this approach was able to represent microstructure and grain rearrangement, it was limited because the important grain-scale micromechanical processes could not be included in the model. To overcome this limitation, we use the discrete-element method (DEM) to develop a model of dry-snow deformation. The value of this approach is that the DEM can represent individual grains and assemblies of grains, and the micromechanics at grain contacts can be described explicitly. This explicit representation provides a way to identify the important processes that control snow deformation, directly compare physical experiments to simulations, and ensure that RVE descriptions of snow deformation accurately represent microscale processes.

Physics of Deformation in Dry Seasonal Snow

The physics that controls dry snow deformation is well understood in a qualitative sense from studies that date back to 1939 (Reference Bader, Haefeli, Bucher, Neher, Eckel and ThamsBader and others, 1954). The type of snow crystals that are initially deposited, the snow density, the state of sintering (the process of forming bonds between snow grains through the molecular diffusion of water) and the internal microstructural geometry are all known to be important (Reference Bader, Haefeli, Bucher, Neher, Eckel and ThamsBader and others, 1954; Reference BucherBucher, 1956; Reference YosidaYosida, 1956; Reference KingeryKingery, 1960; Reference KuroiwaKuroiwa, 1962; Reference Hobbs and MasonHobbs and Mason, 1964; Reference Ketcham and HobbsKetcham and Hobbs, 1969; Reference WilkinsonWilkinson, 1988; Reference ColbeckColbeck, 1997, Reference Colbeck1998).

The rate at which bonds form and metamorphose is controlled by temperature, temperature gradient, grain contact pressure, grain-size, and grain boundary geometry. Bond formation or sintering rates increase with temperature (Reference KeelerKeeler, 1969; Reference GublerGubler, 1982). Metamorphosis in the presence of small temperature gradients produces a relatively random distribution of grain bond sizes and locations, while large temperature gradients produce stiff, brittle, vertically oriented skeletal structures with large grains connected by small bonds (Reference Trabant and BensonTrabant and Benson, 1972; Reference ArmstrongArmstrong, 1980; Reference ColbeckColbeck, 1997). Melting and refreezing produce a snow cover with solid clusters of extremely large grains with large connecting bonds (Reference ColbeckColbeck and others, 1990). The relatively small grain-size and the high homologous temperatures (a normalized temperature scale with reference temperatures of absolute zero and the melting temperature of a given material) of seasonal snow initially produce a high bond growth rate that slows rapidly as the bond size increases. Molecular diffusion under equilibrium conditions (small temperature gradient conditions) moves molecules from small, high-surface-energy grains to larger grains having less surface energy. Over time, this reduces the number of snow grains and increases the sizes of grains and connecting bonds.

Low-density snow deforms primarily by grain rearrangement that occurs when stress concentrations rupture bonds or cause creep within grain contact bonds (Reference Wakahama and ŌuraWakahama, 1967; Reference Hansen and BrownHansen and Brown, 1988).

Confinement pressure has little effect on deformation, brittle failure or the yield stress of snow below a relative density of about 0.3–0.4 (Reference MellorMellor, 1975; Reference Scapozza and BarteltScapozza and Bartelt, 2003) because volumetric deformation can be accommodated within the pore space of the snow. This is a well-known behavior for many types of foam materials (Reference Gibson and AshbyGibson and Ashby, 1997). Test results that show a relationship between the confining pressure and the yield strength for low-density snow appear to be the result of modifications of the test sample introduced during experimental set-up procedures, not actual pressure dependence. These tests are typically performed by first increasing the confining pressure around the sample and then increasing the sample axial load. When the applied confining pressure exceeds the snow strength, bond rupture, grain rearrangement (densification) and an increase in grain contacts and sintering can occur, causing the snow bulk density, viscosity and strength to increase before the axial load is applied (Reference Lang and HarrisonLang and Harrison, 1995; Reference Scapozza and BarteltScapozza and Bartelt, 2003).

The number and size of grain bonds increase during natural snow densification, and movement of snow grains becomes restricted as they interfere with each other. At a relative density of about 0.6, snow reaches a density equivalent to random, close-packed particles that prevents further particle rearrangement (Reference Anderson, Benson and KingeryAnderson and Benson, 1963; Reference Ebinuma and MaenoEbinuma and Maeno, 1987; Reference GermanGerman, 1989). At this point, the dominant process of snow densification changes from rearrangement of particles by frictional shear (for loose grains) or shear deformation of ice concentrated near grain contacts (for bonded snow) to dislocation creep processes in ice with no particle rearrangement (Reference Kojima and ŌuraKojima, 1967; Reference Mellor, Smith and ŌuraMellor and Smith, 1967; Reference AlleyAlley, 1987; Reference WilkinsonWilkinson, 1988; Reference Spencer, Alley and CreytsSpencer and others, 2001; Reference Scapozza and BarteltScapozza and Bartelt, 2003). The increase in snow coordination number (the average number of contacts per grain) occurs primarily during particle rearrangement deformation. Once the critical density has been reached, dislocation processes acting within ice crystals typically produce an increase in bond size with little additional increase in coordination number (Reference Abele and GowAbele and Gow, 1975; Reference AlleyAlley, 1987).

When deformation rates exceed the ability of creep processes to relieve stress concentrations, the frequency of bond rupture increases (Reference Kinosita and ŌuraKinosita, 1967; Reference Abele and GowAbele and Gow, 1975; Reference Hansen and BrownHansen and Brown, 1987). Grain bonds rupture when stresses caused by shear, bending or stretching exceed the ice strength. Bond stresses increase due to elastic deformation, while creep deformation and grain rotation act to reduce those stresses. Intergranular friction provides resistance to movement after grain bonds rupture.

Snow deformation seldom occurs under conditions of pure shear or normal loading at grain boundaries. Off-axis loading creates moments about grain contacts that cause the grains to rotate relative to one another. Grain rotation is intimately linked to the process of deformation and must occur to allow bending or twisting stresses to develop and, subsequently, relax in grain contact bonds (Reference FaradayFaraday, 1859). Grain rotation is most prevalent at 0°C but has been observed at temperatures as low as −7°C. The frequency of occurrence decreases with temperature, which implies an increase in rotational creep viscosity as temperature decreases (Reference Nakaya and MatsumotoNakaya and Matsumoto, 1954).

A mechanistic model of snow deformation must account for the observed effects of internal snow structure at every stage of deformation and the variety of micromechanical and sintering processes that act at grain contacts. These include grain bond formation, deformation, and rupture and grain rearrangement. These processes depend on temperature, deformation rate and a grain-size distribution that is constantly evolving. The DEM addresses the complexity of the problem by providing a way to incorporate individual mechanisms into the model. The physical or structural effects can be examined and tested individually while their combined influence on snow deformation is automatically managed within the DEM construct.

The Discrete-Element Method

The DEM is a technique for explicitly modeling the dynamics of assemblies of particles. It is particularly useful when a material undergoes large-scale discontinuous deformations that depend on microscale contact processes, internal breakage of contact bonds, and compaction of broken fragments. The DEM approach allows for the use of complex particle contact physics. The DEM stores the particle shapes, velocities and locations; finds contacts; calculates forces and moments at each contact; calculates the conditions of contact bond formation, growth and rupture; and calculates the movement of each particle within the aggregate (Reference Cundall and StrackCundall and Strack, 1979; Reference HopkinsHopkins, 2004).

The Cold Regions Research and Engineering Laboratory (CRREL) microstructural snow model (µSNOW) is based on a DEM approach that uses axisymmetric particle shapes such as the cylinders with spherical ends shown in Figure 1. Other particle shapes such as ellipsoids may also be used. The cylindrical particles are formed by dilating a straight line by a sphere. In the dilation process in mathematical morphology (Reference SerraSerra, 1986), a straight line of length L is transformed into a cylinder with diameter 2R and length L + 2 R by placing the center of a sphere with radius R at every point on the line. The aspect ratio of the particle can be varied by changing L and R. The line of length L at the core of each particle is called a constraint line. Each point on the surface is a distance R from the constraint line.

Fig. 1.

Contact between two cylindrical particles and the vectors that define the contact.

Contact detection, the crux of any discrete-element code, is handled by a new iterative method (Reference HopkinsHopkins, 2004). When two particles are found to be in proximity (by standard grid methods), a vector is arbitrarily placed with its head on the constraint line of one particle and its tail on the constraint line of the other particle. This vector is modeled as an elastic band, whose ends are constrained to remain on the two constraint lines. Pulled by its elasticity, the head and tail of the vector move iteratively to locations on the constraint lines that define the shortest distance between the two particles, which is perpendicular to the external surfaces of the two particles and defines the normal to the contact surface (Fig. 1). If the length of the vector is less than the sum of the two dilating radii, the particles overlap.

Grain Contact Force Algorithms

Various contact models are used to simulate unfrozen and frozen contacts between pairs of grains and to govern when and how quickly freezing occurs. The contact models operate on a circular contact plane centered at the point of contact at X in Figure 1. The assumption of a circular contact plane is a simplification of the many possible complex contact geometries that depend on particle shape and contact orientation (Reference JohnsonJohnson, 1987). The unfrozen- contact model simulates viscous–elastic compressive forces in the direction normal to the plane of contact, and frictional sliding in the plane of contact. The frozen-contact model is also viscous–elastic and simulates tensile and compressive forces in the direction normal to the plane of contact, shear forces in the plane of contact, a twisting moment about the normal to the plane of contact, and a bending moment about an axis in the plane of contact. The frozen-contact model permits viscous creep to operate on all deformation mechanisms, and brittle fracture to operate on the shear, tension, bending and twisting mechanisms; failure in compression is through power-law creep.

The compressive force model operates in the direction normal to the contact plane. We use a simple linear viscous–elastic normal contact force model in which the overlap between particles is interpreted as elastic compression. Elastic deformation of unfrozen contacts is damped by combining the elastic force with a linear viscous element in parallel. The normal force component for unfrozen contacts is

where δ is the overlap, A _c is the unfrozen-contact area, and the contact stiffness is k _ne ¼ E/2 R, where E is the material Young’s modulus. The contact viscosity k _nv is a fraction of the critical damping value calculated from k _ne and the particle mass, and acts to attenuate oscillations, helping to maintain stability. The superscript i denotes the current timestep. The subscripts n, e and v denote the normal direction, elastic component and viscous component, respectively.

Friction between grains provides shear resistance. The frictional force model loads incrementally through the relative velocity between particles in the plane of contact until the force reaches the Coulomb limit, whereupon the particles begin to slide. The elastic component of the frictional force is

where k _te is the tangential elastic stiffness, Δt is the simulation time-step, and is the tangential component of the relative velocity of particle 1 with respect to particle 2 at the contact point, including rotational components. The tangential force has a viscous component provided by a dashpot in parallel with the spring (Reference WaltonWalton, 1980). We use a ratio of tangential to normal stiffness of 0.6. The magnitude of the combined tangential force is scaled so that it does not exceed the Coulomb limit where µ is the friction coefficient and is the combined viscous–elastic tangential force.

Frozen-Contact Force Model

In frozen contacts, deformation is measured by the relative displacement between circular contact patches tangent to the surface of each particle. The circular contact patches are connected by a viscous–elastic glue. Initially, with the contact in an undeformed state, the two patches are coincident (that is, the normal and tangential vectors that define the contact patches for the particles are coincident). The locations of the frozen-contact planes are stored in the body coordinate frame of each particle and move with them. As the particles move relative to each other, the glue between the patches is deformed, creating forces and moments. Four modes of deformation are modeled: tension and compression in the direction normal to the contact plane; shear in the plane of contact; twisting about the normal to the plane of contact; and bending about an axis in the plane of contact.

A close-up view of the circular contact patches, with radii r_b, in a deformed state is shown in Figure 2. Deformation causes the contact patch of one of the particles to displace and/or rotate relative to the other, causing their descriptive vectors to change as well. The vector ñ₀, which has a magnitude of zero in the undeformed state, represents the tensile displacement (stretch) between the centers of contact patches 1 and 2. Bending, which takes place about the axis in the plane of contact, is described by the angle α_n between , the unit normal to contact patch 1 and the unit normal to contact patch 2. Twisting, which takes place about the normal to the plane of contact, is determined by the rotation angle, α_t , between the unit vectors that lie in contact patches 1 and 2, respectively.

Fig. 2.

A close-up view of the frozen-contact point, the two contact patches, the vector that connects the center of the two contact patches, the normal vectors and the angle α _n, and inplane unit vectors and and the rotation angle α_t between and . (The projection of onto the contact patch for particle 1 is shown by a dashed line vector.)

The normal component of the elastic force is

where is the area of the frozen-contact patch and δ_nf is the normal component of the average strain defined as The subscript f refers to the frozen contact, and is the unit normal, to the contact patch on particle 1 in Figure 2. The superscript i and other subscripts are as above. The tangential component of the elastic force is

where is the shear component of the average strain where The normal and tangential components of the viscous force are

where is the relative velocity between particles 1 and 2.

If we imagine contact patches 1 and 2 in Figure 2 are connected by viscous–elastic fibers, the bending described above creates regions of tension and compression between the contact patches. We calculate the moments produced by this distribution of forces by integrating the forces and their moment arms produced by deformations in terms of the bending angle α _n, the rate of change of and the contact patch radius r _b:

Similarly, twisting about the axis normal to contact patch 1 creates a distribution of shear forces between the contact patches. We calculate the moments produced by this distribution of forces by integrating the forces and their moment arms produced by deformations in terms of the twisting angle α _t, the rate of change of and the contact patch radius r _b:

where r _b is the radius of the contact patch, is the angle between the normal vectors to the frozen-contact planes (Fig. 2), and is the angle of rotation of one plane relative to the other about the normal, defined in terms of and , unit vectors in the respective contact patch. Vectors and are initially coincident. The forces and moments are summed to find the resultant on each particle.

Frozen-Contact Fracture Model

Bonds fail by brittle rupture when the tensile stress in the outer fibers of the bond equals the temperature-dependent tensile stress of ice. Until microscale ice-strength data become available, we determine ice strength by a fit to available strength tests on polycrystalline ice, given by T _ice ¼ 0.99648–0.045904T, where T is temperature in degrees Celsius and stress is in MPa (Reference ButkovichButkovich, 1954; Reference Hawkes and MellorHawkes and Mellor, 1972; Reference HaynesHaynes, 1978; Reference Cuda and AshCuda and Ash, 1984). The influence of grain-size on the ice tensile strength was not taken into account, but it is less than the scatter of test data used in our curve fit and is a secondary influence to the temperature dependence (Reference Currier and SchulsonCurrier and Schulson, 1982).

Tensile stresses develop in tensile and bending modes of loading, and apparent tension develops during shear and twisting from the elastic tangential forces described in Equations (4) and (9). Brittle failure of a frozen contact is modeled using the rate-independent strain-softening constitutive model shown in Figure 3. In the figure, the stress at a point follows the initial loading path that has a slope of k ₁ = k _ne. When the stress at a point reaches the tensile limit, σ_t = σ_ice, corresponding to a strain of δ ₁, further strain causes the stress to decrease along the unloading path with slope k ₂. The stress during unloading is σ = σ_ice – k ₂(δ – δ ₁). When the strain reaches δ₂, the joint is broken at that point. Rigorously implementing this model over the circular frozen-joint contact patch would require specifying the strain as a function of position and separately integrating the stress over the loading, unloading and broken regions of the circular contact area. This would yield much more complicated expressions for the force and moments than those in Equations (3–10). In the constitutive model in Figure 3, the material is weakening and failing while the area of the contact remains constant. Alternatively, we can hold the material strength constant and reduce the area by allowing the bond radius r _b and the contact area to decrease to produce a weakening of the bond that is at least qualitatively equivalent to the unloading path in Figure 3. We can rewrite the stress equation above as a force equation in terms of the bond area that can be solved for the new bond radius r _b, given by

Fig. 3.

Constitutive model that defines the tensile stress at a point, σ, in terms of the strain, δ.

where r _b0 is the bond radius at the onset of failure. The two approaches are equivalent when the strain δ is uniform over the circular contact patch, as it is in the case of the tensile and shear failure modes expressed in terms of average normal and tangential strains δ _n and δ _t. However, when there are non-zero bending and twisting moments, the strain is not uniform over the contact plane and the approximation is not exact. Nevertheless, the simple approach preserves the strain-softening failure of the constitutive model in Figure 3 while greatly simplifying the integration.

To reiterate, failure begins when the strain exceeds δ₁ in Figure 3, and continues until the strain exceeds δ₂, whereupon the bond radius r _b = 0 and the joint is broken. To bring the various strain components into the failure calculation, the maximum in-plane and out-of-plane strains are calculated, and the larger of the two is substituted for δ₁ in Equation (11) to calculate r _b. The maximum out-of-plane and in-plane strains are, respectively,

Combining all of the strain components to obtain the maximum combined strain on the contact plane is too difficult. At present there is no compressive failure mechanism.

Grain Boundary and Power-Law Creep Models

The idea behind the grain boundary creep model is that creep on the grain boundary of a frozen contact acts to relax the stress in the bond. Power-law creep exponents of 1–1.8 have been observed for grain boundary sliding (Reference Langdon, Whalley, Jones and GoldLangdon, 1973; Reference Goldsby and KohlstedtGoldsby and Kohlstedt, 1997). For the current model, we treat grain boundary sliding as Newtonian with a contact viscosity v _s, which is treated as a free parameter.

When the frozen bond is created, a circular contact patch is placed tangent to the surface of each particle. The locations of the contact patches are stored in the body coordinate frame of each particle. Relative motion between the particles moves the patches out of their initial alignment, creating the forces and moments described above. The creep process acts by bringing the circular contact patches into alignment, reducing the forces and moments.

The creep induced by the elastic component of the shear force F_te (Equation (4)) is in the direction of the in-plane component of the vector that connects the centers of the contact patches. The plane of contact is arbitrarily assumed to be the circular contact patch attached to particle 1. The shear creep displacement of the contact patch Γd ₁ is

Two grains in frozen contact will creep rotationally in response to a bending moment at the contact. This behavior is well documented experimentally (Reference FaradayFaraday, 1859; Nakaya and Matsumoto, 1954). The creep displacement produced by the bending moment M _ne (Equation (7)) is in the direction of the difference in the two normal vectors, The creep displacement Γd ₂ produced by the bending moment is

The displacements Γd1 and Γd2 , with the appropriate direction vectors, are used to calculate the new positions of the contact patches attached to each particle. The twisting moment M _te (Equation (9)) produces creep in the direction of the difference in the in-plane vectors attached to each particle. The twisting creep displacement Γd ₃ is

The displacement Γd ₃, with the appropriate direction vectors, is used to calculate the new position of the inplane vectors in the contact patch attached to each particle.

Power-law creep can act in shear in the grain contact region or in compression. At this stage of model development, shear power-law creep is not implemented, as it is considered to be secondary to grain boundary creep at the stress levels of interest. The compressive power-law creep induced by the normal force F _ne (Equation (3)) acts in a direction normal to the plane of contact. The only way to allow creep to relax the normal force is to allow the frozen- contact point on each particle to migrate into or outside the respective particle. This can be done by allowing the contact point to be more or less than the dilation radius away from the constraint line. To accomplish this, we introduce two new variables, that function as effective dilation radii. The creep displacement Γd ₄ produced by the normal force is

where m is the power-law creep exponent. The displacement Γd ₄, with the appropriate sign, is added to but their magnitudes are limited so that (1/2)R > R′ < (3/2)R for practical reasons. The radius r _b and the area A _f of the frozen contact change as well.

Sintering and Particle Contact Bond Growth

Experimental and theoretical efforts to determine the dominant sintering mechanisms in snow have yet to yield definitive results. Reference KuczynskiKuczynski (1949) proposed that sintering mechanisms could be identified by observing the rate of sintering. Using this approach, Reference KingeryKingery (1960) concluded that the dominant mechanism for sintering of ice grains was surface diffusion. Reference KuroiwaKuroiwa (1962) analyzed his experimental data using Kuczynski’s approach and concluded that volume diffusion was the dominant process. Reference Hobbs and MasonHobbs and Mason (1964) conducted experiments that they interpreted as pointing to the sublimation, transport and deposition of water molecules through the vapor phase as the dominant sintering mechanism. Reference ColbeckColbeck (1997, Reference Colbeck1998) reconsidered the experimental interpretations based on Kuczynski’s approach. In his review of sintering, Reference ColbeckColbeck (1997) suggested the possibility that different mechanisms may dominate at different stages of the sintering process, and described a sintering theory developed from crystalline materials as possibly applying to snow. Later, Reference ColbeckColbeck (1998) observed that the existence of grain boundary grooves at grain bonds implies a stronger role for grain boundary diffusion than previously thought. He cited the work of Reference Zhang and SchneibelZhang and Schneibel (1995), who described the sintering of grains joined by grain boundary grooves to support the hypothesis that sintering is controlled by the grain boundary groove angle and that grain boundary diffusion is the dominant mechanism. When grains first come into contact, the grain boundary groove angle will generally be significantly different than the equilibrium dihedral angle (145° for pure ice and vapor), so the initial rate of growth for the bond radius will be relatively high. As the bond groove angle approaches the value of the dihedral angle, the growth rate will decrease. Once the bond groove angle equals the equilibrium dihedral angle, any further increase in grain contact bond radius requires that the grain-sizes increase proportionally to maintain the equilibrium dihedral angle at the grain boundary groove (Reference GermanGerman, 1996).

A further problem is that the observed rates of grain bond growth in laboratory experiments are less than those observed in natural snow (Reference KeelerKeeler, 1969). Reference ColbeckColbeck (1997) speculated that the increased bond growth rate in nature is due to temperature gradients in the snow cover, driving water vapor through the snow. The increased vapor transport is thought to cause faster grain growth, hence faster bond growth between grains.

At present, we use available empirical data and formulations to treat sintering because of the lack of a definitive physical understanding of the sintering process for snow. Our empirical sintering formulation uses the results of Reference GublerGubler’s (1982) measurements of grain contact bond strength as a function of contact time, temperature and contact force. By assuming the strength of the bond contact is equal to the tensile strength of ice, it is possible to calculate a rate of bond growth from Gubler’s results. This is done by using Reference WilkinsonWilkinson’s (1978) power-law creep model of sintering of two particles (Fig. 4) to treat the effects of pressure sintering, and Gubler’s results, with the effects of contact force removed, to treat the effects of pressureless sintering. The rate of change of the grain bond radius from sintering, based on the particle contact geometry shown in Figure 4, is (Johnson and Hopkins, unpublished internal CRREL report, ‘A discrete element model of snow deformation for dry snow’, 2004)

Fig. 4.

Ice particle bond contact geometry (Reference WilkinsonWilkinson, 1978).

where F _c is the compressive contact force between particles, T_K is temperature in Kelvin, A _{0 ice} = 0.00554458 s⁻¹ Pa⁻³, and Q = 70 × 10³ J M⁻¹ is taken as the most representative value over a temperature range of −0.2 to −20°C for single-crystal ice deformation (Jones and Brunet, 1985). The power-law coefficient m for singlecrystal ice is temperature-dependent, with a mean value of m = 2, given by m = 1.98488–0.00607767T, determined from a curve fit to data from Reference Jones and BrunetJones and Brunet (1978). The radius of sintered ice grains as a function of the contact duration t, determined from Reference GublerGubler’s (1982) data, is

where F _t is the tensile rupture force for two ice grains in contact, time t is in seconds, q = 1 s^–β, and α and β are functions of temperature and geometry, determined by a best fit to Reference GublerGubler’s (1982) data, given by

The first term in the numerator of Equation (17) is derived from Reference GublerGubler’s (1982) data and includes the influence of the contact pressure. The second term in the equation is the derived contribution from the loading pressure used in Reference GublerGubler’s (1982) experiments, and the third term adds the contribution to any applied force normal to the grain contact. The complete derivation of the sintering model is given in Johnson and Hopkins (unpublished internal CRREL report, ‘A discrete element model of snow deformation for dry snow’, 2004).

DEM Simulation of Snow Creep Settlement

To demonstrate the capabilities of the µSNOW model, we performed two series of simulations of snow creep settlement under gravity loading. The first series of four simulations was done using grain contact viscosities in the range 0.0016–1 MPa s and was used to examine the overall performance of the model (Table 1). The second series of simulations used grain contact viscosities that bracketed those of the first test series, to examine the ability to scale creep simulations as a way to improve computational efficiency.

Table 1.

Parameters used in the settlement simulations

A snow sample containing 1000 particles was constructed with randomly oriented particles surrounded by walls on the x and z boundaries (Fig. 5). The sample was constructed by placing particles with random size and orientation in a low- density array at the nodes of a three-dimensional cubic grid. Following a method developed by Reference Cundall and StrackCundall and Strack (1979), space was contracted until the particles reached the desired density. To simplify the simulations, the temperature- and pressure-dependent sintering model described above was turned off. Instead, whenever a particle contacted its neighbor, a frozen bond was instantaneously created with a randomly chosen contact area uniformly distributed between 0.03125πR ² and 0.09375πR ². Creep mechanisms were turned on to relieve stress concentrations. The y boundary was open and periodic, that is to say, the particles at one y boundary interacted with the particles on the opposite boundary. The parameters used in the first series of simulations are listed in Table 1. We placed a 0.125 kg mass on the lid and gradually increased gravity from zero to the final scaled value in Table 1.

Fig. 5.

Initial configuration of the settlement box simulation. The box contains 1000 particles. The x, y and z dimensions of the box are 0.019, 0.016 and 0.019 m, respectively. The lid at the top of the box is free to move vertically.

Fig. 6.

Evolution of test sample stress and particle contact conditions during creep settlement for four grain contact viscosities as functions of time. (a) Vertical lid bulk low stress, P_zz , and sidewall bulk stress, P_xx ; (b) snow density; (c) grain coordination number; (d) the number of broken grain bonds; and (e) the number of new frozen bonds.

The model has four important parameters that govern settlement. Three are material parameters – contact stiffness k _ne, tensile strength σ_ice and contact creep viscosity ν _s – while the fourth parameter is the weight of the lid on the box. The contact stiffness, defined above, is the modulus divided by a particle length scale over which the elastic contact deformation is assumed to take place. The simulation time-step is proportional to the square root of the particle mass divided by effective stiffness. For realistic contact stiffness and particle mass for a millimeter-diameter particle, the time-step is about 3.3 × 10⁻⁷ s, which produces excessively long calculation times. To speed up the simulation, we reduced stiffness, tensile strength and load together to increase the time-step. A reduction in stiffness by a factor of 1000 allows the time-step to be increased to 1 × 10⁻⁵s, reducing the simulation run time by a factor of about 30. The speed of the simulation is about 10⁵ particle time-steps s⁻¹ on a 3 GHz Pentium 4. Therefore, a 1 s simulation using 1000 particles requires about 1000 s.

The simulations were conducted to examine the controlling factors for snow settlement, including the effects of internal structure, coordination number and grain contact creep viscosity. The simulations began with the initial configuration shown in Figure 5. The grain contact viscosity coefficient ν _s was systematically increased from 0.0016 MPa s to 1 MPas in four steps. The loading stress P _zz , the confining stress P_xx (equal to P_yy ), the bulk density, the coordination number, the number of broken bonds and the number of newly frozen bonds are plotted in Figure 6 as a function of time. The loading stress P_zz and the confining stress P_xx exhibit a rapid rise, followed by variations uncorrelated with the contact viscosity (Fig. 6a). The confining stresses, P_xx and P_yy , are due to interparticle forces acting on the lateral walls and across the periodic boundaries. The rate of increase of density in Figure 6b (equivalent to the settlement rate) increased with decreasing contact viscosity and gradually slowed with time as the coordination number and number of frozen bonds increased (Fig. 6c and e). Most bond breakage occurred in the initial loading period (Fig. 6d).

In Figure 7 the results from Figure 6 are recast as functions of density. Since the duration time of each simulation is the same and the settlement rate varies inversely with the contact viscosity, the change in density during each simulation decreases as contact viscosity increases. Within the density range of the simulations, the coordination number vs density curves (Fig. 7b) and the frozen-bond contact number vs density curves (Fig. 7c) for the four simulations are essentially the same. The relationships between the number of broken bonds and density for each of the simulations are similar and appear to be primarily associated with the process of placing the load on the samples at the start of simulations (Fig. 7a and d). The increase in the coordination number appears to be a linear function of density, but this number should eventually decrease at higher densities as particle rearrangement slows. As settlement progresses, the confining stress remains nearly constant over the entire density range. This implies that, at these densities, the available pore space is sufficient to accommodate particle rearrangement during creep. The functional similarities of the coordination number and frozen-bond number curves imply that particle rearrangement paths in each simulation are the same irrespective of the contact viscosity magnitude.

Fig. 7.

Test sample stress and particle contact conditions as a function of density for four grain contact viscosities as functions of snow density. (a) Vertical lid bulk stress, P_zz , and side-wall bulk stress, P_xx ; (b) grain coordination number; (c) the number of new frozen grain bonds; and (d) the number of broken bonds.

The microscale creep and bulk creep for the model snow, calculated during the settlement simulations, are shown in Figure 8. The average creep rates for each of the three micromechanical creep mechanisms that operate are shown separately in Figure 8a–c. Shear creep is the movement of the frozen contact in response to the contact shear force described in Equation (13). Bending creep describes rotation of the grains about an axis in the plane of contact in response to a bending moment (Equation (14)), a mechanism first identified by Reference FaradayFaraday (1859). Twisting creep occurs when grains rotate about the axis normal to the plane of particle contact in response to a twisting moment (Equation (15)). Twisting creep is an important, previously unrecognized mechanism that is necessary to allow particles to reconfigure themselves in response to contact moments. In general, if one or two creep mechanisms are turned off or missing, creep taking place through the remaining operational mechanisms causes stress to increase in the modes of deformation whose creep mechanisms are missing. This eventually produces an explosive, unrealistic sample failure.

Fig. 8.

Micro- and bulk creep as a function of density for four grain contact viscosities. (a, b) The mean creep displacement per grain contact point for all grains in shear (a) and bending (b); (c) the rotation of grains around their contact points with other grains without displacement; and (d) the bulk snow axial strain rate.

The bulk creep rate or lid settlement rate is shown in Figure 8d. The bulk creep rate initially increases as the sample equilibrates to the weight on the lid. It then decreases slowly with increasing density as grains pack together and new frozen contacts form, exhibiting the same general trend shown by the individual creep mechanisms.

The bulk (Newtonian) viscosity, obtained by dividing P_zz by the bulk creep rate, is shown in Figure 9 for the original four simulations, along with three additional simulation results at higher values of contact viscosity. This allows us to further investigate how well creep problems can be scaled. The bulk viscosity for all simulations increases with the contact viscosity. The bulk viscosity also increases with density as the number of grain contacts increases, shown by the coordination number (Fig. 7b) and the number of frozen bonds (Fig. 7c). The range of measured bulk viscosities for natural snow is shown at the top of Figure 9.

Fig. 9.

Bulk viscosity as a function of density for different grain contact viscosities compared to the bulk viscosity for natural snow (Reference Kojima and ŌuraKojima, 1967).

The discrete values of the bulk viscosity obtained using contact viscosities of 1 MPa s and greater were determined in the same way as were the continuous values of bulk viscosity using low contact viscosities, i.e. by dividing the bulk creep stress by the bulk creep strain. Conducting continuous simulations at these high contact viscosities, as was done for the low-viscosity simulations, would have required extremely long run times. We reduced run times significantly by creating the starting samples for the discrete simulations at specified densities using intermediate stages of deformation from the continuous simulation that used a contact viscosity of 0.0016 MPa s. The run times were further reduced by conducting the simulation on these samples for only the length of time needed to establish the creep deformation trends. Thus, only a limited number of points (discrete) for bulk viscosity at high contact viscosities were used in Figure 9 to determine the bulk viscosity trend lines.

The separation between the bulk viscosity curves in Figure 9 is a function of the grain contact viscosity, while the slope of the bulk viscosity curves is a function of the grainscale creep deformational mechanisms (Equations (13–15)) that operate during densification. In the model, all three creep processes depend on the same grain contact viscosity.

The parallel bulk viscosity vs density slopes for the different contact viscosities shown in Figure 9 are an indication that particle rearrangement paths, and the relative interaction of the three creep mechanisms, are unaffected by parameter scaling. As a result, it may be possible to simulate creep settlement problems that are controlled by the three linear creep mechanisms with reasonable computational time by carefully scaling snow material and loading parameters. Further work will be required to determine the ability to scale high-density snow creep when the power-law creep mechanism (Equation (16)) becomes more important.

Our ability to select an effective grain contact viscosity to simulate the bulk viscosities vs density implies that microscale contact viscosities can be estimated from simulation of macroscale experiments. However, the ability to make such an estimate is limited by the fact that the effective contact viscosity (Equations (13–16)) is determined by both the contact viscosity and the bond area. Further complications arise because of the influence of snow structure geometry and the mean number of contacts per grain. This means that microscale tests will be needed to fully define the micromechanical properties used in the model.

A linear scaling relationship between the grain contact viscosity and bulk viscosity for a snow density of 300 kg m⁻³ was obtained by scaling the contact stiffness k_ne, the tensile strength σ_ice, the contact viscosity v_s, and the lid mass by successive orders of magnitude (Fig. 10). This relationship is empirical and valid only for the mean particle contact area used in our simulations. For example, a larger mean particle contact area would result in lower contact viscosities to achieve the same bulk viscosity. The influence of non-linear, power-law creep on the scaling relationship was not studied but needs to be investigated to determine the difficulties associated with using scaled parameters to reduce the computation time for snow creep problems.

Fig. 10.

Calculated bulk snow viscosity for a snow density of 300 kg m⁻³ as a function of the grain contact viscosity. The open marker is the predicted contact viscosity required to achieve the same bulk viscosity as natural snow.

Conclusions

The µSNOW model, developed using a new discrete-element technique that employs grain-scale force models and an explicit snow microstructure composed of assemblies of grains, provides a physically accurate representation of the mechanisms that control dry-snow deformation. The deformation behavior of the model is determined by the distribution of the axisymmetric particle sizes and orientations, the distribution of contact bond areas, and the magnitude of the contact viscosity.

Simulations of creep settlement using 1000-grain model snow samples and scaled values of contact viscosity are able to replicate the observed dependence of the snow bulk viscosity on density. In the simulations, the dependence of snow bulk viscosity on density is caused by an increase in coordination number (which causes the number of frozen bonds to increase), the mean area of particle contact bonds, and the viscosity at the contacts. The dominant micromechanical deformation processes at densities below 420 kg m⁻³ were particle rearrangement through shear creep at particle contacts, grain rotation (bending), and grain twisting around grain contact points using a linear contact viscosity. Power-law creep should become important as the snow grain packing restricts the ability of grains to rearrange themselves.

Many problems of interest in snow geophysics and engineering at micro- and RVE scales can be addressed directly using the model. As computational capabilities improve, the model can be used to determine local-scale snow properties based on accurate microscale processes. The results of the simulations can be used to test and develop the continuum descriptions of snow behavior needed to treat landscape-scale problems.

Simulation of snow creep deformation, using realistic snow grain stiffness values, produces excessively long calculation times. The current simulations indicate that for linear creep mechanisms, the contact viscosity can be scaled linearly, which allows creep problems to be simulated in much shorter times than using actual snow creep viscosities. Further work is needed to determine if non-linear mechanisms, such as power-law creep, can be similarly scaled.

Quantitative values for many of the parameters used in the µSNOW model are lacking because of a dearth of high-quality microscale experiments on snow deformation and sintering processes. In some cases, these parameters can be estimated by comparing the simulated bulk mechanical properties of snow with measured values, but these are effective values only. For example, the creep of grain contacts will depend on both the contact viscosity and the contact bond radius. In our simulations we set the range of contact bond radii and varied the contact viscosity. We could just as easily have held contact viscosity constant and varied contact bond radius. This means that to accurately determine the microscale parameters needed in snow simulations, grain-scale experimental measurements of deformation, failure and sintering over a range of temperatures, grain-sizes and loadings are needed. In addition to micromechanical data, information about snow microstructure from both laboratory and in situ measurements is needed to construct realistic model snow samples and to assess model snow-structure evolution and metamorphosis.