Simulation of Snow-Avalanche Flow in Run-Out Zones

Abstract

The authors of this paper have previously, in 1987, outlined the development of a set of constitutive equations representing the behaviour of avalanching snow, and the resulting equations for fully developed steady shear flow. The present paper considers the development of non-steady two-dimensional shear flow together with the use of a finite-difference programme to calculate snow-avalanche velocities and flow heights in the run-out zone. The numerical results thus obtained are compared with full-scale experimental data. These comparisons indicate that front velocity and run-out distance are simulated well by the model, and that the predicted snow deposits are fairly well in agreement with those actually observed.

Introduction

Since 1981, the Norwegian Geotechnical Institute (NGI) has carried out full-scale experiments in the Ryggfonn avalanche path in western Norway, recording avalanche velocities and impact pressures. The practical arrangements employed and the results obtained have been reported by Norem and others (1985, 1986, 1988a, b); evaluation of the results of these experiments has clearly indicated the need for a better physical understanding of avalanche flow, especially in the lower part of the avalanche path and in the run-out zone.

A theoretical study of avalanche flow was started in 1985, with the following aims:

To give a physical description of the mechanics of the flow, and the retarding mechanism.

To identify the physical parameters which are most important in defining avalanche flow.

To predict the velocities, flow heights, and stresses in the lower part of an avalanche path, and also the run-out distance, with a reasonable degree of accuracy.

The preliminary results of this project have already been given by the present authors (Norem and others, 1987) in a report which includes the development of a set of constitutive equations representing the behaviour of avalanching snow, and the resulting equations for fully developed steady shear flow. The present paper deals with the development of non-steady two-dimensional shear flow, and with the use of a finite-difference programme to calculate snow-avalanche velocities and flow heights in the run-out zone. The numerical results are compared with recorded data from the Ryggfonn project.

Introduction to Continuum Model

Constitutive equations

Norem and others (1987) proposed that, where viscoplastic behaviour is predominant, snow in avalanches should be modelled as a granular material. In a simple steady shear flow, the general constitutive equations yield the total stresses, τ_xy, yz, zx (Fig. 1)

(1)

(2)

(3)

(4)

where

a is the cohesion (N/m²), b is the material friction parameter, p_e is the effective pressure (N/m²), p_u is the pore density (N/m²), ρ is the density (kg/m³), m is the shear viscosity (m²);

= dv_x/dy = shear velocity (1/s),

V₁, V₂ = normal stress viscosities (m²).

The shear stress, Equation (1), is caused by a combination of cohesion, Coulomb friction, and viscosity.

Fig. 1. Definition of normal and shear stresses, and velocity distribution.

Fig. 2. Definition of steady-flow geometry.

The normal stresses, σ_y, Equation (2), are divided into three separate parts: effective pressure, pore pressure, and a dispersive pressure. The viscous part of the shear stress and the dispersive pressure are assumed to be due to collisions of particles moving relative to each other.

Steady shear flow

In a fully developed steady shear flow with a velocity field given by v_x = v_x(y) and v_y = v_z = 0 (Fig. 2), where ρ is assumed constant, the stresses at a height, y, are given by

(5)

(6)

where ϕ is the angle of inclination; g is the acceleration due to gravity. At the surface of a dense flow, y = h, the shear stress is zero, and the pressure is equal to the atmospheric pressure p₁ which is assumed to be equal to the pore pressure within the flowing snow.

For most avalanches, some slip velocity at the bed surface is likely to exist. Norem and others (1987) proposed the following slip condition

(7)

where v₀ = slip velocity; s = a constant, whose value is dependent on the nature and roughness of the material of the solid surface. From Equations (1), (2), (5), (6), and (7) we obtain expressions for the slip velocity, v₀, the maximum velocity, v₁ and the velocity distribution

(8)

(9)

(10)

where

(11)

Non-Steady Flow in Run-Out Zone

Introduction

In the run-out zone of a snow avalanche, the linear flow velocity, v_x, the flow height, h, and the plug-flow thickness. h_p. will be functions both of time and of position along the avalanche path. The linear flow velocity will be represented by the maximum velocity, v₁, together with the slip velocity, v₀. The solution of the flow problem consists of introducing four equations for the four functions, h, h_p, v₁ and v₀, of time and position along the avalanche path, and then solving these equations numerically. So far we have not formulated a plug equation, and we therefore leave out the plug-flow height, h_p, and the cohesion parameter, a, in the following discussions. A more detailed exposition of the method of solution has already been presented by Irgens (1988).

