Introduction

Dry flowing-snow avalanches are made up of particles of snow and ice. These particles may be individual snow or ice crystals or larger clumps of sintered snow. A slab avalanche, for example, starts out as large slabs and big chunks of snow, and as the avalanche descends the slope, the slabs and chunks get broken into smaller pieces. Particularly near the bottom of the avalanche, where the shear and normal stresses are greatest, the snow is ground into ever smaller pieces which may approach the size of the individual ice grains that make up the snow fabric. Also augmenting the mechanism of grinding at the bottom of an avalanche is a particle-sorting process whereby, in a gravity field, smaller particles fall to the bottom and larger particles rise when a granular material is sheared (Reference Haff and WernerHaff and Werner, 1986; Reference SavageSavage and Lun, 1988). On the other hand, near the top surface of an avalanche, where stresses and deformation may be small, the size of the snow chunks can remain large.

Evidence from the field (Reference Dent and LangDent and Lang, 1982; personal communication from H. Gubler) indicates that in a flowing avalanche moving over a relatively flat slope the majority of deformation takes place near the base of the avalanche, resulting in large velocity gradients at the basal surface. Above this active region of shearing, deformation may occur, but its magnitude is much reduced from that which takes place below. The dynamic characteristics of this bottom shear region, made up mainly of small ground-up ice grains, in large part determines the dynamic friction force acting on the bottom of an avalanche.

The dynamic properties of this active shear region at the base of an avalanche are investigated using a very simple model. The model consists of a two-dimensional array of round uniform spheres undergoing rapid shear ( Fig. 1). The vertical depth of the array corresponds to the depth of the high shear-rate region at the base of the avalanche in the slope-normal direction, and the horizontal width of the array corresponds to the slope-parallel length of the avalanche. Three-dimensional effects, such as path constrictions and lateral slope variations, are necessarily neglected by the two-dimensional nature of this model. Furthermore, only steady,uniform flow conditions are considered. The weight of the avalanche mass above the high shear region is applied to the shear layer as a normal stress N, pushing down on the top of the particle array. This normal stress is proportional to the total depth of the avalanche. The motion of the avalanche mass is imparted to the array as a horizontal velocity, U, applied to the top layer of particles. The bottom layer of particles in the model is fixed in order to simulate the stationary snow surface over which the avalanche moves. The array of particles, driven by the moving top layer, and constrained by the applied normal force, is allowed to expand or contract in the vertical direction in response to the internal forces generated by the shearing motion of the array and the applied normal force.

Fig.1 . Two-dimensional shear-layer model shown here with four layers of particles between stationary bottom layer and top layer moving horizontally with velocity U. Array is compressed by normal stress, N, which is due to weight of overburden.

This simple two-dimensional array of particles, modeling the shear layer at the base of an avalanche, is used to determine how the shear force is transmitted from the top to the bottom of the shear layer. This shear force is the model’s representation of the effective dynamic friction force acting at the base of an avalanche. In particular, the ratio of this shear stress to the applied normal stress is calculated as a function of the velocity of the top layer of particles, the normal stress applied to that top layer, and the material properties of the spheres that make up the shear layer. This ratio is equivalent to an effective coefficient of dynamic friction for the shear layer. The behavior of the dynamic friction acting at the base of this very simple avalanche model may be used to extrapolate to the behavior of real snow and real avalanches. That is, the form of the solutions found here may help established how, as a function of the depth, speed, and the material properties of the snow particles, the shear layer of snow at the bottom of a dry-snow avalanche behaves.

Computer Model

