Introduction

A two-dimensional linear viscous model has been used to simulate the mechanics of flowing snow (Dent and Lang 1980), It was found that for small test slides of less than 20! m s⁻¹ the model, with several modifications, provided reasonable simulation. One modification to the linear viscous model that was found necessary was the inclusion of a friction boundary layer between the flowing snow and the stationary snow surface. This friction condition introduced a second parameter into the model. By adjusting the friction coefficient and the viscosity coefficient, the model was able to simulate the snow tests.

A major failure of the linear viscous model for flowing snow was its response to low stresses. In the modeling, it was necessary to halt the computations at the point where the leading edge of the flow fell below an arbitrarily small velocity. The continuation of the calculations beyond that point would have eventually allowed the fluid to deform until the depth of the material was reduced to zero and the horizontal dimension had become infinite. This motion exemplifies one obvious difference between flowing snow and a linear viscous fluid. Snow is seen to come to rest with a finite depth where the stress is nonzero. This indicates a threshold stress state in snow that must be overcome in order for deformation to take place. This property is due in part to the cohesion of the individual snow particles, but is mainly a result of the granular nature of the material.

The simplest continuum mechanical model to exhibit this locking property is called a Bingham material. The constitutive equation for a Bingham material is made up of two parts. First, if the stress intensity is below a threshold value T₀, no deformation takes place. Second, if the stress intensity is above this value, deformation takes place, and is proportional to the amount that the stress level exceeds τ 0. In two dimensions, for positive shear stresses τ, the constitutive relation can be expressed by

and

where u and v are the components of the velocity in the x and y directions and u is a constant parameter similar to the viscosity coefficient for a linear viscous fluid. This is a simplified version of the general Bingham equation that can be found in Reference MalvernMalvern (1969), for example. In that form the equations must be expressed in a manner that is frame-invariant. This requires the yield stress t0 to be expressed as a function of the stress invariants. The detailed development of this and subsequent equations in a general two-dimensional form can be found in Dent (unpublished). Well-known materials of the Bingham type include paints, greases, concrete, and toothpaste.

In addition to modeling the locking property of flowing snow, a Bingham model contains the necessary stress-deformation components to model the boundary layer that was treated as a friction force In the linear viscous model. In part, this is due to the additional parameter T0 involved in the Bingham equations. But also the very nature of the Bingham model, being physically more accurate, allows a more realistic representation of the motion, and, as will be seen by the results, provides a very good fit to the data.

Biviscous Model

The implementation of the Bingham model proved to be a difficult task. Of primary importance to this model is the location of the yield surface that separates the two flow regimes. On one side of this surface, the material is locked and behaves as a rigid body. On the other side, the constitutive equation, when substituted into Cauchy's equations, provides the Mavier-Stokes equations. In both cases, existing methods allow the governing equations to be integrated. Unfortunately, the calculation of the location where τ = τ₀ is not a simple procedure, since the constitutive equation does not define the stress when it is less than t0. For this reason, an alternate formulation to perform essentially the same task was considered.

This new approach allowed small deformations to take place according to a linear viscous flow law in the locked portion of the flow. The viscosity used in this region is taken so high that the resulting deformation can be neglected relative to deformations outside the region. The small deformations and linear viscous flow law allow calculation of stress values from the constitutive equation. Location of the yield surface, T = T₀, is then easy to find. Outside the region of small deformations, as in the pure Dinghan model, the flow is still 1 i near-viscous, but with a different viscosity. This two-viscosity system was dubbed the biviscous modified Gingham model. Figure 1 shows a one-dimensional characterization of this flow law along with the corresponding pure Bingham relationship. The parameter v, analogous to the kinematic viscosity, has been introduced in this figure and is equal to μ/ρ, where p is the mass density of the material. The mathematical representation of the biviscous model is

and

Fig. 1. One-dimensional form of the Bingham and biviscous constitutive laws.

