Introduction

Natural and mechanical compaction of snow can significantly affect the thermal and physical properties of snow, the estimation of the age of ice cores used to determine past climate conditions and cold-regions engineering operations. During compaction at low density, individual ice grains rearrange within the void space as the connecting bonds brake or undergo viscous or plastic flow due to stress concentration effects. At higher compaction rates, this process often proceeds through a series of stick/slip events as particles lock up and bonds fail with little deformation of ice grains. Eventually, ice grains become too closely packed for further movement and grain re-arrangement ceases. The critical density at which grain re-arrangement ends depends on grain shape, rate of deformation, temperature and moisture content and can vary, from about 200 kg m⁻³, for shock compaction to about 600 kg m⁻³, for natural compaction of snow on temperate glaciers (Reference Johnson, Solie, Brown and GaffneyJohnson and others, 1993). However, critical densities of about 550 kg m⁻³ are most commonly observed for natural and quasi-static mechanical snow compaction (Reference Anderson, Benson. and KingeryAnderson and Benson, 1963). Compaction of the snow beyond its critical density requires that ice grains deform, filling in pore space and increasing the number of particle contacts and their contact area. For natural snow compaction (low strain rate), power-law creep and diffusion (vapor, surface, grain boundary and lattice) are thought to dominate (Reference Maeno and Ebinuma.Maeno and Ebinuma, 1983; Reference Ebinuma and Maeno.Ebinuma and Maeno, 1987; Reference WilkinsonWilkinson, 1988). At higher strain rates (greater than about 0.001 s⁻¹), snow that is well-bonded and with densities generally greater than about 300 kg m ⁻³ may exhibit an initial elastic deformation to a yield stress and then follow along separate deformation paths which asymptotically approach the density of ice (Fig. 1). This rather complex compaction process results in a wide scatter of measured properties for different snow types with common density values even when other test conditions are similar (Reference Mellor.Mellor, 1975). It is widely recognized that much of the scatter in measured snow properties as a function of density is due to differences in snow structure (Reference KryKry, 1975; Reference Mellor.Mellor, 1975; Reference Voitkovsky, Bozhinsky, Golubev, Laptev, Zhi-gulsky and Slesarenko.Voitkovsky and others, 1975; Reference SalmSalm, 1982). Features of snow structure that may be important for compaction include the number and cross-sectional area of grain contacts, grain-shape and size distribution, and the geometry of connections. Efforts to measure these parameters and to relate them to snow-deformation behavior have proven to be difficult, making the development of representative constitutive models problematic.

Fig. 1. Uniaxial strain compaction of snow for compaction rates from 2.7 × 10⁻⁴ to 0.4 m s⁻¹. Data from Reference Abele and Gow.Abele and Gow (1975,Reference Abele and Gow.1976). Abele and Gow's data indicate little dependence on compaction rate for rates in this range.

A better understanding of the relationship between snow structure and its deformational behavior is needed before significant improvements to current models of snow deformation can be made. It is to this end that a study to investigate the feasibility of using numerical methods to help determine the relationship of internal snow structure to its compaction under uniaxial strain conditions was undertaken. It is an effort to understand better the role that structure plays in the deformation of snow at the unit-cell scale and to determine the important structural parameters that can guide future experimental and continuum-scale modeling efforts.

Numerical Modeling Methods for the Uniaxial Strain Compaction of Snow

A three-dimensional fully dynamic simulation of the uniaxial strain compaction of a collection of roughly spherical particles, assuming elastic-plastic constitutive behavior, was constructed using the PRONTO3D transient solid-dynamics program (Reference Taylor and Flanagan.Taylor and Flanagan, 1989). Recent enhancements to the PRONTO3D program enabling robust particle-contact algorithms, adhesion between particles and the modeling of large-scale shear deformation made the feasibility of successfully modeling of microscale processes in a realistic manner a possibility (Reference Heinstein, Attaway, Sweyle and MelloHeinstein and others, 1993; Reference Swegle and Attaway.Swegle and Attaway, 1995).