The primary means of analyzing the shear-layer model described above is the construction of a two-dimensional computer simulation of the interacting particles that make up the shear layer. The simulation is similar to computer models that have been constructed by Reference Savage and LunWalton (1986), Reference Campbell and BrennenCampbell and Brennen (1985), and others, and is a direct descendent of the model described by Reference DentDent (1986). The simulation operates by integrating the equations of motion of each particle in the shear array over a succession of small time increments, in order to compute the velocity and position of each particle as a function of time, while the particles collide and interact with one another. The particles are modeled as uniform spheres that are rough, deformable and inelastic. A normal contact force and a tangential friction force are employed whenever two particles interact. Interactions are considered to occur whenever the boundaries of two particles overlap. The normal contact force, F_n , is determined linearly (with spring constant k) from the amount of overlap, δ between the two colliding spheres, and this force is then applied to the center of each sphere:

F_n = kδ.

Energy dissipation in the collisions is accomplished by using two different spring constants: k₁ as the two spheres approach each other, and k ₂ as the two spheres rebound. It can be shown that the square root of the ratio of these two spring constants corresponds to the coefficient of restitution for the collision:

During a collision, a friction force acts tangentially to the surface at the point of contact of each impacting sphere. This interparticle Coulomb friction force, F_t , is taken to be proportional to the normal contact force found from the overlap of the two spheres and is always applied in a direction to oppose the relative tangential motion between the contact points.

where μ is the interparticle Coulomb friction coefficient. Finally, although these forces are not considered in this analysis, a gravitational-body force and a cohesive-contact force may also be specified for each of the spheres.

Using these mechanisms, all the particles in the shear layer can interact with one another over discrete time increments. The forces on each sphere at the beginning of each calculational time step are used in Newton’s laws of motion to determine the two components of linear acceleration and one component of angular acceleration for each particle. The accelerations are then integrated over the time step to produce new positions and velocities for each particle. From the new positions of the particles, overlaps are recalculated, and new forces are determined to start a new time step. This cycle of calculations is then repeated. In this way the forces acting on the spheres and the motion of all the particles in the shear layer may be computed as a function of time.

The simulation is typically taken to be only several particle layers deep. For most of the results reported here, the shearing is started with an array of particles that contains four layers of closely packed spheres in addition to the fixed stationary bottom layer and the rigid velocity-driven top layer. The spacing and size of the spheres in the two boundary layers may be varied to simulate different roughness conditions at the boundaries, but the results reported here are for spheres of the same size as the particles in the array, with no space between adjacent boundary particles.

Numerically, only a finite number of particles can be manipulated. For this reason, and also to be able to produce quasi-steady uniform motion, the number of particles in each layer of the array is limited by specifying periodic boundary conditions in the horizontal direction. These boundary conditions, at the left and right ends of the computational domain, are constructed so that when a particle leaves the boundary through the one side, it reappears at the other boundary moving with the same velocity and at the same height. This boundary is denoted in Figure 1 by the two dotted vertical lines at the left and right sides of the particle array. In Figure 1, particles denoted a and b are images of the same particle as it passes through the periodic boundary. The same is true for the single particle shown as c and d. This boundary condition effectively makes the shear layer infinite in the direction of the applied motion, and it also imposes a uniformity on the motion of the array. Motion repeats itself horizontally with a period equal to the width of the periodic boundary condition. In effect, it is as if the 2-D computational plane curves back around on itself to form a cylindrical domain, with the circumference of the cylinder equal to the length of the periodic boundary. Periodic boundaries also permit uniform flow to be achieved without having to deal with leading or trailing edge effects.

The computer simulation begins with the particles in the system arranged in a close-packed array and with each particle possessing a random velocity. The top boundary layer is set in motion horizontally at a constant specified speed, and all the particles, including those in the top and bottom boundary layers, are allowed to interact through normal contact and tangential Coulomb friction forces. The top boundary layer remains rigid while moving horizontally, and is allowed to respond vertically as a unit to the applied normal force and the forces exerted on it from collisions with the particles in the array. The motion of the system of particles evolves as each particle bounces against other particles and the two boundary layers. The particles acquire a fluctuational kinetic energy and assume some mean spatial distribution which depends upon the top boundary speed, applied normal force, and particle properties. The mean particle spacing determines the total dispersion in the array, which is measured by the average distance between the top and bottom boundary layer of particles.

