Introduction

Snow is a granular geologic material with ice as the matrix material. The particles have varying shapes, crystallographic orientations and numbers of bonds connecting each particle to neighboring particles. In its natural state, snow exists at close to its melting point and can readily undergo a variety of metamorphic processes including sintering and bonding of grains, temperature-gradient metamorphism, melt/freeze processes, and heat and mass transport due to inhomogeneities in physical properties. The process leading to the formation of a rounded equilibrium form is called equitemperature metamorphism, for it proceeds in bodies which are not far from a uniform temperature. The overall strength of snow is also affected during this metamorphic process. The deposition of ice at contact points produces bonds or necks between adjoining grains in a process known as sintering. As the process proceeds, the necks become larger and the snowpack strengthens.

An externally applied pressure serves to increase the rate of sintering. The mechanism of pressure sintering is analogous to creep kinetics for polycrystalline metals. A number of processes, including lattice diffusion, grain-boundary diffusion and dislocation creep, have been used to explain the phenomenon. The regime of dominance of each of the above mechanisms depends upon the snow microstructure, the applied pressure and temperature. Reference Maeno and EbinumaMaeno and Ebinuma (1983) used pressure-sintering diagrams to show these various regimes in snow samples with densities varying from 600 to 900 kg m⁻³. Depending on the radius of the particles, at very low pressure lattice diffusion or boundary diffusion may be the predominant mechanism for sintering. At higher pressures dislocation creep takes over. However, if pressure is very high, it is likely to lead to neck fracture, resulting in interparticle slip.

In this study, snow is modeled as an assembly of spherical particles joined by areas of much smaller cross-section called bonds or necks. Because of their smaller cross-section, the necks are subjected to much higher stresses than are the ice particles and therefore undergo large deformations. It is the deformation of necks which is responsible for deformation of snow with densities in excess of 250 kg m⁻³, particularly at low stresses. At higher stresses, the necks may fracture and interparticle slip becomes a significant deformation mechanism. It is this relative displacement between the particles which is largely responsible for strains at large stresses.

In this study we attempt to determine the relative contribution of each of these mechanisms (dislocation strain, superplastic strain and interparticle sliding) to the overall deformation of snow. We will restrict ourselves to medium- to high-density snow and therefore utilize an approach that does not need to consider chains of ice grains as done by Reference KryKry (1975). For snow with densities greater than 250 kg m⁻³, the coordination number (number of bonds per grain) is equal to or larger than the value of 3. Therefore, the role of chains is minimized relative to deformation taking place in the necks. Modeling all of the processes mentioned above and incorporating them into a macroscopic constitutive relation is not an easy task. The results we present here are not final but do provide a measure of the potential of constitutive equations developed in terms of microstructural processes.

Constitutive Formulation

Reference OdaOda (1974) and Reference Kanatani, Jenkins and SatakeKanatani (1983) used probability principles to determine stresses at contact points of sand grains. The same principles can also be applied to ice particles in snow. Suppose the material, composed of rigid spheres, is subject to a macroscopically uniform stress. The contact forces vary from particle to particle and with the position of the contact on the grain surface. All contact forces are superposed on a hypothetical representative particle whose radius is the average particle radius.

If the number of particles is sufficiently large, the contact-force distribution on the representative particle is approximated by a continuous function, D(n), of the contact direction determined by n, the outward unit normal vector at the contact point. The direction of n is determined by the two angles, α and β, shown in Figure 1. The unit normal, n, is actually the same vector as the vector ν ₃ shown in Figure 1. Since D(n)dn is the number of contact points on all particles in the differential solid angle, dn, divided by the number of particles, we see that

(1)

where Ν is the coordination number (number of contacts per grain). Let Ν. f_i (n)D(n)dn be the components of the total force acting in the differential solid angle, dn, then by equilibrium of force and moments, one can show that

(2)

(3)