Fig. 3. Geometry of avalanche path. Velocity profile.

Geometry of avalanche path: velocity profile

The avalanche is considered as a two-dimensional flow-in a vertical plane, as shown in Figure 3, for which x and y are local Cartesian coordinates originating from any point (X,Y) along its path. The unknown velocity distribution for the non-steady flow case is taken to be represented by the expression

(12)

where maximum velocity v₁, slip velocity v₀, and flow height h are functions of the horizontal coordinate, X, and of time, t, such that

(13)

The velocity profile Equation (12) is identical in form with the steady-shear flow profile Equation (10) for non-cohesive materials.

In accordance with standard thin shear-layer assumptions, the velocity component in the y-direction can be neglected. Furthermore, it may be shown that in general constitutive equations only terms involving the gradient dv_xdy need be considered. The Cauchy equation of motion in the y-direction will then reduce to

(14)

in which the terms on the left-hand side represent a normal acceleration due to the curvature, κ, of the avalanche path. Integration of Equation (14) yields

(15)

and

(16)

The shear and normal stresses are now given by Equation (15) together with the following constitutive equations

(17)

(18)

(19)

(20)

Equation (17) for the effective pressure, p_e, is obtained by comparison of Equation (15) with the constitutive equation for σ_y.

Fig. 4. Differential control volume.

Momentum and continuity equations

The equation of motion in the direction of linear flow is replaced by a one-dimensional momentum equation. This is obtained by application of the von Kármán integral method, the basic idea being to develop the momentum equation for a differential fixed control volume, hdx, as shown in Figure 4. The momentum theorem and the principle of conservation of mass, applied to the control volume, yield the following integral equations

(21)

(22)

Differential equations

The shear stress at the bed surface, y = 0, is assumed to satisfy the slip condition of Equation (7). Equating this expression for the shear stress with that of the constitutive Equation (18), we get the following ratio between the maximum velocity, v₁ and the slip velocity, v₀:

(23)

The velocity profile Equation (12), the normal stress Equation (19), and the shear stress Equation (7), or alternatively Equation (18) for the no-slip case, are substituted into the integral momentum and continuity Equations (21) and (22). The final forms of the integrated equations are the following partial differential equations.

(24)

(25)

The flow problem in the run-out zone is thus simplified to the solution of the three Equations (23), (24), and (25) for the unknown functions v₁(X,t), v₀(X,t), and h(X,t).

Numerical Solution

Definition of the model avalanche-path profile

The two-dimensional run-out zone of the avalanche path is presented in Figure 3. The Cartesian coordinates (X,Y) of n characteristic path points Q₁, …, Q_n are chosen to describe the main topography of the path. In order to obtain realistic intial values for the velocity, the segment between the two first points, Q₁ and Q₂, is always taken to be straight. The slope of this segment and the slope at the end point, Q_n, have to be specified as input data. Between any two neighbouring points, Q₂, …, Q_n, third-order parabolas are introduced; these parabolas satisfy continuity of slope and curvature at the Q points, and an actual avalanche path is thus replaced by a cubic spline. The segments between the Q points are further subdivided into an arbitrary number of sub-segments by the station points, P₁, …, P_m. A computer sub-routine evaluates the coordinates (X_i, Y_i), the slope, ϕ_i, and the curvature, κ_i, at each station, P_i.

Computer program

To start the computation, an “avalanche” is introduced covering a chosen number of sub-segments on the straight segment between Q₁ and Q₂. The flow height, h, has the same constant value on all sub-segments, and the avalanche is assumed to have obtained its terminal velocities, v₀ and v₁, computed from Equations (8) and (9) on entering the straight segment. Any other set of initial conditions, however, may also be introduced into the computer programme.

In order to simplify comparison with other avalanche simulations, four different sub-routines are written for the solution of the flow Equations (23), (24), and (25). These are

S.1 = Slip: v₀(X,t); v₁(X,t); h(X,t).

S.2 = No-slip: v₀ = 0; v₁(Χ,t); h(X,t).