The forces exerted on the top and the bottom boundary layer due to collisions with the particles in the shear array are computed and averaged over time. The simulation is allowed to proceed until the fluctuations in the calculated time averages become small. From the average vertical force on each boundary, the vertical normal stress is computed and compared to the applied normal stress, N, to ensure that steady flow has been achieved and transient effects have averaged out. The average horizontal force on each boundary is used to calculate the average shear stress transmitted by the array of particles. The shear stress is then divided by the normal stress to determine the effective dynamic friction coefficient. This dynamic friction coefficient is found to depend upon the speed of the top boundary layer; the normal stress applied to the top boundary layer; and the properties of the particles in the array, such as their mass, coefficient of restitution, and interparticle Coulomb friction coefficient. As an aid in visualizing the results, a graphics package was written to display the animated motion of the particles.

Results

A dimensional analysis can be used to simplify the relationships between the various factors that determine the dynamic friction coefficient. If it is assumed that the dynamic friction coefficient may be a function of the shear speed, U, the applied normal stress, N, the particle coefficient of restitution, e, the coefficient of interparticle Coulomb friction, μ, the diameter of the particles, d, and the mass density of the particles, ρ then dimensional arguments require that the dynamic friction be a function of e and μ and the dimensionless parameter [Inline 1] This result indicates that equally varying U or the square root of the particle density, √ρ, or inversely varying the square root of the normal stress, 1/√N, should all have identical effects on the dynamic friction coefficient.

Figure 2 shows the simulation results of finding the dynamic friction coefficient as a function of the dimensionless number The xs are results for varying U while holding ρ, Ν and d, the particle diameter, constant. The circles are results for varying Ν while holding ρ, U and d constant. The triangles are for U and ρ held constant while varying Ν and d. Each data point in Figure 2 represents about 8 h of computer time on an IBM RS6000 model 530h. Typically, each simulation uses about 30 × 10⁶ time steps of 1 μs, yielding 30 s of model time, which is required in order for the fluctuations in the forces to average out and for the results to converge to three significant figures. The material parameters for the spheres used in the simulations that produced Figure 2 were chosen to conform to particles in a physical experiment that is under construction. The spheres were typically given a diameter of 8 mm, a density of 1231 kg m⁻³, a coefficient of restitution of 0.7, and an interparticle friction coefficient of 0.7. The results from Figure 2 indicate that the dynamic friction increases with increasing shear speed, U, and decreases as the normal stress, N, is increased.

Fig.2 . Numerical-modeling results for the shear-layer model showing the dynamic friction coefficient as a function of the nondimensional shear speed, using an interparticle Coulomb friction coefficient of 0.7.

Shown in Figure 3 are simulation results for the dynamic friction coefficient using two other values of the interparticle Coulomb friction coefficient, μ = 0.0 and μ = 0.2. Even for the case of zero interparticle friction, there is a substantial dynamic friction force over much of the velocity range. Since there is no interparticle Coulomb friction, this must be due solely to momentum transfer by the normal collision process. In this case the dynamic friction comes from the horizontal component of the normal collision force between particles.

Fig.3 . Numerical-modeling results for the shear-layer model showing the dynamic friction coefficient as a function of the nondimensional shear speed, using an interparticle Coulomb friction coefficient of 0.0 and 0.2.

In all the simulations, it was found that the dispersion of the array of particles, measured by the average distance between the top and bottom boundary layers, was an increasing function as was the average kinetic energy acquired by the particles. For all three interparticle Coulomb friction conditions shown in Figures 2 and 3, at the highest shear rates the dynamic friction coefficient begins to approach the value of the interparticle friction coefficient. This is particularly evident in Figure 3 for the zero particle friction case, where at high shear rates the dynamic friction approaches zero.

Figure 4 shows how the dynamic friction coefficient varies with the number of layers of spheres in the shear flow. The simulation was run first with one layer of particles between the top and bottom boundaries, then in each subsequent run the number of layers was increased by one. The data plotted is in each case for a top boundary speed of 20 m s⁻¹ and an applied normal stress of 2000 Pa. As can be seen from the figure, the dynamic friction is essentially independent of the number of layers of particles in the shear flow.