where μ and μ are the viscosities in the two flow regions. The term T₀/μ in the second equation is a small correction to account for the fact that the velocity gradient is not zero at the point where the stress intensity becomes equal to T₀ (see Fig.l). Again these equations are a simplification of the general two-dimensional constitutive equation (Dent unpublished). In the general form, the yield surface represented by T₀ must be found according to a frame-invariant yield condition. The simplest method, and the method used in this model, is the equivalent of the Treska yield criterion found in solid mechanics (Reference MalvernMalvern 1969).

The general two-dimensional biviscous constitutive equation was implemented by using the framework of a numerical code to solve the two-dimensional incompressible Navier-Stokes equations (Amsden and Harlow 1970). The flow was assumed to be incompressible, which simplifies both the analysis and the resulting equations. This assumption is dubious at best, but little data are available to check it. The resulting code utilizes a marker and cell method which finite differences the governing equations. Stresses are calculated at the cell nodes using the constitutive equation and and the kinematic flow field. These stresses are then used in a finite-difference approximation to the momentum balance equation to determine the advanced time flow field. The exact implementation of this procedure is again detailed in the thesis by Dent (unpuhli shed).

Biviscous Modeling Results

The biviscous model was used to simulate the tests on flowing snow described in Dent and Lang (1980, 1982). These tests decelerated 2.2 m³ of snow from 18 m s⁻¹ to rest, on a level runout of packed snow. Data on the position of the leading edge of the snow, velocity versus depth, and final distribution of debris were collected. Also, qualitative information on the mechanics of the flawing snow was gathered.

Since the flow entered the runout area from an essentially friction-free polyethylene surface, it was allowed that the initial configuration would be a mass of material moving at constant speed on a horizontal friction-free surface. The initial velocity of this material was taken to be 17 m s⁻¹, which was derived from the initial slope of the curve relating position to tine. The spatial dimensions of this material were determined from film footage taken of the test and are illustrated in Figure 2.

Fig. 2. Input flow confiquration for computer simulation of snow flow.

The numerical modeling commenced with the flow of the material off the frictionless surface onto a surface employing a no-slip boundary condition. The computational grid consisted of an area 23 n long and 50 cm high. The horizontal dimension was divided into 140 cells, each 0.20 in long, and the vertical dimension divided into 10 cells each 0.05 m high. This proved to be about the minimum grid size that was economically feasible. A smaller cell size was tried for an abbreviated run and the results showed little overall variance from results of a similar test on the 0.20 × 0.05 m grid. A larger-celled grid was, however, deemed inappropriate since the boundary layer at the bottom of the flow was of the order of 5 cm. Cells with vertical dimensions larger than 5 cm would be unable to resolve this layer. The horizontal dimension was then chosen to provide reasonable resolution in that direction and to maintain an aspect ratio between the cell dimensions of no more than 5 to 1.

The three-program modeling parameters, T₀, ν = μ/ρ, and ν’ = μ’/ρ, were then adjusted so that the computed flow conformed to the observed motion of the test. It was found that the parameters T₀ and ν were principally responsible for the motion of the leading edqe and the distance of total runout. However, many different combinations of T₀ and ν produced the same runout. Flow velocities were not large enough to provide definite distinctions between combinations of these parameters. Figure 3 shows several one-dimensional equivalent constitutive relations involving combinations of T₀ and ν that gave good results for leading-edge motion.

Fig. 3. Constitutive relations that modeled the motion of the leading edge of snow.

The velocity profile measured in the window test provided another criterion to be satisfied by the numerical simulation. It was found that the comnutergenerated velocity profile was also principally a function of the two parameters T₀ and ν. As T₀ was increased, and v decreased, to maintain the same leading-edge characteristics, the velocity profile became sharper, with larger gradients near the surface and smaller gradients above. Conversely, combinations of snail T₀ and large v produced gradients more closely resembling the parabolic shape expected for pure viscous fluids. Matching the shape of the velocity gradient provided the necessary information to define the two parameters -c0 and v uniquely. These two values were found to be: for T₀, expressed in units of stress per unit density, 2.20 m² s⁻²; and for v, the kinematic viscosity, 0.002 s² s⁻¹. It was also noted that these values provided the best comparisons of leading edge versus time with the experimental snow test. This comparison is illustrated in Figure 4. Figure 5 shows examples of the velocity gradient calculated by the computer model, corresponding to the location of the data acquired in the snow tests. Also shown is the profile found from the snow tests, which is plotted on a velocity scale twice that of the other plots because the velocities measured fron behind the window were about half those measured for the motion at the center of the flow. It is believed that this is due primarily to the boundary drag exerted on the edge of the flowing snow. The velocities of the flow measured from the window were about 7.0 m s⁻¹, at the leading edge. Meanwhile, at the center of the flow, the leading edge was found to be moving at nearly 16 m s⁻¹.