S.3 = Uniform profile: v₀ = v₀; v₁(X,t); h(X,t).

S.4 = Constant height: v₀ = 0; v₁(X,t); h = constant.

The partial differential Equations (24) and (25) are solved by a finite–difference scheme, with central differences in the stream-wise direction, and by a fourth-order Runge–Kutta procedure with respect to time. The finite–difference scheme is Eulerian, requiring a special procedure to make the avalanche progress along the path. Figure 5 illustrates this procedure; the volume of snow passing through the front section at station P_f fills the down-stream sub-segment and when the accumulated volume exceeds the value h_fΔx_f+1 it is assumed that the avalanche front has advanced one sub-segment. A similar procedure is applied to the tail of the avalanche; when the accumulated volume of snow passing through the section at station P_t+1 exceeds the volume of snow originally on the tail segment, it is assumed that the tail segment has emptied, and that the tail has therefore advanced one sub-segment.

Fig. 5. Flow through the elements. The indices, t and f, refer to tail and front of avalanche.

Discussion of the Physical Parameters

Viscosities and friction parameters

The sub-routine programme for no-slip conditions is based on the selection of the three viscosities, m, v₁ and v₂, and for the calculations we set values for m, m/v₂, and v₁/v₂. The term v₁/v₂ is assumed to be equal to 10, in accordance with experimental results for non-Newtonian fluids, although for our project the numerical results turn out to be insensitive to the selection of v₁/v₂ values.

The term m₁/v₂ has been studied theoretically by Savage (1984), and experimentally by Savage and Sayed (1984), who found it to be dependent mainly on the coefficient of restitution. Norem and others (1987) have assumed the value of m/v₂ for slow particles to lie in, or close to, the range 0.8–1.0, and also shown that for slope angles, ϕ, where ϕ = tan⁻¹m/v₂, the expression for the terminal velocity for no-slip conditions, could be reduced to

(26)

and that this expression predicts for steep slopes that the effect of the frictional parameter, b, is reduced, and that of the material parameters only the shear viscosity, m, influences the terminal velocity, v₁. This means that if the terminal velocity in the steep part of the path, and the flow height, are known, the value of m can be calculated. A probable range of values of m is 0.001–0.01 m².

The only dynamic parameter used in the sub-routine programme for uniform profile is the dynamic friction parameter, s. To obtain terminal velocities of 30–70 m s⁻¹, the value of s must lie in the range of 0.5–2.0 kg m⁻³. The friction parameter, b, has the same characteristics as the static internal friction parameter of granular materials. According to Savage and Sayed (1984), the value of b is very close to the tangent of the internal friction angle and is dependent on the size, roughness, and moisture content of the particles. Lang and Dent (1983) found the value of b to be close to 0.42. The authors of this paper also refer to experiments by Inaho, who found values within the range of 0.40–0.55 for b.

Lied and Bakkehøi (1980) have surveyed major avalanches in Norway and have found that the angle from the starting point to the run–out point is on average 33°, although values as low as 18° and as high as 49° have been found. The limits for b may thus be between tan 18° (= 0.32) and tan 49° (= 1.15). The numerical results cited show that, for avalanches having the same terminal velocity, the run-out distance will vary considerably with the values of b. Our calculations are based on b values of 0.3 and 0.4.

Length and height of avalanche

Data from the Ryggfonn experiments show that impact pressure rises very rapidly when a dense flow hits a solid structure, and that maximum pressure will be obtained 3–10 s after the first impact. Pressure is then reduced gradually as the avalanche tail passes the structure.

Radar experiments by Gubler (1987) indicate that the maximum velocities within an avalanche are found close to its front, and that the highest velocities are seldom found for a period of more than 10 s at any one point. The findings of the radar experiments and the Ryggfonn data thus coincide.

Front speeds of 30–60 m s⁻¹ are typical for major avalanches, and lengths of 100–400 m seem to be reasonable estimates. An exact value for the flow height is difficult to determine, as the cross-section of an avalanche is usually of an irregular shape and there are no clear distinctions between the dense flow and the accompanying snow clouds. Impact-pressure recordings from Ryggfonn indicate a height of dense flow of between 1 and 3 m. The impact pressures on the upper load cells of the recorder also indicate that maximum flow heights are found only for a short time and are sited in the front of the avalanche. The assumption that an avalanche initially has a given length and constant height and velocity in each sub-segment are thus simplifications. Probably, a more accurate assumption would be that of a varying flow height distribution, with the maximum height fairly close to the front of the avalanche. The model allows such height distributions have been taken into account but, due to the limited number of measurements available, they have not been considered in this paper.