Finally, it must be mentioned that changes in the length of the periodic boundary condition in the simulations had a slight affect on the quantitative results that are reported here. As a result, for all the simulations, a constant periodic boundary length was used. The shapes of the curves reported in Figures 2–4, and the qualitative results discussed below, should therefore not be influenced by the periodic boundary length.

Fig.4 . Numerical-modeling results for the shear-layer model showing the dynamic friction coefficient as a function of the number of layers of particles between the top and bottom boundary layers. The same flow conditions and material properties are used in each calculation.

Discussion

Understanding the physical processes that are producing these dynamic-friction results begins with a consideration of the data plotted in Figure 4. The dynamic friction for a single layer of particles is the same as the dynamic friction for many layers of particles between the boundaries. Also, in observing the animated motion of the shear layer, it is seen that the shearing motion at any instant in time is taking place primarily between only two of the shear layers. A shear plane or slip plane forms momentarily between two of the layers, and the majority of the shearing in the array at this instant takes place here. It is found that the location of this slip plane moves sporadically and randomly throughout the shear layer. The movement of the slip plane may be enhanced somewhat by the absence of a gravity force acting on the particles in the simulation. The weight of the particles would produce a normal stress gradient that increases the normal stress in the shear flow from the top to the bottom of the array. Since for steady motion the shear stress must remain constant over the depth of the array, the ratio of shear to normal stress would decrease toward the bottom of the array. In the absence of other factors, this lower dynamic friction at the bottom would tend to concentrate the shear plane at the bottom of the array.

The formation of the shear plane seen here may in part be attributed to the imposed periodic boundary condition. By having only a small number of particles across the shear flow (typically 12 in these calculations) the slip plane is formed by only that number of particles and is relatively easy to develop. The periodic boundary condition then effectively replicates that slip plane so that it becomes infinite. However, slip-plane formation has also been observed in experimental shear devices. One device which is being built by this author, designed to produce two-dimensional shear flow similar to that modeled here, and another annular shear cell built by Reference Hanes and InmanHanes and Inman (1985) to produce three-dimensional shearing of glass beads have both observed the formation of slip planes.

The layered structure and the formation of slip planes reported here contrasts with observations of shallow-chute flows of granular materials. For shallow flows, with no effective overburden weight, it is found that the shear motion is distributed throughout the granular flow regime, producing a velocity gradient over the entire depth of the flow (Savage, 1979; Reference DrakeDrake, 1990). In addition, since there is no overburden material to constrain the flow, the dispersion of particles in shallow flows is much greater than observed for the flow conditions modeled here, and the layered structure is not seen. Essentially, for shallow flows the layered structure, and as a result slip-plane formation, is lost. Only when the flow becomes sufficiently deep, or sufficiently slow, does the structured motion begin to be seen.

The formation of a slip plane allows the behavior of the dynamic friction coefficient at high shear rates to be interpreted by considering the impact between the two layers of particles adjacent to the slip plane. Consider the impact angle, θ, for two of these particles as shown in Figure 5. θ is defined as the angle between the line connecting the centers of the impacting particles and the horizontal. First, consider the effects of the normal impact force between the colliding particles. For a given impact velocity, as θ is increased toward 90°, the horizontal momentum transferred in the collision is decreased while the vertical momentum transferred is increased. Momentum transfer across the slip plane results in a force being transmitted across the slip plane. The ratio of horizontal to vertical momentum determines the ratio of the shear to normal force transmitted across the shear plane for this collision, which in turn determines the dynamic friction coefficient. Thus, for each impact, the dynamic friction due to the normal impact is decreased as θ increases toward 90°.

Fig.5 . The definition of the impact angle, θ, between colliding particles in the two layers of particles adjacent to the shear plane.