The simulation consists of a confinement vessel containing 113 randomly distributed 1 mm diameter elastic-plastic spherical particles with frictionless, adhesionless contact surfaces that are compacted using a piston moving at constant speed (Fig. 2). The diameter of the container equals seven particle diameters, which is sufficient to provide a unit-cell scale behavior and to eliminate any influence of the container walls for frictionless spherical particles in a frictionless cylinder. For particles with friction, adhesion or non-spherical geometries, a container diameter of ten or more particle diameters is needed to obtain unit-cell scale behavior. Each particle, the container and the piston are contact materials and have separately defined constitutive descriptions. A contact material may contact any other contact material without restriction so that particle compaction can proceed in a natural fashion without concern about where individual particles will travel. Once contact stresses exceed the yield stress, particles can deform to further reduce porosity, simulating the compaction of deformable particles. Limited computer capacity necessitated that each particle be constructed with a coarse mesh containing 48 three-dimensional quadrilateral finite elements (FE) or smoothed-particle hydrodynamic elements (SPH). Tests of adhesion were done using SPH while all other model runs used FE methods. SPH methods allow for contact definitions of either strain averaging (used to simulate adhesion) or strain isolation (transfer momentum only).

Fig. 2. Cut-away view of the simulated uniaxial strain-loading device including the confinement vessel, elastic plastic rough spherical particles and loading piston.

Results and Discussion

Model calculations were done for uniaxial strain compaction at constant compaction rates of 0.4, 1 and 20 m s⁻¹ and for two-particle systems to compare SPH adhesion and FE frictionless interactions. These were used to investigate the limit at which quasi-static loading breaks down and to examine how inertial stresses produce rate effects in a granular medium composed of a rate-independent matrix material.

Simulation results indicate that quasi-static conditions prevail for both the 0.4 and 1 m s ⁻¹ compaction rates. Inertial stresses that increase with compaction speed are evident at the beginning of compaction when the loading piston first contacts the particles. However, particle accelerations diminish quickly, so that the bulk stress at the top and bottom of the sample are the same and are unaffected by the compaction rate (Fig. 3). Compaction at 20 m s⁻¹ is dominated by inertial effects with a wave of deformation propagating from the piston face into the aggregate of particles. Particle deformation begins next to the loading piston, because of the large inertial stresses so that little re-arrangement of particles occurs before they become locked in place. Hence, a higher stress is required for the 20 m s⁻¹ compaction than for quasi-static compaction to achieve the same compaction relative density. The differences in microstructural evolution for a compaction rate of 20 m s⁻¹ as compared to 0.4 or 1 m s ⁻¹ are responsible for the strain-rate dependence even though the material from which the particles are composed is rate-independent. Simple two-particle interactions using SPH methods demonstrate that contact strain averaging can be used to model the influence of adhesion. This opens the possibility of examining the deformation of sintered particles in a future study.

Fig. 3. Calculated normalized axial piston stress and normalized containment vessel side-wall stress vs relative density from the uniaxial simulation. The maximum piston stress was used as the normalization factor. The relative densities of particle packings are indicated, for reference, by vertical arrows for cubic-packed loose-packed spheres (CLPS), random loose-packed spheres (RLPS), close-packed spheres (CPS) and cubic loose-packed oriented disks (CLPD).

The quasi-static compaction evolved in four stages identified by slope changes in the axial stress/relative density curve for a 0.4 m s⁻¹ compaction rate shown in Figure 3. The microstructural state of the aggregate at the end of each compaction state is shown in Figure 4. The main features are that, except during elastic deformation, the number of particle contacts increase and, for plastic deformation, the area of contacts increases throughout the compaction process. During the first compaction stage, particles re-arrange at low stress. Stresses appear to be centered around zero until the very end of this stage when particles start to lock up. Second-stage compaction is elastic and occurs after the particles have reached a stable-packing arrangement at a relative density of P _r = 0.55, slightly lower than that for the random loose packing of spheres. Third-stage compaction progresses through the plastic yield of the particles and begins at the particle/loading piston contacts, where stress concentrations are largest, then continues at particle contacts. Deformation is greatest along the axis of principal stress which produces a flattening of particles. Compaction continues in this manner until ρ _r = 0.75 which is between the relative density of close-packed spheres (ρ _r = 0.74) and cubic loose-packed oriented disks or columns (ρ _r = 0.785). Beyond this density, further flattening of particles to fill available void space becomes more difficult because of lateral constraint, causing radial stresses to increase so that deformation becomes more uniform and particles begin to indent each other. The stable packing of the deformed particles and increasing radial stress produce higher mean stresses in the particles which then require a larger axial stress to achieve the same deviatoric stress condition needed to continue the compaction process as before the packing arrangement was reached. This is conveyed as the increased slope of stage IV compaction (Fig. 3).