Numerical Results

Approaches to practical problems io snow-avalanche engineering

In evaluating the safety and protection in populated areas, and considering structures in the path of an avalanche, one needs to know:

The run-out distance for a particular avalanche path, as a function of a given size of avalanche.

The effect of varying path profiles in the run-out zone, with a particular emphasis on the stopping effect of collecting dams.

Avalanche velocity at any point in the run-out zone.

Impact pressures on structures, and shear stresses at ground level.

In order to evaluate the accuracy and usefulness of the proposed model, we have used it as the basis for comparison with some of the numerical calculations for the Ryggfonn avalanche profile, and also for several idealized avalanche profiles.

Comparisons with Ryggfonn data

The numerical results generated by the model are compared with data from three major dry-snow avalanches in the Ryggfonn area. The experimental arrangement in Ryggfonn consisted of a 25 m Y-tower, a concrete structure

Fig. 6. Time of impact and surveyed avalanche deposits compared with numerical results of simulated avalanches.

with load cells, several geophones, and a 15 m high collecting dam with a mast on top of it. The profile of the lower section of the avalanche track, and the location of the arrangement described, are shown in Figure 6. The complete sets of recorded data from the three avalanche events have been presented in Norem and others (1986, 1988a, b). Figure 6 shows the time of impact at the different experimental sites, and the surveyed avalanche deposits. It also shows the calculated flow of three simulated avalanches in the lower section of the track, and the calculated snow deposits on both sides of the dam. The height and length of the model avalanches were selected by use of the impact-pressure recordings, and the friction parameter, b, was set as b = 0.4 in every case. The shear viscosity, m, was selected to obtain similar front velocities between the Y-mast and the load cells for both the simulated and the recorded avalanches.

The numerical results for avalanche No. 1 show that the time calculated as being necessary for it to reach the top of the dam is slightly less than the actually observed time for the avalanche on 13 February 1985. The snow cloud of this avalanche had a front velocity of 15 ms⁻¹when passing the dam, but only a small amount of the dense flow did pass the dam. The calculations show that the model avalanche stops only 2.0 m below the top of the dam, which is in reasonable agreement with the experimental observations. However, the calculated snow deposits extended somewhat higher than the surveyed snow depths, and the tail of the avalanche was not included in the model at all.

The avalanche of 11 April 1988 had a speed of 43 ms⁻¹ in the area of the Y-mast. Only a small proportion of the dense flow material passed the dam, and this

Fig. 7. Calculated run-out distances for paths with different run-out slopes.

material was finally deposited up to 55 m behind the dam. The numerical results for avalanche No. 2 agree very well with the observed results for the front velocity, the run-out distance, and the height of the deposits, except at the tail of the avalanche.

The avalanche of 28 January 1987 was a major natural avalanche, having a run-out of between 100 and 150 m behind the top of the dam. Unfortunately, there had been another large avalanche on the previous day and the surveyed snow depths are therefore the sum for the two avalanches. The major part of the deposit, however, is considered to derive from the later avalanche, since the duration of impact of the preceding day’s avalanche was only short. It is thought probable that the effective dam height was reduced when the avalanche occurred on 28 January, because of previous avalanche deposits; however, since there is so little information about these snow depths, we have based our calculations on the original terrain profile.

The calculations for simulation avalanche No. 3 have been made for an avalanche length 50% longer than the previous ones and for an initial velocity 3% higher than for avalanche simulation No. 2. The numerical results correspond well both with the recorded times of impacts and with the snow deposits in front of the dam. The main deviations are found on the lee side of the dam, where the calculated snow depths are between 1 and 2 m greater than predicted, and the run-out distances are from 20 to 60 m shorter.

Effect of slope at end of the run-out zone

