List of Symbols
Centre-of-Mass Model
-
g x, g y, g z, g Gravitational force per unit mass
-
M/D Mass-to-drag ratio
-
R Radius of curvature of the centre-of-mass path line along the wall
-
r Vertical run-up height above the base line of the deflecting dam, measured in a vertical cross section perpendicular to the horizontal projection of the base line of the deflecting dam
-
ν Velocity of centre-of-mass
-
ν T Centre-of-mass velocity immediately before impact
-
y Run-up height above the base line of the deflecting dam, measured along the terrain in a cross section perpendicular to the base line of the natural deflecting dam
-
α s Inclination of line of steepest descent along the upper wall of the dam
-
β Angle of inclination of plane terrain
-
γ Angle between the centre-of-mass path tangent line and the base line
-
μ Dry-friction coefficient
-
ψ Angle of inclination of the wall relative to the terrain
One-Dimensional Continuum Model for Three-Dimensional Avalanche Flow
-
a Cohesion parameter
-
b XZ Centrifugal force per unit mass in the horizontal plane
-
b xy Centrifugal force per unit mass in vertical tangent plane
-
g Gravitational force per unit mass
-
g e Effective gravitational force per unit mass
-
h Flow height
-
m Shear viscosity
-
m 1 m 2 Visco-clasticities
-
n Material power law exponent
-
P i Station points on the centre line
-
P e Effective pressure
-
p u Pore pressure
-
Q i Path points defining the centre line of the avalanche terrain
-
R Radius of computational cross sectional profile
-
ν a Average velocity through flow cross section
-
ν x, ν y, ν 1 Velocity components
-
w Computational flow width
-
XYZ, xyz Cartesian coordinates
-
α X Slope of horizontal projection of centre line
-
θ Angle of inclination of the free surface of the avalanche
-
k XZ Curvature in the horizontal plane
-
k xy Curvature in the vertical tangent plane
-
μ Dry-friction coeffcient
-
ρ Mass density
-
σ, τ Normal stress and shear stress
-
ø Slope of centre line in the vertical tangent plane
1. Introduction
Increased human activity in mountain regions, deforestation from pollution, forestry and ski resorts, as well as a reduced acceptance of living in regions exposed to snow avalanches have caused a growing need for protection against avalanches. Such protection is more and more often obtained by constructing retaining dams to influence the dense-snow avalanche course. Better knowledge and understanding of terrain deflection of dense-snow avalanches will improve the design of deflecting dams.
The first attempt to formulate a general theory of dense-snow avalanche motion was made by Reference VoellmyVoellmy (1955) and this theory is still widely used. Both statistical and comparative models for run-out distance computations, as well as dynamic models for avalanche motion simulations, are now developed. However, no universal model has so far been made. The dynamics of avalanches are complex, involving both fluid, particle and soil mechanics. The limited amount of data available from real events makes it difficult to evaluate or calibrate existing models. Often several models with different physical descriptions of avalanche movement can all fulfil the deficient observational records.
In spite of the uncertainties with which the existing models are encumbered, simulations of avalanches impacting retaining dams at oblique angles of incidence have to take into account that real avalanche flows are three-dimensional. Neither the run-up heights pioneered by Reference VoellmyVoellmy (1955) nor the leading-edge models by Reference Hungr and McClungHungr and McClung (1987) or Reference Takahashi and YoshidaTakahashi and Yoshida (1979) (Reference TakahashiTakahashi, 1991) contain deflection effects. This is also the case for the paper by Reference Chu, Hill, McClung, Ngun and ShcerkatChu and others (1995) on experiments on granular flows to predict avalanche run-up. Reference NohguchiNohguchi (1989) has developed a three-dimensional dense-snow avalanche centre-of-mass model based on the equations of Reference VoellmyVoellmy (1955), while Reference Sassa and BonnardSassa (1988) has developed a geotechnical quasi-three-dimensional continuum model. Reference Lang and LeoLang and Leo (1994) developed a quasi-three-dimensional dense-snow avalanche model. However, according to the originators it is still unknown whether the model can represent naturally occurring events.
Effect of natural deflecting dams on reported avalanches are described by Reference Harbitz and DomaasHarbitz and Domaas (1997). The observations indicate that the height difference between the gully floor and the upper limit of extension on deflecting terrain formations might exceed 60 m for large avalanches (estimated volume ¾100 000 m3) reaching a velocity of more than 40 m s-1 in the avalanche track.
In this paper, we present two models of snow avalanches impacting on retaining dams at oblique angles of incidence. The first is a centre-of-mass model based on a Voellmy type of resistance force. In the second model, we extend the simulation model of Reference Norem, Irgens and SchieldropNorem and others (1989), for one-dimensional avalanche flow parallel to a vertical plane, to avalanche flows following a three-dimensional channel with varying width. The material model used was first presented by Reference Norem, Irgens and SchieldropNorem and others (1987) and a discussion on the importance of the physical parameters has been given by Reference Irgens and NoremIrgens and Norem (1996).
2. Centre-of-Mass Model for Avalanche Motion Along the Side of a Retaining Dam
A centre-of-mass model for avalanche motion along the side of a retaining dam was developed by B. Schicldrop (personal communication, 1996) in co-operation with the Norwegian Geotechnical Institute. Strictly speaking the centre-of-mass is representative of the frontal part of the slide, projected onto the terrain (the total avalanche centre-of-mass may not even reach the dam). As in the model of Reference NohguchiNohguchi (1989) for centre-of-mass motion on a three-dimensional surface of arbitrary configuration, the equations are derived from classical mechanics, including a resistance force represented by a dynamic drag and a Coulomb friction (as in the Reference VoellmyVoellmy (1955) model). However, a lumped-mass consideration does not comprise any dynamic effects of the avalanche extension. Hence, the model results will anyhow be encumbered with obvious restrictions. For these reasons, it was preferable to perform a simplified-geometry study of the influence of avalanche-impact velocity, terrain inclination, dam configuration, and dam orientation on avalanche-course deflection and run-up height along a deflection dam. An additional advantage of a simplified-geometry study is that the deflecting dam does not have to be superposed on a complex digital terrain.
The simplified geometry consists ofa plane terrain of inclination b and the upper plane wall of the deflecting dam, oriented by its angle relative to the terrain, y, and the angle between the base line of the wall (the x axis) and the terrain contour lines, j, Figure 1.
Fig. 1. Simplified geometry configuration for centre-of-mass model.
The tangential and normal components of the centre-of-mass momentum equations are
and
respectively, where ν is the centre-of-mass velocity at time t, g x = g sin β sin φ g y = −(cos β sin ψ − sin β cos φ cos ψ), and g z = −g (cos β cos ψ + sin β cos φ sin ψ) are componentS of the gravitational force per unit mass, g, in the upper wall plane, g x along and g y normal to the base line respectively, while g z is the component normal to the wall plane. γ is the angle between the centre-of-mass path tangent line and the base line, μ is the dry-friction coefficient, M/D is the mass-to-drag ratio described by Reference Perla, Cheng and McClungPerla and others (1980) and R is the radius of curvature of the centre-of-mass path line on the wall. By means of the kinematic condition ν 2/R = −ν(d/dt) and the transcription (dv/dt = (dv/dγ)(dγ/dt), Equations (1) and (2) can be combined into:
which is solved numerically by a fourth order Runge Kutta procedure. The angle γ is reduced by constant increments dγ throughout the simulations. For each new pair of (v,γ) values, the centre-of-mass is moved a distance ds = ν dt along the upper wall in the direction determined by the value of γ. The time increment dt = ν dγ/(g x Sin γ + g y cos γ) is found by combining Equation (2) and the kinematic condition above.
Also the effects of energy loss due to impact may be investigated. Without any loss, initial values are γ 0 = π/2 − φ and ν 0 = ν T, where ν T is the centre-of-mass velocity immediately before impact. If the centre-of-mass velocity component normal to the upper wall is lost during the impact, initial values are γ 0 = tan −1 (cos ψ/tan φ) and ν 0 = ν T (sin2 φ + cos2 φ cos2 ψ)½
3. One-Dimensional Continuum Model for Three-Dimensional Avalanche Flow
3.1 Physical and Mathematical Model
The avalanche channel path is approximated by a set of volume elements with varying widths, compensating for converging and diverging effects in a real avalanche flow. Furthermore, horizontal centrifugal effects due to the curvature of the horizontal projection of the path are taken into account. A preliminary version of the simulation procedure, in which the centre-line of the avalanche path had to be specified prior to the simulation, has been presented by Irgens (19th ICTAM, Japan, August 1996 to be published). The main feature of the model is the fact that the centre-line of the avalanche is a space curve, which is determined by the terrain and in the present improved version also by the dynamics of the flowing material. This is in contrast to the three-dimensional model of Reference Sassa and BonnardSassa (1988) and Reference Lang and LeoLang and Leo (1994), where the centre-line is a curve in a vertical plane. Tentative two-dimensional models for three-dimensional dense snow avalanche flow are being developed at the Swiss Federal Institute for Snow and Avalanche Research (personal communication, U. Gruber, 1997) and CEMAGREF, France (personal communication, M. Naaim, 1997).
A representation of the three-dimensional avalanche topography is shown in Figure 2. The geometry of an avalanche channel is defined by a preliminary centre-line space curve and terrain profiles in cross sections perpendicular to this line. The centre-line is specified by a selected number of path points, Q 1, . . , Q n, at the bottom of the profiles and defined by Cartesian coordinates (X,Y,Z). The projections of the centre-line in the XY- and XZ- planes are replaced by cubic splines. The centre-line is then subdivided into a chosen number of subsegments by station points. P 1, . . , p m.
Fig. 2. Centre-tine of three-dimensional avalanche. Q1, …, Qn are the path points defining the centre-line. XY is a vertical plane, and XZ is a horizontal plane. The x axis is in the direction of the flow. The 
plane is vertical. ø is the angle of flow inclination with respect to the horizontal plane, αx is the slope of the projection of the avalanche centre-line in the XZ plane with respect to the X axis. g is the gravitational force per unit mass.
The cross-sectional terrain profile at each path point Q i is approximated by a circle of radius R i (X) as shown in Figure 3. The radii R 1 . ., R m of similar profiles at the station points, P 1, . ., P m are found from a cubic spline through points with Cartesian coordinates (X, R) for the path points, Q 1, . ., Q n. By this procedure the real avalanche channel is replaced by a set of elements between the cross-sectional profiles. The avalanche of the flowing material is defined by a subset of these elements filled with snow. The height of the snow is given by h i(X, t) at the station points P i. The circular terrain profile shown in Figure 3 represents the terrain profile both at the path points, Q, . . , Q n and the station points, P 1, . ., P m.
Fig. 3. Circular segment cross-sectional profile at a path point Q ora station point P. 
is the adjusted station point. θ is the angle of inclination of the profile, bXZ is the centrifugal force per unit mass in the XZ plane, h is the height of the flow and is determined In the flow. w is the computational width of the corresponding rectangular profile, and is determined by h and the radius R. w1 is the width of the circular segment.
The profile of the flowing material through the cross section is approximated by a circular segment. Due to centrifugal forces the trace of the free surface in the cross section will be inclined with respect to the horizontal plane. The angle of inclination θ defines the origin of the local coordinate systems xyz and 
as shown in Figure 3. The x axis is tangent to the path curve, and the 
plane is vertical. The origins are taken to be adjusted station points 
with new global co-ordinates (X
i,Y
i, Z
i) for the calculated centre-line of the avalanche. The circular-segment cross section of the flowing material is further replaced by a rectangular cross section of height h and computational width w, and with the same cross-sectional area A as the circular cross section. The assumption of a circular-segment profile implies an inter-dependence between the flow width w
1(X, t) and the flow height h(X, t)�. The slope αX and the curvature k
XZ in the horizontal XZ plane, and the slope ø and the curvature 
in the vertical 
plane are all computed from the coordinates of the stations P
i and based on central-difference formulas.
The projected curved motion of the flowing material in the XZ plane is responsible for a horizontal centrifugal force component b XZ per unit mass.
where ν
a is the average velocity through the cross section of the flow. Due to this centrifugal force the free surface of the flowing material will be inclined with respect to the horizontal 
axis at each station profile. The gravitational force g per unit mass has a driving component g sin ø in the x direction and a component g cos ø in the 
direction. To the latter component we add a centrifugal-force component in a vertical plane
The effective gravitational force, g
e, in the 
plane is the resultant of these forces
This body force defines the angle of inclination θ of the free surface of the avalanche, which is determined from
This angle determines the directions of the coordinate axes y and z in the 
plane. The flow is now considered to be two-dimensional with the velocity field given by the two components ν
X(X, y, t) and ν
Y(X, y, t).
For the sake of simplicity, the complete version of the simulation model presents two special options as alternatives: (1) For highly cohesive material extensional flow with a uniform streamwise velocity ν x = ν o(X, t) is assumed. The constitutive equations contain terms representing active- and passive-pressure contribution. (2) When cohesion may be neglected, shear flow and the no-slip condition ν x = 0 on the bed surface y = 0 are assumed. On the free surface, y = h(X, t), the normal stress must be equal to the atmospheric pressure, which is assumed to be equal to the pore pressure p u, and the shear stress must be zero. The constitutive equations do not produce active- and passive-pressure terms in this case.
3.2 Constitutive Model
A general discussion of the constitutive model may be found in Reference Norem, Irgens and SchieldropNorem and others (1987, Reference Norem, Irgens and Schieldrop1989). For a two-dimensional steady gravity-driven shear flow the equation of motion in the streamwise direction yields the velocity field
where ν 0 is the velocity at the bed surface, ν 1 is the velocity of the free surface and n is a material power-law exponent. The relevant shear and normal stresses for this flow in the x and y directions are
The material properties of the model are specified by the cohesion parameter a, the coefficient of dry friction μ the shear viscosity m, the exponent n, and the two visco-elasti-cities m 1 and m 2, representing the effect of normal stress differences. p e is the effective pressure and ρ is the density of the snow.
3.3 Numerical Simulation
It is assumed that the equation of motion in the y direction may be approximated by an equilibrium equation, which may be integrated to
An expression for the effective pressure P e is obtained by a comparison of Equations (11) and (10). The velocity in the streamwise direction ν x(X,y,t) is assumed to be approximately given by the steady-state function in Equation (8), where the velocity of the free surface ν 1 is now assumed to be a function of X and t. The equation of motion in the streamwise direction and the continuity equation are integrated in the y direction. The assumed velocity profile ν x(X, y, t) is substituted into the integral equations. From previous two-dimensional simulations, it is known that with reasonable parameter values, the results are more or less similar for case 1 and case 2, referred to in section 3.1. For case 2, which is more easily implemented in the numerical model, the following differential equations are obtained
where
and
To determine the motion, the equation of motion (12) together with the continuity equation (13) and the geometrical relation between the flow height and the flow width provide three equations for the three unknowns ν
a, h, and W. The partial differential equations are solved by a finite difference scheme with spatial central differences in the streamwise direction, and by a fourth order Runge-Kutta procedure with respect to time. The finite-difference scheme is Eulerian, requiring the following special procedure to make the avalanche progress along the path: The volume of snow passing through the front section at station 
fills the downstream subsegment, where the subscript f denotes the number of the contemporary front station. When the accumulated volume exceeds the current value of the volume h
f·w
f·Δx
f+1, where Δx
f+1 is the distance between the stations 
and 
, it is assumed that the avalanche front has advanced one subsegment. A similar procedure is applied to the tail of the avalanche.
Initially the snow is assumed to fill a certain number of volume elements and is at rest. The origins of the local coordinate systems xyz are located at the station points P. As the motion starts and new volume elements are filled with snow, the origin 
of the next front station is determined by the angle θ from Equation (7). The result is an avalanche that finds its path according to the terrain and the dynamics of the flow. The width w
1 of the circular segments in Figure 3 are found from a geometrical formula relating h, R, and W
1. The computer program is developed for personal computers.
4. The 1986 Vassdalen Avalanche
4.1 The Observed Avalanche
Extreme snow fall combined with strong winds and a cold period during the first part of the winter followed by temperature variations were the main triggers for the dry-snow avalanche in Vassdalen, Narvik, northern Norway on 5 March 1986. 16 soldiers were killed in the avalanche, which was therefore reported extensively (Reference LiedLied, 1988). The mapping of the avalanche was accomplished shortly after the event, and was based mainly on location of snow deposits and injuries on the birch forest. The deposited snow masses, the severe deflection of the avalanche course, the eye-witnesses accounts, and the limited extension of the injuries all indicate that this was a dense dry-snow avalanche with an insignificant powder-snow cloud. For these reasons, the 1986 Vassdalen avalanche serves as a well-defined full-scale experiment regarding avalanche travel path and extension.
The fracture line of the avalanche was located 475 m a.s.l. The approximately 100 m wide release zone has an upward concave transversal terrain profile and an average inclination of 35.5° above 375 m a.s.l. (Fig. 4). The terrain inclination along the base of the natural deflecting dam is approximately β = 26.5°. The angle between the base line of the wall and the terrain contour lines is estimated to φ = 50°. In the region where the avalanche obtains its maximuni run-up height, the line of steepest descent along the tipper wall of the dam has an inclination of α s = 33°. The deflecting terrain in the lower half of the track caused an absolute maximum run-up height above the base line of the wall of about 25 m measured in a vertical cross section perpendicular to the horizontal projection of the base line of the natural deflecting dam.
Fig. 4. Simulations of the 1986 Vassdalen avalanche by three-dimensional continuum model (dotted line) compared with the observations made (solid line). All angles and lines are horizontal projections. Location of observed maximum run-up height, measured in a vertical cross section perpendicular to the horizontal projection of the base line of the natural deflecting dam, indicated. Contour-line interval is 2 m.
4.2 The Simulated Avalanche
4.2.1 Simulations by Centre-of-Mass Model
Input values β = 26.5°, φ = 50°, α s = 33° and ν T = 28 m s−1 (where the velocity is deduced from the three-dimensional simulations described below), give ψ = 44.6°. Hence the simulation of the centre-of-mass motion along the wall including energy loss due to impact starts with deduced values ν 0 = 25.0 m s−1 and γ 0 = 30.9°. With μ = 0.2 and M/D = 500 m (chosen clearly within the documented possible range of these values), the maximum calculated run-up height of the centre-of-mass above the base line ot the deflecting dam, measured along the terrain in a cross section perpendicular to the base line of the natural deflecting dam. is y max = 16.0m (Fig. 5). This corresponds to a maximum vertical run-up height above the base line of the deflecting dam, measured in a vertical cross section perpendicular to the horizontal projection of the base line of the natural deflecting dam. of r max = 7.7 m, where
and
This is the run-up height estimated from maps and intuitively pointed out in the terrain, and also the most convenient height in dam design.
Fig. 5. Principal sketch of run-up heights on deflecting dam. Carved dotted line indicates outer extension of avalanche flow on the wall of the dam. y is the run-up height measured along the terrain in a cross section perpendicular to the base line of the dam, while r is the vertical run-up height measured in a vertical cross section perpendicular to the horizontal projection of the base line of the deflecting dam.
Without energy loss due to impact, ν 0 = 28 m s−1 and γ 0 = 40.0°. The corresponding maximum run-up height is now y max = 29.5 m or r max = 14.1 m.
4.2.2 Simulations by Three-Dimensional Continuum Model
The input values to this model were chosen to give reasonable agreement between the simulated avalanche and the recorded path and the run-out distance of the real avalanche. The material parameters are: ρ = 300 kg m −3, μ = 0.4, m = 0.00146 m2, m 1 = 0.0144 m2, m 2 = 0.00144 m2 and n = 2.0. Ten path points were chosen to describe the centre-line of the preliminary path. The radii R of the computational cross sections near the initiation of the avalanche and the run-out zone had to be estimated on the basis of the recorded width of the real avalanche. The fracture width and the curvature of the uppermost computational cross section implicitly define the fracture-height area.
5. Calibration and Comparison of Model Results
For the centre-of-mass model a natural first assumption is that the centre-of-mass should reach half the run-up height observed in the terrain. This assumption is valid when the whole avalanche is climbing the wall of the deflecting dam with the masses equally distributed on both sides of the centre-line. Relatively good agreement with these conditions was one of the reasons why the Vassdalen avalanche was chosen for back calculations.
As described in section 4.2.1., the maximum run-up height calculated by the centre-of-mass model without energy loss due to impact is 14.1 in, while the corresponding observed maximum run-up height was 25m (Fig. 4), i.e. close agreement with the assumption above. Including energy loss due to impact reduces the calculated run-up height to 7.7 m. Probably the correct maximum centre-of-mass run-up height is somewhere in between these extreme limits. The avalanche masses were probably distributed more towards the base of the deflecting terrain than indicated by a centre-line, hence reducing the upper value of 14.1 m. On the other hand, due to snow compressibility and smooth terrain conditions, the lower extreme of 7.7 m based on total neglect of the centre-of-mass velocity component normal to the upper wall during impact, is probably too small.
Figure 4 shows the result of the simulation with the three-dimensional continuum model. In comparison with the recorded real avalanche the simulated avalanche is too sensitive to centrifugal effects. This may be due partly to the fact that a very crude approximation to the avalanche chan-nel is used. The simulated avalanche channel in Figure 4 demonstrates the main features of the model: the centre-line is determined by the centrifugal effects, and the width of the avalanche is governed by the topography of the terrain. The maximum flow height is initially 3 m, and increases to 4.5 m at the major bend, and becomes finally 3.5m in the run-out zone.
6. Conclusions
The 1986 Vassdalen dense dry-snow avalanche is back calculated both by a centre-of-mass model applying the Voellmy-Perla equation, and a quasi-three-dimensional continuum model applying one-dimensional depth inte-grated equations of mass and linear momentum. The centre-of-mass path line calculated by the centre-of-mass model is in close agreement with the observations when energy loss due to impact is neglected. The model is applicable for studying the influence of terrain inclinations, impact velocity and dam configurations on the avalanche centre-of-mass path line. However, the model results are still encumbered with uncertainties, even though the parameters describing the features above are included. Whether the whole avalanche is elimhing the wall depends on the same features. If a considerable part of the avalanche is not climbing the wall, the simulations should be made by consideration of a characteristic width of the slide representing the climbing fraction of the slide only. However, it is cumbersome to detinc an objective criteria for calculation of such a width. Besides the interaction with the masses not climbing the wall, being deflected or not, should also be taken into account. Hence, for future design of deflecting dams, there is hardly any alternative to (quasi) three-dimensional simulations.
The three-dimensional continuum model has in its present version several deficiencies. It may be difficult to model the terrain satisfactorily using the proposed avalanche channel. In future work a better approximation to the cross-sectional profiles will be attempted. From the present simulation it is clear that the model is too sensitive to centrifugal effects. The cubic-spline approximation of the centreline has unwanted consequences, as it sometimes introduce ripples in the centre-line which give rise to unrealistic curvatures and centrifugal effects.
Acknowledgements
Thanks are due to colleagues at the Norwegian Geoteclmical Institute for their technical advice and discussions.