The collision rate between particles in the two slip layers is proportional to the slip velocity between the two layers. The greater the slip velocity, the less time it takes for a top-layer particle to traverse the bottom layer and encounter the next bottom-layer particle. As the time between consecutive collisions becomes smaller, the top layer has less time to respond vertically to the applied normal force, and the distance that the top layer of particles is able to penetrate into the space between adjacent bottom particles is decreased. As the shear velocity is increased then, the collision angle between particles in the two layers adjacent to the slip plane increases toward 90°, and the dynamic friction coefficient due to normal impacts decreases. At very high speeds only glancing collisions occur at impact angles near 90°. The normal impact contribution to the dynamic friction coefficient then becomes very small. For the case of zero interparticle friction, where normal impacts are the only mechanism to transfer momentum, the dynamic friction must tend toward zero at high shear speeds.

When an interparticle Coulomb friction force is introduced, a collision angle near 90° results in the Coulomb friction force becoming the primary contribution to the horizontal shear force transmitted across the slip plane, and therefore the major factor contributing to the dynamic friction coefficient. At high shear speeds the dynamic friction coefficient therefore approaches the value of the interparticle friction coefficient. At lower speeds, as the collision angle decreases, the dynamic friction is a composite effect of the normal impact force and the tangential Coulomb friction force.

Recall from the dimensional analysis that increasing the normal stress applied to the top of the shear layer has the same effect as decreasing the shear speed. Reverse arguments can therefore be made for the effect of the normal stress on the dynamic friction coefficient. A large overburden load, which produces a large normal force, compresses the entire shear layer. The particles in the two layers adjacent to the slip plane are forced closer together. Between consecutive collisions a particle from the top layer is forced deeper into the space between adjacent bottom particles. This decreases the collision angle and has the same effect that occurs when the shear speed is decreased.

Conclusions

The computer study reported here is a highly idealized simulation of one aspect of the motion of dry flowing-snow avalanches. Many mechanisms that are important to avalanche dynamics have been ignored to focus on a simplified picture of the mechanics of the shear layer at the base of an avalanche. Snow entrainment and air drag are just two of the effects that have not been considered here. In addition, the simple shear-layer model employed here ignores many significant features of snow and avalanche motion: (1) snow is not made up of round uniform particles with uniform properties; (2) the surface over which an avalanche flows is not a smooth plane; (3) avalanche motion is neither steady nor uniform; and (4) in many avalanches three-dimensional effects are important. Other important dynamic mechanisms, such as cohesion effects, might also have to be eventually considered. The importance of each of these effects is difficult to evaluate quantitatively and under certain circumstances any single mechanism could be dominant. But as a starting point for understanding the behavior of a complicated phenomenon like an avalanche, there is utility in simplicity. If an understanding of simple models can be accomplished, then other mechanisms may be studied and added to produce a model that more closely resembles the real phenomenon.

What has been presented here are results for a very simple flow model representing the shear at the base of an avalanche. What has been found is that for this simple model, for reasonable values of the interparticle friction coefficient (0.7–1.0 for rough, non-uniform grains) and the coefficient of restitution (0.5–0.8 for ice), the dynamic friction coefficient in the shear layer increases with the speed of the avalanche. Also, since increasing the depth of the avalanche increases the normal stress acting on the shear layer, the dynamic friction coefficient decreases as the depth of the avalanche increases. Disregarding other effects, this implies that large avalanches should move faster and run out farther than small avalanches.

Additionally, the results presented here provide an explanation of why, even in a cohesionless granular material, the dynamic shear layer should be located at the bottom of the flow and not elsewhere in the flow, or not equally distributed throughout the flow. The normal stress is greatest at the bottom of a deep flow of granular material due to the overburden weight of the material above. The dynamic friction coefficient has been shown to decrease when the normal stress is increased ( Fig. 2), and as a result, the dynamic friction coefficient for steady uniform flow is a minimum at the bottom of the shear flow. So slip or shear occurs preferentially at the base of the flow.