The profiles of six idealized run-out zones are shown in Figure 7. The first section is the same for all profiles, consisting of a steeper and a more gentle part, the final sections of the run-out zones have slopes varying gradually from zero at the start to an end value of between 5.7° and −33.7°. The numerical results show that the run-out distance falls off considerably when the end value of the slope varies from 5.7° to −5.7°, and that for slopes steeper than −16° the reduction in run-out distance is slight; in this case the run-up height is a more significant factor.

Also presented in Figure 7 are the calculations made for the sub-routine programme, for constant height and with the same input data. The calculated run-out distances vary only slightly with the slope at the end of the run-out zone, because in this case the retardation is calculated by using the averaged slope angle over the whole section covered by avalanching snow.

Fig. 8. Calculated run-out distances for varying path profiles.

Effect of different path profiles

When evaluating the run-out distance of avalanches, it is necessary to take into account the effect of the shape of the avalanche profile in the run-out zone. For instance, Martinelli (1986) observed different run-out distances in what he called hockey-stick profiles compared with those having a more parabolic shape.

The results presented in Figure 8 are based on simulations for six different idealized avalanche paths. The profiles are named straight line, parabolic, hockey-stick, concave, 5 m wave, and 10 m wave, respectively. The first and final sections are common to all profiles, while the shape differs for 200 ≤ x ≤ 450 m. The numerical results for the sub-routine programme with varying heights indicate

Fig. 9. Calculated run-up heights on collecting dams, as a function of velocity at front of dam.

that the run-out distance varies considerably with the shape of the path profile. The longest run-outs, x = 478–484 m, are found for profiles that start to level off at x = 200 m. The shortest run-outs are found for the hockey-stick profile, for which x = 421 m, where the avalanche enters the slopping area with a considerable speed. It is interesting to note the differences in the calculated run-out distances and snow deposits for the 5 m wave and the 10 m wave profile. The height differences in the 5 m profile have obviously not been significant in altering the flow, whilst the second peak in the 10 m profile has stopped most of the snow, with snow deposits up to 8.0 m deep. Only a minor part of the avalanche passed this peak, and it then stopped in the first gentle-slope section.

Run-out distances have also been calculated with the sub-routine programme for constant height. The calculations show that this programme predicts almost the same run-out distances as are observed on gentle slopes, but is less sensitive to variations in path profiles. The standard deviations for the calculated run-outs are 27.7 and 10.7 m for varying and constant height, respectively.

Calculated run-up heights on collecting dams

Collecting dams have for a long time been used to protect houses and roads against avalanches. An equation derived by Voellmy (1955) has been the basis for most practical calculations. When the effect of dynamic friction is omitted, and only the loss of potential energy and Coulomb friction are taken into account, the Voellmy equation yields the following formula for the run-up height

(27)

although Hungr and McClung (1987) have presented a calculation for the run-up height that gives approximately twice the value for ΔH of that given by Equation (27).

Figure 9 presents the numerical results of several calculations made by the present model. Most test runs have been made for the profile with a −33.7° slope at the path end, and these are presented in Figure 7, but some were also made for the Ryggfonn profile, and these are shown in Figure 6. The calculations are made for a constant flow height of 2.0 m and for varying avalanche lengths. The initial velocity is varied by selecting different shear viscosities. Calculations indicate that there are distinct differences in the run-up heights for avalanches having a short length (L < 100 m) compared with longer ones. For instance, the calculated run-up heights for a shear viscosity m = 0.006 m² and L = 96, 144, and 192 m, the respective run-up heights are 7.0, 10.4, and 11.5 m. However, when the results are presented, as in Figure 9, with run-up height and front velocity at the front of the dam along the axes, the difference is not so distinct. This is because the shorter avalanches have a lower front velocity when entering the dam. The curve in Figure 9 is the best fit for the results for the avalanches with L = 144 m, and is in close agreement with Equation (27).

The calculation made for the Ryggfonn profile shows higher values than those obtained for the profile shown in Figure 7. This is probably due to the presence of a steep slope just in front of the dam area for the Ryggfonn profile. The velocity difference between the front and the rear will, in this case, be higher than expected and that in turn will influence both the velocity and the run-up height of the front element. Preliminary results of calculating the effect of collecting dams indicate that the necessary height of the dam is dependent not only on the Coulomb friction and the front velocity, but also on the volume of the avalanche and on the shape of the area in front of the dam.

Evaluation of Present Simulation Model