Fig. 4. Position of the leading edge versus time: comparison between experiment and computer model.

Fig. 5. Velocity profile comparison between snow test and computer model. Calculations for various combinations of model parameters.

The magnitude of the third parameter v' was found to have very little effect on the motion of the leading edge. The velocity profile, however, was affected by this parameter, though small adjustments of t_o and ν could be made to compensate. It was also found that v' had a pronounced effect, with T₀ and ν, on the final distribution of debris. As ν’ increased, deformation in the upper regions of the flow decreased. This resulted in less total deformation of the initial flow configuration. A value of v' equal to 0.10 m² s⁻¹, combined with the previously specified values of x0 and u, provided the best comparison of final depth profiles of the debris. This result is is plotted in Figure 6, as well as the results of a simulation with ν = 0.20 m² s⁻¹. Figure 7 shows a full time series of particle plots for this simulation. In these plots, the friction-free surface extends from the left boundary to the 8.00 m mark. From there onward the surface is no-slip. The vertical dimension (labeled depth) is plotted on a scale exaggerated by a factor of 4 over the horizontal scale.

Fig. 6. Final depth profile of debris; comparison between experiment and computer model.

Fig. 7. Time sequence particle plot of computer simulation of snow test.

As can be seen from examining Figures 4, 5 and 6, the modeling results, with T_O/ρ = 2.2 m² s⁻², ν = 0.002 m² s⁻¹ , and ν' = 0.10 m² s⁻¹, model closely those of the snow experiment. Moreover these parameters form a unique set in which variation in one parameter will degrade the modeling results, whatever adjustments may be made in the other two coefficients. Additional validity to the values of these parameters is obtained from other experiments. The work of Maeno and Nishimura (1979) and of Maeno and others (1980) on snow suspended by air to form a fluidized bed, produced measurements of kinematic viscosity of the order of 0.001 m² s⁻¹ for incompletely fluidized snow. Bûcher and Roche (1946), in measuring the frictional resistance of hard wet snow for speeds between 0.2 and 2.4 m s⁻¹, found that the linear fit to their data yielded a constant of proportionality of 475 (N-S)/m. If it is assumed that there was a 2 mm layer of granulated snow of density 300 kg m⁻³ between the sliding surfaces and that the velocity gradient was linear in this region, then the viscosity in this layer would be about 0.003 m² s⁻¹. Similar tests by Dent and Lang (1982), using hard sintered snow over the velocity gradient range 50 to 300 m s⁻² yielded a viscosity coefficient of 0.004 m² s⁻¹ and a T_Q/P value of 1.8 m² s⁻². These values are for a very narrow range of slow speeds and probably differ at higher speeds, but do serve as order of magnitude values.

A last observation is that the tangential boundary condition used in the modeling at the bottom boundary was the no-slip condition. The quality of the modeling results lends credence to the hypothesis that this boundary condition is appropriate for flowing snow.

In carrying out the computer modeling the time-step between calculation cycles was chosen such that the maximum distance traveled by any part of the fluid was less than 0.1 of a cell dimension. Using the cell dimensions previously described and this tine-step criterion, no numerical instabilities were encountered for the range of parameters involved in this modeling. To generate the results exhibited in Figure 7, each modeling run required about 1 000 calculation cycles, taking, for the 1 400 cell computational system, about 30 ruin of CPU time on the system used.

Conclusions

