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 Reference Norem, Kvisterøy and EvensenNorem 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 (Reference Norem, Irgens and SchieldropNorem 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
Reference Norem, Irgens and SchieldropNorem 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)
where
a is the cohesion (N/m2), b is the material friction parameter, pe is the effective pressure (N/m2), pu is the pore density (N/m2), ρ is the density (kg/m3), m is the shear viscosity (m2);
V1, V2 = normal stress viscosities (m2).
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 vx = vx(y) and vy = vz = 0 (Fig. 2), where ρ is assumed constant, the stresses at a height, y, are given by
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 p1 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. Reference Norem, Irgens and SchieldropNorem and others (1987) proposed the following slip condition
where v0 = 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, v0, the maximum velocity, v1 and the velocity distribution
where
Non-Steady Flow in Run-Out Zone
Introduction
In the run-out zone of a snow avalanche, the linear flow velocity, vx, the flow height, h, and the plug-flow thickness. hp. will be functions both of time and of position along the avalanche path. The linear flow velocity will be represented by the maximum velocity, v1, together with the slip velocity, v0. The solution of the flow problem consists of introducing four equations for the four functions, h, hp, v1 and v0, 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, hp, 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
where maximum velocity v1, slip velocity v0, and flow height h are functions of the horizontal coordinate, X, and of time, t, such that
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 dvxdy need be considered. The Cauchy equation of motion in the y-direction will then reduce to
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
and
The shear and normal stresses are now given by Equation (15) together with the following constitutive equations
Equation (17) for the effective pressure, pe, 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
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, v1 and the slip velocity, v0:
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.
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 v1(X,t), v0(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 Q1, …, Qn 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, Q1 and Q2, is always taken to be straight. The slope of this segment and the slope at the end point, Qn, have to be specified as input data. Between any two neighbouring points, Q2, …, Qn, 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, P1, …, Pm. A computer sub-routine evaluates the coordinates (Xi, Yi), the slope, ϕi, and the curvature, κi, at each station, Pi.
Computer program
To start the computation, an “avalanche” is introduced covering a chosen number of sub-segments on the straight segment between Q1 and Q2. The flow height, h, has the same constant value on all sub-segments, and the avalanche is assumed to have obtained its terminal velocities, v0 and v1, 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: v0(X,t); v1(X,t); h(X,t).
S.2 = No-slip: v0 = 0; v1(Χ,t); h(X,t).
S.3 = Uniform profile: v0 = v0; v1(X,t); h(X,t).
S.4 = Constant height: v0 = 0; v1(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 Pf fills the down-stream sub-segment and when the accumulated volume exceeds the value hfΔxf+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 Pt+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, v1 and v2, and for the calculations we set values for m, m/v2, and v1/v2. The term v1/v2 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 v1/v2 values.
The term m1/v2 has been studied theoretically by Reference Savage and SayedSavage (1984), and experimentally by Reference Savage and SayedSavage and Sayed (1984), who found it to be dependent mainly on the coefficient of restitution. Reference Norem, Irgens and SchieldropNorem and others (1987) have assumed the value of m/v2 for slow particles to lie in, or close to, the range 0.8–1.0, and also shown that for slope angles, ϕ, where ϕ = tan−1m/v2, the expression for the terminal velocity for no-slip conditions, could be reduced to
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, v1. 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 m2.
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−1, the value of s must lie in the range of 0.5–2.0 kg m−3. The friction parameter, b, has the same characteristics as the static internal friction parameter of granular materials. According to Reference Savage and SayedSavage 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. Reference Lang and DentLang 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.
Reference Lied and Bakkehøi.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−1 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 Reference Norem, Kristensen and TronstadNorem 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−1when 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−1 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, Reference MartinelliMartinelli (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 Reference VoellmyVoellmy (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
although Reference Hungr and McclungHungr 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 m2 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 Reference VoellmyVoellmy model (1955) or the Reference Perla, Cheng and McClungPerla, 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 Reference Trunk, Dent and LangTrunk and others (1986) and by Reference Hutter and SavageHutter 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 Reference Trunk, Dent and LangTrunk 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 Reference Hutter and SavageHutter 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.