The sub-routine programme for constant flow height is the simplest of the four programmes developed. When selecting a very short avalanche length, this programme relates to the Voellmy model (1955) or the Perla, Cheng, McClung (PCM) model (1980). On open slopes, the programme for constant height may give the same run-out distance as the programme for varying height, the latter programme being less sensitive to the shape of the avalanche path, and its results tending to give unrealistically high run-up heights on collecting dams. Thus we conclude that the imposition of constant height on the model gives rise to inaccurate data about avalanche flow.

There are three sub-routine programmes for varying flow height: one for a no-slip condition, one for a uniform velocity profile, and one which is a combination of these two. Test calculations have shown that when selecting similar initial conditions, the three programmes give almost the same results; the run-out distances are virtually identical, but the avalanche deposits are somewhat less for the uniform profile programme than those observed. Our knowledge of the avalanche-velocity distribution and of the magnitude of the slip velocity is very limited, and because of this we suggest that at the present time it is impossible to make a choice between these three alternatives. The model has the capacity for further refinements, especially of the boundary conditions at the bed. These refinements will need to take into account the effects of entrainment from a soft bed, and of the loss of mass if there is a no-slip condition at the bed. Other recognized weaknesses are the two-dimensional flow assumption, and the failure properly to simulate the tail of the avalanche. However, for major avalanches with a considerable width, the assumption of two-dimensional flow is probably less constricting than other simplifications made. The simulation of the tail can easily be included in the computer programme to improve the agreement between the shapes of simulated and observed snow deposits. These refinements of the model will increase computing time, and are assumed by us to be of only minor practical interest.

In recent years two interesting models for simulating rock slides have been presented by Trunk and others (1986) and by Hutter and Savage (in press). Both models are based on theories of general granular flow and offer an opportunity to take into account the effect of the volume of snow involved in an avalanche and the variation in flow height along the path profile.

The model of Trunk and others (1986) needs two viscosity parameters for the material properties, and the flow is calculated in a grid system. Flow height and velocity are calculated as functions of both time and of distance along the path. It is not known whether the model has been adapted for snow avalanches.

The model by Hutter and Savage (in press) also gives flow velocity and height as functions of time and distance. The only material parameter in their model is a friction parameter, and the variation of flow height is controlled by passive and active earth pressures. It is our opinion that snow avalanches attain a terminal velocity, and that this terminal velocity can be regulated only by some kind of a kinetic friction parameter or by viscosity. The Hutter and Savage model is thus probably more generally applicable to major rock avalanches, that will not attain a terminal velocity, than to snow avalanches. It must at this point be admitted that their model probably has a more realistic physical simulation of the triggering and start of granular flows than is implied by our model, and it may therefore be appropriate to combine some elements of these two models in order to improve the calculations for the whole avalanche track. However, the programmes are not available to the authors, and it has thus not been possible for them to make a numerical comparison between the different models at the time of writing.

Conclusions

Based on published data for snow avalanches, practical experience gained by the NGI group, and data recorded from the Ryggfonn experiments, we assume the present model of snow-avalanche flow to have an accuracy sufficient for practical use. The velocities along the avalanche path, and also the run-out distances, seem to be especially well simulated, and the experimentally recorded snow deposits are fairly well estimated by the model. The model probably takes into account most of the parameters which are important in describing avalanche flow, since the parameters used all have a direct physical explanation and can be measured experimentally. It is thus our hope that the model can be adapted to apply to the granular flow of other materials, and to model experiments. The selection of parameter values is currently based on qualified guessing, although it is hoped that further research, both with full-scale and with model experiments, will increase our knowledge of granular flow. The research results will also be helpful in reaching an understanding of the importance of the presented parameters, and in selecting more reliable values for these parameters. The next step in the Ryggfonn project will be the analysis of the dynamic pressures recorded as being exerted on different constructions, and we have assumed that the model will be a useful practical tool in the evaluation of snow-avalanche impact pressures.

Acknowledgements

This project is financed mainly by the Royal Norwegian Council for Scientific and Industrial Research, for which the authors express their gratitude. The authors are also grateful to the Norwegian State Power Board and NGI, who have financed the Ryggfonn experiments, making it possible to compare simulations with data from real avalanches. Finally, we also thank the staff at NGI for fruitful discussions.