For snow flow in the speed range <20 n s⁻¹ the biviscous model has provided satisfactory results. The overfall motion of the snow as depicted in the motion of the leading edge and the final distribution of the depth of the debris were well simulated. In addition, details noted in the snow tests were reproduced by the computer model. Quantitatively, the velocity as a function of depth was accurately modeled. Qualitatively, the formation of the boundary layer can be seen in the time-seauence particle plots in Figure 7. The particles near the front of the flow at the bottom are retarded as the upper part of the flow proceeds over them. These particles are seen to be left in a layer along the bottom boundary, just as the dye placed originally in the front of the flow in the snow tests was seen to be distributed as a layer over the entire runout area (Dent and Lang 1982), Examination of the motion of the marker particles in the upper portions of the flow shows that little deformation is taking place in this region. This motion is confirmed by observations (Dent and Lang 1982). Another aspect of the Flow seen in both the snow tests and the computer modeling is the surging motion of the leading edge. Although it is not shown clearly in the particle plots, the front of the flowing mass was continually breaking over the slower moving flow near the surface. This motion showed up most strikingly by monitoring the velocity at the leading edge. It was found that this velocity was not a smooth function of time but exhibited large variations around the average decaying velocity. Figure 8 is a plot generated by the computer at the time of execution showing the speed of the leading edge versus time. The speed plot shows this surging motion clearly. This motion was also seen when reducing the velocity of the leading edge in the tests of snow flow from 16 mm film. It showed up as anomalous measurements of the velocity of the leading edge at sporadic times in the flow. It could also be seen viewing the motion-picture film, as surging or jetting of the leading edge, rauch like the motion of water waves shoaling on a beach after breaking.

Fig. 8. Velocity of the leading edge versus time (computer).

At speeds above 20 m s⁻¹, much conjecture still exists as to the behavior of flowing snow. Mellor (1963), Perla (1980), and others have speculated that the flowing snow enters a turbulent flow regime at high speeds. This transition point must be a function of speed and type of snow In the avalanche. So far, there is no documentation on when avalanching snow enters a turbulent flow regime. The motion of the leading edge of the avalanche or the powder dust cloud is freguently cited as evidence of turbulence. However, for those avalanches with a central core, the motion of that mass of material does not necessarily have to be turbulent. As seen in the experiments on flowing snow (Dent and Lang 1982), the largest velocity gradients, and hence the greatest dissipation, is at the base of the flowing snow. It is also at this location that gravitational forces produce the largest normal stresses. In order for this area to become turbulent these normal stresses must be overcome by granular interaction. This may happen if the avalanche speed is fast enough, but this will be the last area to become turbulent. The conditions necessary for this to happen are unknown. However, for speeds below 20 m s⁻¹, no evidence of turbulent motion was observed in the tests on snow flow (Dent and Lang (1932)). Marker dye placed in the snow was not seen to diffuse, but remained in localized regions, deforming by what appear to be streak-lines. As well as turbulence in avalanches, there must also exist other velocity-squared forces. Air drag, ploughing, and entrainment are all effects that at some point need to be considered.

The exact constitutive relation for flowing snow is extremely complicated. On physical grounds, a Bingham-type relation seems reasonable at low speeds, with the addition of velocity-squared mechanisms at higher speeds. The transition point is unknown. Investigations of flow velocities and density distributions in higher speed avalanches need be made to clarify this transition and the mechanics involved. Once these data are collected, flow mechanisms may be evaluated. The computer model described in this paper is easily modified to simulate other flow laws. It is a simple matter to generalize from a biviscous formulation to a tri-viscous formulation or a flow law that involves more viscosities. In this way the velocity-squared forces could also be approximated.

The computer simulation methodology, particularly the multi-viscosity approach, has proved to be suitable for solving problems concerned with flowing snow. The excellent internal consistency shown by the bi-viscous modeling of the tests on flowing snow inspires a great deal of confidence in the method. In addition, the ease with which the model may be generalized to include more complex constitutive laws indicates a very good prospect for its use as additional information about the mechanics of the flowing material is learned.