Fig. 4. Selected isolated particles taken from the center-line cross-section through the aggregate of particles at zero stress (a), critical density 1 (b), critical density 2 (c) and critical density 3 (d).

The stable-packing arrangements of undeformed or deformed particles that separate compaction stages define critical densities and constrain further particle movement or deformation. As a consequence, the axial stress to density ratio must increase after a critical density has been reached for compaction to continue at the same rate. This implies that an additional critical density may occur when the relative density of close-packed oriented disks is reached (0.91). Evidence for such a packing in the densification of polar ice is discussed in the next section.

Critical Densities in Pressure-Density Profiles of Polar Snow

Three critical densities have been identified in pressure density profiles of polar snow (Fig. 5) (Reference Maeno and Ebinuma.Maeno and Ebinuma, 1983). The first critical density occurs at a relative density of ρ _r = 0.6 (550 kg m⁻³) (Reference Anderson, Benson. and KingeryAnderson and Benson 1963), which is identical to the relative density of random loose-packed spheres. The second critical density ρ = 0.8 (730 kg m ⁻³) was described by Reference Maeno and Ebinuma.Maeno and Ebinuma (1983) as being the density of optimum bonding and packing of constituent snow particles and exhibits a cylindrical pore geometry at the intersection of several grain boundaries. Gow (in Reference Ebinuma and Maeno.Ebinuma and Maeno, 1987) has indicated that at some locations the change in pressure-density slope for densification of polar snow may be attributed to shear deformation as the principal cause rather than optimal particle packing. The third critical density occurs at ρ _r, = 0.89-0.92 (820-840 kg m ⁻³) and is associated with a change in pore geometry from cylindrical to spherical.

Fig. 5. Pressure vs density for natural snow compaction at Byrd Station and Little America V, Antarctica. Data taken from Reference GowGow (1968). The relative densities of particle packings are indicated, for reference, by vertical arrows for random loose-packed spheres (RLPS), cubic loose-packed oriented disks (CLPD) and close-packed oriented disks (CPD).

Critical densities are readily apparent in the pressure density profile for Byrd Station (Fig. 5), where critical densities 1 and 2 for the Byrd Station profile correspond to the critical densities 1 and 3 of the simulation. The pressure density curve for Little America V shows a different behavior than for Bvrd Station or the simulation. Between critical densities 2 and 3 the pressure-density ratio decreases. This is consistent with the action of shear-enhanced compaction processes and indicates that both uniform compaction and shear-enhanced compaction may occur in polar snow. The shear-enhanced compaction appears to be noticeable only between the second and third critical densities (Fig. 5). While possible, it is unlikely that shear stresses are coincidentally confined only to the depths that lie between the second and third critical densities. A more feasible explanation is that the geometrical arrangement of the ice grains prior to the second critical density is such that shear does not significantly increase the ease at which deforming ice fills available voids. After the second critical density, however, the constraining influence of the optimal snow-grain packing causes the mean stress to increase more rapidly for a given applied axial stress so that deformation of ice grains due to deviatoric stress is reduced significantly, so that continued shearing may enhance material flow into available pore space. After the third critical density is reached, the effects of shear-enhanced compaction may decrease as the pressure of air trapped in the pores increases sufficiently to restrict additional creep deformation.

Conclusions

A dynamic finite-element computer code has been used to simulate successfully the micromechanical quasi-static compaction of roughly spherical elastic-plastic particles. Simulation results indicate that stable-packing arrangements of undeformed and deformed particles that constrain further movement of particles or flow of material occur at distinct densities (critical densities). The additional constraint caused by the stable packing of particles causes the mean stress to increase more rapidly Ibr a given axial deformation which requires that the stress-density ratio increase (in the absence of shear) to maintain a constant compaction rate. In the simulation, critical densities occurred at the end of particle re-arrangement and during plastic flow. The results of the simulation, excluding clastic response, are similar in appearance to observed pressure-density profiles for polar snow (in the absence of shear). Shear stresses appear to enhance snow compaction between critical densities 1 and 2 but do not appear to affect significantly the other compaction stages. This study has demonstrated the usefulness of using numerical simulations to understand better the role of microstructure on continuum-scale deformation.