where f[ _i (n)n_j ] is the skew part of f_i (n)n_j ]. In other words, the matrix of components of f[ _i (n)n_j ] are just f[ _i (n)n_j = (f _in_j − f_jn_i /2. This last relation results directly from requiring that the forces applied to the grain satisfy equilibrium of moments about the grain center. f( _i (n)n_j ) = (f_in_j + f_jn_i )/2 is the symmetric part of f_i (n)n_j .

Fig. 1. Coordinate system. x_i global coordinate system entered at the ice particle center. The x_iʹ neck coordinate system is centered at the center of neck.

Reference Kanatani, Jenkins and SatakeKanatani (1983) imposed a virtual displacement which distorts the spherical particles into an ellipsoid. Even if the particles are assumed to be rigid, virtual deformations are hypothetical motions and can be imposed (Reference Kanatani, Jenkins and SatakeKanatani, 1981). The virtual work done by the contact forces on the representative particle is equated to the virtual work per unit volume done on the virtual strain due to the applied stress ty. After some algebra, this yields

(4)

where a is the particle radius and γ is solid volume fraction. We see here that the stress tensor t_ij is determined only in terms of the symmetric part of f_in_j , where we have dropped the n in these terms for the sake of brevity. The contact-force density f_iD can be expanded into a series of spherical harmonics and truncated to the first two terms since |n_i | ≤ 1:

(5)

The term B _[ij] is just the skew part of the matrix B_ij . Substitution of this form into Equations (2) and (3) results with

(6)

and Equation (5) reduces to

(7)

Substituting Equation (7) into Equation (4), one obtains

(8)

(9)

Here f_i is the force per contact on the contact with normal n_i . In the development presented later this is replaced by (σⁿ _i · Α_TΝ) where σⁿ _i is the stress vector acting on a grain at the contacts, A _T is the total area of the necks or contacts, and Ν is the coordination number. A _T/Ν is therefore the average neck cross-sectional area. In the above the superscript “n” is used to denote the stress vector acting on a grain surface with a unit outward normal vector n. Subscripts are used to denote coordinate components of vectors and tensors.

The distribution function, D, can also be written as D = Ρ · Ν, where Ρ is the probability that two adjacent particles will form a contact. For the isotropic case, Ρ is independent of α and β and therefore has the value Ρ = 1/4π. Substituting this into Equation (9), we get in vector notation

(10)

This relates the stress vector on the ice grains to the macroscopic stress vector in the snow.

A local coordinate system, x_i ʹ, can be set up at the contact point and the stress vector at the contact point or the ice neck can be resolved into three components along the coordinate axis of this local coordinate system. This coordinate system is illustrated in Figure 1. The unit vectors ν _i are defined relative to the spherical coordinates as shown if Figure 1. The vector ν₃ is the same vector as n.

We now denote the components of stress on the face with normal n

(11)

These components act in the directions of ν _i , i.e. the directions of the x_i ʹ coordinates. The component σ₃₃, for instance, is the normal stress acting perpendicular to the grain surface and in the direction parallel to the neck center line. The average value of σ₃₃ over a section of the grain surface between the angles α and β to α + Δα and β + Δβ can be shown to be

(12)

Similar expressions can be found for the other two components, σ13 and σ23.

In Equation (12) the factor a ²/γA _T reflects the influence of the radius of the neck, density ratio and area of the ice particle. Rather than treating all of these as separate variables, this equation shows that they can be grouped as one parameter. Snow samples having different densities, grain and neck radii can still end up having the same stress at the contact point, provided a ²/γA _T is the same for the different samples. This term is a fundamental dimensionless parameter relating the microstructure to many properties of the material.

The strain in snow is primarily due either to the strain in the necks or to relative sliding of particles with respect to each other once the necks have been broken. Experiments (Reference KryKry, 1975; Reference BrownBrown, 1980) show that snow does not collapse even when the stresses are high enough to cause fracture of some of the necks and subsequent intergranular sliding. From this it was concluded that even when sliding of necks is taking place, the original skeletal structure still remains largely intact. We now determine the different deformation mechanisms taking place in the necks under varying load conditions in order to incorporate the mechanisms into a constitutive formulation for snow.

Inelastic Deformation of Necks

At stresses lower than 0.004 MPa, insignificant neck fracturing can be assumed. In this case the deformation of snow can be attributed to deformation of intergranular necks, as has been suggested earlier. For a stress of 0.004 MPa applied to snow, the principal stresses in the necks are less than 0.7 MPa, as can be shown by using Equation (12). In ice it has been observed that crack nucleation does not take place until the stress has exceeded a limiting value of 0.5 MPa (Reference GoldGold, 1972). Reference SinhaSinha (1984) cited the experimental results in which specimens were seen to last for 4 d without fracturing. He has conjectured the possibility of existence of super-plasticity in ice due to favorable conditions such as high temperature, low stress and fine grain-size.

If the stresses in the neck are of a moderate value, superplastic deformation is not significant, and ordinary constitutive laws should be applicable for modeling both steady-state and transient deformation of ice. Deformation at this intermediate stress level is due primarily to dislocation processes. When strains in the neck reach a critical value, fracturing can begin to occur. Thus, in a snow sample subjected to high stress levels there may be necks which are oriented such that they are subjected to high stresses and undergo fracturing and subsequent sliding. On the other hand, there may be suitably oriented necks which are subjected to a low stress and therefore undergo either superplastic deformation or dislocation-dominated creep.

To describe the behavior of ice necks under combined compressive and shear stress, a multi-axial constitutive law is required. A number of models exist to describe this behavior. The model used in this paper is based on work of Reference Szyszkowski and GlocknerSzyszkowski and Glockner (1986, 1987). In this model ice is treated as an isotropic, non-linear viscoelastic material. The heredity effects have been included with the use of a Volterra integral in the following equation of Szyszkowski and Glockner:

(13)

where

(14)

ν is Poisson’s ratio, E is Young’s modulus, j is the creep compliance, and ν₁ and ν₂ are material constants. s_ij is the deviatoric stress given by the relation s_ij = σij − σ_kk δ_ij. The stress components shown above are all taken relative to the coordinate directions, x_i ʹ, of the neck.

Instead of directly solving the integral in Equation (13), Szyszkowski and Glockner have approximated the strain ė_ij ^c represented by the integral in Equation (13) by a nonlinear viscoelastic model of a generalized Kelvin body in series with a nonlinear potentiometer. This is shown in Figure 2. The equations for the Kelvin body are written as

(15)

(16)

(17)

is the effective viscous stress tensor defined above,

is the component of effective stress in the spring and

is the component of effective stress in the potentiometer. The second term involving the constant c on the lefthand side in Equation (16) is included to lend numerical stability to the constitutive equation at low stresses. Without it, the

curve becomes infinitely steep as

approaches zero, and this produces stability problems and computational difficulties in numerical schemes. Szyszkowski and Glockner avoided this problem by adding the second term to Equation (16). The value of c was generally between 0.1 and 0.3 and it had little effect on the stress behavior at intermediate to high stress levels. Accordingly, since this term is a computational convenience, a physical meaning cannot be attached to it. With this second term added to Equation (16) the strain e_ij ^c associated with the integral term was accurately represented. For a more detailed discussion the reader is referred to Reference Szyszkowski and GlocknerSzyszkowski and Glockner (1986, 1987).

Fig. 2. Spring potentiometer model for constitutive relation of ice.

Equations (15)–(17) are solved to find the strain rate ė_ij to replace the integral term in Equation (13), which can then be integrated to determine the strains in the necks.

As indicated earlier, the strain tensor components e_ij are defined relative to the x_i ʹ coordinates of the neck. The components E_ij relative to the global coordinates can be calculated by using the usual equation for transforming a second-order tensor from one coordinate system to another. This equation is

(18)

where Q is an orthogonal transformation matrix carrying the local coordinate system x_i ʹ centered at the neck to the global coordinate system x_i centered at the ice-particle center.

Interparticle Sliding

For an ice specimen in tension, failure is often defined as the instant when nucleation of microcracks is initiated. In compression, on the other hand, it is the stress required for the propagation of microcracks. The prediction of first crack occurrence (crack nucleation) under uniaxial compressive loading is based on the hypothesis that the crack nucleates due to lateral tensile strain resulting from the Poisson effect of elasticity and material incompressibility. The first crack in a neck is postulated to occur when the lateral tensile strain equals the strain for tensile fracture at the same instantaneous strain rate. If this limiting strain is reached, the material can still continue to sustain a compressive load but loses its ability to take any lateral tensile loads. The failure strain thus depends on the strain rate and a graphical relation between the two has been given by Reference Ting and SunderTing and Sunder (1985). For strain rates between 10⁻³ s⁻¹ and 10⁻⁶ s⁻¹ the tensile fracture strain varies from 3 × 10⁻³ to 6 × 10⁻³. In a neck there are both shear and compressive stresses. To determine when failure takes place, we calculate the principal tensile strain and, if this exceeds 5 × 10⁻³, we assume the neck has lost the ability to carry load (principal tensile stress) and therefore fractures. After the occurrence of fracture, deformation in the neck is primarily due to intergranular sliding. The tangential and normal velocities,

and ξ, and the displacement d across a neck are calculated using the equation

and

(19)

where

are respectively, the velocity and the component of traction in the direction normal to the grain surface at the neck, η is the tangential velocity, τ is the component of traction in the shearing direction and δ is the relative displacement between the grains. The expression 1/δ ^m appears because the sliding displacement increases as the particles begin to form new contacts with other particles, and this impedes further relative displacement. The constants c ₁, c ₂, n and m are determined from the experimental data with a single experiment. An unconfined compression test was used for this purpose. These constants are adjusted to provide a best least-squares fit to the data after the necks are calculated to begin fracturing in the materials determined by the above fracture criterion. This formulation was found to provide a very good agreement with the data. It should be noted that other material constants such as Young’s modulus, E, were also determined from the same test data.

If

are the components of slip velocity in the neck coordinate system, x_i , then the components l _ij of the velocity gradient, L, relative to the global coordinate system x_i are

(20)

where s is the unit vector directed normal to n and in the sliding direction and a is the radius of the ice particle. The symmetric part gives us the rate-of-deformation tensor which for small strain theory is the same as strain-rate tensor. The unsymmetrical part gives the spin tensor. The strain rate is

(21)

The above equation can be integrated with respect to time to give the strains. On multiplying both sides of Equation (21) (after time integration) with sin β dα dβ and integrating over α and β, we obtain

(22)

where E_ij on the righthand side is obtained from

Equation (21) or (18), depending on whether intergranular sliding or plastic deformation is taking place.

Comparison with Data

For the uniaxial compressive-stress case, the results for various values of axial stresses are shown in Figures 3–6 and compared to creep data from experiments run at Montana State University. Figure 7 illustrates a comparison of theory and data for uniaxial creep deformation of snow. For a grain radius of 0.508 mm, neck radius of 0.053 mm, coordination number equal to 2.5 and density ratio equal to 2.75, a Young’s modulus of 270.63 MPa was found to give good agreement with the data. For the same density rado, Mellor (1974) gave Young’s modulus of snow obtained from dynamical testing to range in value from 100 to 225 MPa.

Fig. 3. Strain vs time for a uniaxial compressive stress of 0.004 MPa. Test temperature = −10°C; a = 0.254mm; Α_T = 0.022 mm²; γ = 2.72.

Fig. 4. Strain vs time for a uniaxial compressive stress of 0.008 MPa. Test temperature = −10°C; a = 0.254 mm; A_T = 0.022 mm²; γ = 2.75.

Fig. 5. Strain vs time for a uniaxial compressive stress of 0.012 MPa. Test temperature = −10°C; a = 0.254 mm; A_T = 0.028 mm²; γ = 2.62.

Fig. 6. Strain vs time for a uniaxial compressive stress of 0.028 MPa. Test temperature = −10°C; a = 0.254 mm; A_T = 0.024 mm²; γ = 2.67.

Fig. 7. Strain vs time for a uniaxial compressive stress of 0.012 MPa. Test temperature = −10°C; a = 0.225mm; A_T = 0.020 mm²; γ = 2.67.

The ratio of lateral strains to axial strains, in most of the experimental data, varies from 0.15 to 0.2 for uniaxial compression tests and from 0.2 to 0.3 for uniaxial tension tests. The theoretical results give the value of ratio of lateral strain to axial strain as approximately 0.18 for compression and 0.23 for tension. An interesting feature of the lateral-to-axial strain ratio is that it does not change much even when the underlying deformation phenomenon is entirely different. At a compressive stress of 0.004 MPa when superplastic deformation of necks is the sole deforming mechanism, we have almost the same lateral-to-axial strain ratio as for a snow sample subjected to 0.028 MPa with interparticle sliding as the only deformation mechanism. This is because at low stresses the lateral deformation of the necks is into the pore space and does not contribute to the snow strains. For sliding particles there are no such lateral strains, and therefore the value of lateral-to-axial strain in this case is also almost the same. A better understanding about this ratio could be obtained if information were available about the plane of fracturing of ice under multi-axial loading.

The major advantage of this formulation is its ability to calculate strains under multi-axial loading without requiring any additional constants other than the ones used for uniaxial tests. This reduces the need to perform extensive experiments normally required to find additional constants arising in multi-axial constitutive laws.

For multi-axial loading, the constants used were those used in uniaxial compression tests. The results from these are plotted in Figures 8–10. Unfortunately no experimental data were available to compare with the theory. In Figure 8 the results from hydrostatic state of stress are shown. A 5.65% change in volume occurs over 30 h. This is in accordance with the compressible nature of snow. Figure 9 shows results from a test in which the stresses are very low and strains are due to superplastic deformation of necks rather than sliding of broken necks.

Fig. 8. Strain vs time for a hydrostatic state of stress. t₁₁ = t₃₃ =−0.008 MPa. Test temperature = −10°C; a = 0.254 mm; A_T = 0.022 mm²; γ = 2.72.

Fig. 9. Strain vs time for t₁₁ = −0.004 MPa; t₁₂ = 0.028 MPa. Test temperature = −10°C; a = 0.254 mm; AT = 0.022 mm²; γ = 2.72.

Fig. 10. Strain vs time for t₁₁ = −0.008 MPa, t₁₂ = t₂₁ = 0.008 MPa. e₃₃ differs significantly from e₂₂. Test temperature = −10°C; a = 0.254 mm; A_T = 0.022 mm²; γ = 2.72.

At very low stresses (0.004 MPa) the strain rates vary almost linearly with stress. Even at higher stresses (0.008–0.028 MPa) the non-linearity is not large as long as the coordination number, density, grain and neck size remain nearly constant.

In Figure 10, where t ₁₁ = −0.008 MPa and t ₁₂ = t ₂₁ = 0.008 MPa, the strains e ₂₂ and e ₃₃ are no longer equal. Although not shown here, when shear stresses t ₁₃ and t ₂₃ are applied along with a normal stress t ₃₃, the shear strain e ₁₂ is observed although no shear stress was applied in this direction. The value of this shear strain is one order of magnitude smaller than the strains in the directions in which stress was applied. In multi-axial tests when both shear stress and normal compressive stress are applied, the relation between deformation rate and stress is of the form

Since the dependence on the third term on the righthand side is rather small, as is apparent from relatively small value of e ₁₂ in the test mentioned above, the relation between deformation rate and stress can be approximated by

In the above equation C ₀ is a function of first and second principal invariant of stress tensor and second invariant of strain tensor. C ₁ is a function of second principal invariant of stress and second invariant of strain tensor. For both C ₀ and C ₁ the dependence on second invariant of strain tensor can be represented quite accurately by a second-order polynomial. However, no simple dependence could be found with respect to the other two variables.

Discussion

A theory has been proposed to explain the multi-axial deformation of snow. The theory is limited to small strains, and much more experimental work needs to be done before it can be extended to large strains. The major advantage of this formulation is its ability to calculate strains under multi-axial loading without requiring any additional constants other than those obtained from one uniaxial test. This is accomplished by means of defining the multi-axial deformation in terms of microstructural processes, i.e. the deformation taking place within the necks connecting the ice grains. This eliminates the necessity of performing extensive sets of experiments to evaluate highly empirical constitutive relations. This formulation is also capable of transcending those regions where different deformation mechanisms are dominant. Of course, significant disadvantages of this formulation are: (1) the need to describe accurately the various deformation mechanisms, (2) the need to define how these microstructural processes can be extrapolated to observed macroscopic results and (3) the complexity of these types of relations. The above work does represent one of the first attempts to rely heavily on processes at the microstructural level to define in detail the properties of snow with as little gross empiricism as possible. The work reported here is far from ideal but it does demonstrate that much can be accomplished with this approach.

Acknowledgement

The work reported here was funded by the Army Research Office under grant number DAAL 03-92-G-0310. The authors wish to express their appreciation for ARO’s support.

Notation

Ν Coordination number

n, ν ₃ Normal at the point of contact dn Differential solid angle

f_i (n) Contact force per contact

t_ij Components of stress on the snow sample

(global coordinate)

σ_i ⁽ⁿ⁾ Contact force per unit contact area

σ _i3 Components of stress tensor at the neck due to σ _i ⁽ⁿ⁾

e_ij Small strain tensor in the ice neck coordinate system (primed coordinate in Fig. 1)

S_ij Deviatoric stress

ν₁ Constant for delayed elasticity in ice

ν₂ Constant for viscous creep in ice

E_ij Strain in global coordinate system (unprimed coordinate in Fig. 1)

Q_ij Rotation tensor

σ In Equation (19), component of traction in normal direction

τ In Equation (19), component of traction in shearing direction

S_i Components of unit vector in direction of sliding of necks.

c₁, c₂, m, n Empirical constants

ξ Tangential component of sliding (tangent to the grain)

η Normal component of sliding