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/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

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. Reference Norem, Irgens and SchieldropNorem and others (1987) proposed the following slip condition

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

where

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

where maximum velocity v₁, slip velocity v₀, 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 dv_xdy 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, 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

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₀:

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 v₁(X,t), v₀(X,t), and h(X,t).

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.

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 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/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

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 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⁻¹ 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.

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⁻¹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, 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 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.