## Abstract

The Savage–Hutter (SH) avalanche model is a depth-averaged dynamical model of a fluid-like continuum implementing the following simplifying assumptions: (i) density preserving, (ii) shallowness of the avalanche piles and small topographic curvatures, (iii) Coulomb-type sliding with bed friction angle *δ* and (iv) Mohr–Coulomb behaviour in the interior with internal angle of friction *ϕ*≥*δ* and an *ad hoc* assumption reducing the number of Mohr's circles in three-dimensional stress states to one. We scrutinize the available literature on information regarding these assumptions and thus delineate the ranges of validity of the proposed model equations. The discussion is limited to relatively large snow avalanches with negligible powder snow component and laboratory sand avalanches starting on steep slopes. The conclusion of the analysis is that the SH model is a valid model for sand avalanches, but its Mohr–Coulomb sliding law may have to be complemented for snow avalanches by a second velocity-dependent contribution. For very small snow avalanches and for laboratory avalanches starting on moderately steep and bumpy slopes it may not be adequate.

## 1. Introduction

The Savage–Hutter (SH) avalanche model (Savage & Hutter 1989, 1991) and its extensions (Hutter & Nohguchi 1990; Hutter & Koch 1991; Greve & Hutter 1993; Hutter & Greve 1993; Hutter *et al*. 1993, 1995; Greve *et al*. 1994; Koch *et al*. 1994; Gray *et al*. 1999; Tai *et al*. 1999; Tai 2000; Pudasaini & Hutter 2003; Pudasaini *et al*. 2003*a*; Wang *et al*. 2004), henceforth also called SH model, is a dynamical fluid-like model, which consists of hyperbolic partial differential equations for the distribution of the depth and the (two) topography-parallel, depth-averaged velocity components in an avalanching mass of cohesionless granules (e.g. sand, grains, rocks and snow). It is designed to predict the motion and deformation from initiation to run-out along a concomitantly determined avalanche track along a prescribed topography. In the past, it has been used to describe flows in straight and curved chutes (Hutter & Nohguchi 1990; Hutter & Koch 1991; Greve & Hutter 1993; Hutter *et al*. 1993, 1995) in channels with plane and parabolic cross-sections and simply curved thalwegs (Voellmy 1955; Greve *et al*. 1994; Koch *et al*. 1994; Hutter 1996; Gray *et al*. 1999; Tai 2000; Tai *et al*. 2001, 2002; Koschdon & Schäfer 2003), but has been extended to flows in corries having arbitrarily curved and twisted thalwegs and arbitrary topographies (Pudasaini 2003; Pudasaini & Hutter 2003; Pudasaini *et al*. 2003*a*,*b*). The basic simplifying assumptions in the various models are mathematically not exactly the same; however, physically they are identical, namely consisting of

the assumption of density preserving (incompressibility);

the assumption of shallowness of the avalanche piles, and of small topographic curvatures;

the assumptions of Coulomb-type sliding with bed friction angle

*δ*;nearly uniform velocity profile through the avalanche depth and

Mohr–Coulomb frictional behaviour in the interior with internal angle of friction

*ϕ*≥*δ*and an*ad hoc*assumption, reducing the number of Mohr's circles in three-dimensional stress states from three to one.

All these assumptions can be justified, but they limit the applicability of the model equations. In the ensuing study we scrutinize the literature and provide documentation under what conditions the SH avalanche equations are likely to be valid. Specifically, it is shown both for laboratory avalanches and for flow avalanches of snow that the above stated five simplifying assumptions are validating the model equations to the extent that they are able to reproduce laboratory and field experiments sufficiently accurately. This is not to say that one or the other assumption (i)–(v) would always be fulfilled, but that its violation may be minor or, if it is large, of short duration or length as compared with the global scales, that it will hardly be identifiable in laboratory or large scale observations. This neither means that there would not be a number of situations in which the model equations fail as an adequate predictive tool. This happens for instance for very small snow avalanches (slab avalanches) or for avalanches starting on slopes with a small inclination angle, or for the motion over very rough beds or abrupt topographic steps. All these cases are mainly short-lived and thus do not reach a significant level. Avalanches of this sort move on inclines provided their height is sufficiently large, and they often come to a premature stop. They have been studied in the laboratory by Pouliquen (1999) and Louge (2003) and by discrete particle methods by Silbert *et al*. (2003).

In §2 the ranges of validity of the aforementioned five assumptions will be identified; this section deals with laboratory experiments and field observations, and the assumptions (i)–(v) will be justified. The performance of the model equations will be discussed in §3; this means that model output will be compared with the corresponding findings from experiments, but it equally also relates to the computational performance of numerical codes in attempts to integrate the free boundary value equations. Finally, §4 brings an outlook of unsolved problems to be handled in the future.

## 2. The simplifying assumptions and their ranges of validity

### (a) Density preserving

There are hardly any measurements available on volume changes in rapid shear flows of snow or granular materials. However, Hutter & Koch (1991) reported measurements of areal changes of granular flows in an exponentially curved chute of 10 cm width by taking fast speed photographs from the side and measuring the areas of the individual pile shapes. These measurements revealed volume changes in sand-avalanche-chute flows initially of 15–20%, i.e. between the pile volume at rest and that of the first photograph in motion, less than approximately 4% during motion, i.e. between any two photographs when the avalanche was moving and inconclusive information during the settling process. This is consistent with what one would expect. A considerable expansion when the granular mass is set into shearing, practically remaining density preserving in rapid flow and experiencing a compaction when the avalanche stops suddenly or gradually.

Field observations on density variations in artificially released flow avalanches of snow were made by Gubler *et al*. (1986) and Gubler (1993). These authors use microwave FMCW techniques. Their measurements provide limited information about the density variation with depth of the snow in motion, but the results are not conclusive enough to infer that density variations would play a dynamic role in the avalanche motion.

Qualitatively, these results are plausible. Snow avalanches in the flowing regime are relatively dense, bouncing is seldom seen, much less than in laboratory sand flows, and when so, only in dry snow avalanches. So, the mean particle distances are probably not much larger than the particle diameters. Moreover, shearing seems to be weak and restricted to the thin fluidized layer at the base, which is most probably small (see also velocity measurements later on). So, the depth changes due to the associated dilatations are small. It follows that the conditions of a Boussinesq medium are satisfied, but since no traces of buoyancy effects have been seen, *the assumption of density preserving is an appropriate approximation*.

The above statements do not conflict with the fact that many real snow avalanches are of mixed type containing disjoint flow and powder avalanche regions as, for example, described by Issler (2003). The flow avalanche part remains particularly density preserving.

### (b) Shallowness of the avalanche piles and small topographic curvatures

Observations of dense flow avalanches show that the moving snow masses are thin, long and wide and have an aspect ratiowhich is very small and is of the order 10^{−3}–10^{−4}. Avalanches with large aspect ratios (*ϵ*∼1) are generally small, travel short distances and practically never cause damage nor constitute a danger. Therefore, they are less important. So, approximating non-dimensionalized equations to terms linear in *ϵ* may be sufficient. This also requires that local radii [*R*] of curvature are large, such that , 0<*α*<1, that the coefficient of bed friction angle is equally , 0<*β*<1, and that dimensionless shearing in the bulk velocity profile is , with *γ*=min{*α*, *β*}. These assumptions are stringent restrictions. Topographic variations must be relatively smooth, basal friction ought to be limited (namely *δ*≤*ϕ*), so bumpy bottom boundaries are excluded. The order of magnitudes as expressed above are mathematical requirements necessary to derive equations (3.1)–(3.3), (below), by a rational scaling procedure. They make the equations appear to be valid only under rather limited restrictions. However, their applicability may be (and is often) stretched beyond these formal limits. For instance, the transition from a steep corrie to a horizontal run-out may be accomplished by a curved region where [*L*]/[*R*]=(1), and computational results may still agree well with experimental findings. Or the sidewise curvature in a corrie may neither be of , *γ*>0 with maintained satisfactory agreement between computational and experimental findings. We have encountered many such situations. In reality, topographic regions with large curvature are often small and then only have a limited effect on the overall flow of the avalanche. When, on the other hand, avalanches shoot over a topographic step and become airborne and go through a ballistic motion, they often turn into a powder snow avalanche after impinging on the ground. For laboratory experiments on this ‘ski-jump’ effect, see Hakonardottir *et al*. (2003). In these situations the model, of course, ceases to be valid. On the other hand, today's geographical information systems (GIS) are usually based on a 25×25 m^{2} grid, while a 5×5 m^{2} grid would be more appropriate; therefore they do not allow for the topography to be sufficiently accurately resolved to account for local bumps and troughs.

Thus, the shallowness assumption is reasonably satisfied and may only occasionally be violated. In circumstances when it is only locally violated the simplified model equations are likely to generate reasonable solutions, see Hutter *et al*. (1995).

### (c) Bed friction only Coulomb-type

All classical avalanche models write the basal shear stress as a combination of a normal-stress-dependent Coulomb term and velocity dependent terms.1in which *δ* is the bed-friction angle, *p*_{⊥} the pressure normal to the base, *c* the velocity-dependent drag coefficient, generally treated as constants but varying from avalanche to avalanche, and * v* is the velocity tangential to the bed, see Salm (1966), Perla

*et al*. (1980), Salm & Gubler (1985) and Voellmy (1955). Hutter and others almost exclusively set

*c*=0 and state that granular avalanches in the laboratory can be well reproduced by the Coulomb term alone (Gubler

*et al*. 1986; Hutter & Nohguchi 1990; Hutter & Koch 1991; Greve & Hutter 1993; Hutter & Greve 1993; Hutter

*et al*. 1993, 1995; Greve

*et al*. 1994; Hutter 1996; Gray

*et al*. 1999), sometimes by accounting for a reduction in the value of

*δ*in the rear portion of the avalanche (Gray

*et al*. 1999; Wieland

*et al*. 1999), because abraded, segregated, and deposited fines reduce the bed friction angle. On the other hand, identifications of

*δ*and

*c*from observations and back calculation of field avalanches, come up with both

*δ*and

*c*(e.g. Gubler

*et al*. 1986; Gubler 1993; Zwinger 2000; Ancey 2001; Zwinger

*et al*. 2003; Ancey & Meunier 2004). For the run-out process, shortly before the avalanche comes to an immediate stop, there were even suggestions to parameterize

*c*as

*c*∝|

*|*

**v**^{−3}, see Schaerer (1974).

It is well known that with a Coulomb sliding law alone and a constant bed friction angle, a steady flow2 of an avalanche down an inclined plane cannot be obtained; the avalanche is continuously accelerating. (See, however, the work of Savage & Nohguchi 1988 for an exponentially curved bed and when *δ* is not a constant; and work by Hutter & Nohguchi 1990; Hutter & Greve 1993; Hutter 1996). Based on this, it is sometimes argued that, if the steady state motion exists, an additional friction process, e.g. a viscous component, is needed. In the field, this question can probably never be settled. Verification is not possible, since avalanche tracks are not plane, not uniformly rough and generally too short in the sense of producing steady state behaviour. At least as far as field avalanches are concerned, an additional friction mechanism only needs to be added parametrically when large discrepancies with observations arise. From this it follows that existence of a steady asymptotic flow in natural events is rather a belief of the scientist demanding it than an established fact. Pouliquen & Forterre (2002) find steady flows in their laboratory experiments on bumpy inclines with inclinations smaller that 30°. Their chute has a length of 2 m, and glass beads were used. Our own experiments done with sand of 0.5 mm nominal diameters were conducted on Plexiglas chutes of 0.65, 1.85 and 2.85 m length, 0.275 m with an outflow orifice height of 25 mm and inclination angles between 30° and 45°. Our experiments showed no conclusive information as to whether steady state flows will asymptotically be reached3 (see Eckart *et al*. 2003). These results led us to conclude that in the range of parameters (*δ*, *ϕ*, inclination) similar to those of naturally released avalanches, extension of a purely Coulomb-type friction law to include a velocity-dependent drag does not seem to be compelling.

In a recent paper, in which Ancey & Meunier (2004) performed a back analysis of 15 documented avalanche events, the bulk frictional force, experienced by an avalanche was computed. They state that “three types of rheological behaviour were identified: (i) the inertial regime, where the frictional force drops to zero, (ii) the Coulombic frictional force, where the force is fairly independent of the avalanche velocity, and (iii) the velocity-dependent regime, where the force exhibits a complicated (nonlinear and hysteretic) dependence on velocity. During its course, an avalanche can experience one or several regimes. Interestingly, the Coulomb model can provide predictions of the velocity and run-out distance in good agreement with field data for most events, even though for some path sections, the bulk frictional force departs from the Coulomb model”.

Results obtained by Zwinger (2000) point in a slightly different direction. He models mixed avalanches consisting of powder and flow avalanche components, in which the SH flow avalanche component is used in combination with a turbulent, particle-laden powder avalanche with a coupling through a saltation layer. In post computations of the Madlein avalanche of 1984 in the community Ischgl (Austria) he found that reasonable agreement with the observed deposition could only be obtained when the sliding law was Coulomb-type at low velocities and viscous-type at large velocities. No extensive parameter study seems to have been made by Zwinger, so that better agreement would also be possible when a careful parameter identification would be made. Besides this, the Madlein avalanche is mixed and thus of a different class.

In situations of water-saturated debris flow, the Coulomb friction law with vanishing velocity dependent component is still an acceptable friction law, but the bed friction angle is made to depend on the pore water pressure, reducing the friction angle with increasing pressure (Iverson 1997; Iverson & Denlinger 2001; Wang *et al*. 2004). A similar dependence has also been observed in dry avalanches on small slopes (Pouliquen 1999), but the bed friction angle is made depth-dependent in that paper rather than pressure-dependent. They only correspond when the hydrostatic pressure assumption holds true. Such dependence is also employed by Mangeney *et al*. (2003).

### (d) Shearing is unimportant except at the very base

Typical internal, *ϕ*, and bed, *δ*, friction angles for real snow avalanches and for laboratory avalanches of quartz, sand, marmor chips and vestolen (plastic beads) moving on Plexiglas chutes or chutes coated with drawing paper and sandpaper (SIA120), respectively, are (Gubler *et al*. 1986; Hutter 1996; Iverson 1997)This implies that the basal surfaces are smoother than the material under motion, and the smallness of the bumpiness makes it unlikely that the bed transmits considerable shearing into the moving pile. Note also that the ranges of *δ* and *ϕ* stated above cover similar intervals, making inferences from laboratory results to large scale snow avalanches possible, since these are the only two phenomenological parameters in the SH model with which the material response of a granular avalanche is described.

Data on measured velocity profiles are scarce for both field and laboratory experiments. Velocity profiles were measured at selected points in artificially released avalanches by Dent *et al*. (1998) for an avalanche in Montana. From a shelter behind a rock nose in the avalanche track the passing avalanche was observed through a window from bottom to top and pure plug flow was observed. Gubler *et al*. (1986) and Gubler (1993) used radar Doppler measurements to determine the depth variation of the downhill velocity in an artificially released flow avalanche at the Lukmanier pass in Switzerland; shearing was observed, but it was small, and a statement that it needed to be accounted for was not conclusive. Kern *et al*. (2003) performed chute experiments at Weißfluhjoch, Switzerland, 2680 m.a.s.l. in a 38 m long rectangular channel with snow and observed that their artificial snow avalanches suffered shearing within the lower most layer, approximately covering 20–30% of the avalanche depth, above which plug flow and below which sliding was observed. Kern *et al*. roughened the bed artificially so that the bed friction angle may have been large, above say 30°, but they did not provide measured values. So, the experiments probably constitute a situation not typical for flow avalanches. In spite of this, we will demonstrate below that up to 30% of shearing velocity and 70% sliding velocity does not (at all) invalidate the SH equations.

Eckart *et al*. (2003) conducted laboratory experiments with sand on Plexiglas and used particle image velocimetry (PIV) techniques to view the moving mass from above, below, and from the side, so measuring the surface and basal velocity as well as its profile at the sidewall of the uniform steady flow. Chutes 65, 185 and 285 cm long and 27.5 cm wide, inclined at 30°, 35°, 40° and 45° were used. Shearing was observed immediately at the outlet of the material, but disappeared after a distance of 50 cm. Below this region of flow establishment, basal and surface velocities were equal with a relative deviation of less than 3.5%. The shearing in the flow establishment region depends on the inclination angle of the chute and the constructive details of the outlet mechanism. In this case the granular material is released from a rectangular tank by quickly lifting a gate, thus freeing the particles along a line perpendicular to the base. Because of the basal friction, the bottom particles are held back, while the top particles are free to move. This automatically introduces a shearing that is weakened as the material moves down the slope as a free-surface flow.

Pudasaini (2003) and Pudasaini *et al*. (submitted) conducted similar experiments with a laboratory chute consisting of an inclined plane that merges via a cylindrical segment into a horizontal plane, but this time the flow is completely three-dimensional, unconfined, unsteady and non-uniform. PIV measurements were made from above and below in the steep portions of the chute close to the transition region with similar results as reported above: the mean deviation of velocities at the base and the free surface is no more than 3%.

### (e) Mohr–Coulomb behaviour

Because in a depth-integrated model only the stress states at the free surface and at the base enter the description, a detailed model for the mechanical constitutive behaviour of the stresses is not needed. It suffices to formulate the constitutive relations for the stresses at the boundary points and to interpolate. Nevertheless, the granular structure of the material is still thought to be important insofar as active and passive stress states are distinguished according to whether the flow is extending or compressing in the direction of the stress considered.

The following approximate description of the Mohr–Coulomb behaviour is probably the most critical assumption of the model. It is based on the observation of real flow avalanches that the motion is primarily unidirectional, and transverse shearing is small. Consider a local Cartesian coordinate system at a basal point, in which *x* is in the downhill direction, *y* is transverse and *z* is orthogonal to the two, figure 1*a*. Moreover, consider an infinitesimal cube (its lower face lies in the sliding surface) with the stress vectors as indicated on the visible faces; *p*_{zz} is the overburden pressure and *τ*_{xz} the shear traction, at the base as exerted on the sliding surface; *p*_{yy} is the pressure exerted on faces of the cube normal to the *y*-direction. It is assumed that the shear stresses *τ*_{xy} and *τ*_{yz} are negligibly small, so that *p*_{yy} is very close to a principal stress, which will be assumed; *τ*_{xy} is small since transverse shearing is supposed to be small—the motion is practically all downhill—and *τ*_{yz} cannot be large, because, by construction, *x* is in the direction of the local velocity. Given the pair (*b* for base), it must lie in the stress plane (*p*, *τ*) on the line through the origin, forming the angle *δ* with the *p*-axis. All stress states on planes at the base perpendicular to the *y*-axis then lie on the active (dashed) or passive (solid) circles through that are tangent to the wedge with vertex angle 2*ϕ*, and the points are obtained by a rotation of 180° on these circles. Of the two possibilities ‘act’ (‘pass’) is selected if the normal strain rate in the *x*-direction, ∂*u*/∂*x* is positive, i.e. extensional (negative, i.e. compacting). The principal stresses and the direction of the elements at which they apply can also be calculated; they are given as indicated by the solid circles in the Mohr diagram in figure 1*b*. Given , all stress states can be expressed in terms of these as well as *ϕ* and *δ*.

There remains the determination of . To this end it will now be assumed that of the three principal stresses two will agree with one another. This *ad hoc* assumption makes the determination of unique and sets it equal to one of the four principal stresses on the two circles shown in figure 1*b*. The sign of the longitudinal normal strain rate, ∂*u*/∂*x*, selects the active or passive circle, and the sign of the transversal normal strain rate, ∂*v*/∂*y*, will determine the smaller (∂*v*/∂*y*>0) or the larger (∂*v*/∂*y*<0) of the two values on the respective circle. Thus,(2.1)where , are the earth pressure coefficients that are functions of *ϕ* and *δ*. This is the least rational assumption of all, and it is only justified by the results it produces. It should, however, also be mentioned that the above *ad hoc* assumption destroys the rotational invariance of the equations about the *z*-direction perpendicular to the tangential plane of the reference surface. The reason for this lies in the omission of the stresses *τ*_{yz} and the preference of the *x*-velocities in the above construction of (2.1). Iverson & Denlinger (2001) in their model do not make this assumption and preserve rotational invariance.

Two additional remarks are in order. First, since the avalanche motion is predominantly downhill with large velocity components in this direction and small velocities sidewise, *xy*-shearing must be small, as will (less convincingly) also be the *yz*-shearing. So, on an element surface perpendicular to the *y*-axis the traction (*τ*_{xy}, −*p*_{yy}, *τ*_{yz}) may be approximated by (0, −*p*_{yy}, 0), making −*p*_{yy} essentially a principal stress. This assumption is probably only violated when the avalanche experiences large *xy*-shearing. Moreover, Vollmöller (personal communication) reports computations performed by using our closure assumption and that of Iverson & Denlinger (2001), respectively, for a flow down an incline into a horizontal plane, which show graphically nearly invisible differences of the avalanche topography.

Second, the earth pressure coefficient arising on the right-hand side of (2.1) suffers jump discontinuities when active-passive flow transitions occur. These jumps are computationally accounted for. However, in one-dimensional flow the value of the earth pressure coefficient at ∂*u*/∂*x*=0 can be determined experimentally. This has been done by Tai *et al*. (2000) in rotating drum experiments, including a proposal of a continuous transition of the earth pressure coefficient as a function of ∂*u*/∂*x* from active to passive flow states such that for ∂*u*/∂*x*→±∞ the constant values *K*_{act} and *K*_{pass} are reached.

At the free surface all stress components must be zero in a cohesionless Mohr–Coulomb material if the traction from above is set to zero. Moreover, since *p*_{zz} and *τ*_{xz} are linearly distributed, it is also reasonable to assume a linear dependence for *p*_{xx} and *p*_{yy}, connecting at the base with (0,0) at the free surface. This is compelling if the earth pressure coefficient does not show a *z*-dependence.

## 3. Performance of the model

### (a) Governing equations (Pudasaini & Hutter 2003)

The SH equations were derived in various different coordinate systems. In appropriately chosen orthogonal curvilinear coordinates they take the following conservative form:

Mass balance(3.1)

Momentum balances(3.2)(3.3)in which(3.4)Here, *x*, *y* are curvilinear coordinates in the directions along and perpendicular to the curved and twisted master curve; *h*, *hu*, *hv* are the avalanche depth and the specific momenta in the *x*- and *y*-directions, respectively; *g*_{x}, *g*_{y}, *g*_{z} define the gravity components in the three orthogonal directions *x*, *y*, *z*; *λκ* is the local radius of curvature of the master curve and *β*_{x}, *β*_{y} incorporate the earth pressure coefficients in the *x*- and *y*-directions, respectively. Moreover, *b*(*x*, *y*) defines the basal surface, i.e. the deviation of the basal topography from the reference surface *z*=0, while *φ* gives the accumulation of the torsion of the master curve from an initial position and *ϕ*_{0} is an arbitrary constant (see Pudasaini & Hutter 2003). The parameters *α*_{xx}, *α*_{xy}, *α*_{yy}, are the so-called momentum correction factors for the velocity distributions with depth and defined aswith , and and being the local velocity components before the depth averaging. The correction factors *α*_{xx}, *α*_{xy}, *α*_{yy}, convert the fluxes of momenta in terms of , to the average values *u*, *v*. They are larger than unity when the depth profile of and deviate from uniform distributions, and they in general differ from one another. For simplicity, one often sets *α*_{xx}=*α*_{xy}=*α*_{yy}=1; see, however, later computational results.

Equations (3.1)–(3.4) form a hyperbolic system of partial differential equations with coefficients that may have jump discontinuities and can be written in standard mathematical form:(3.5)where * u*,

*,*

**f***, and*

**g***are vector-valued functions. For their solution the input quantities are the topography (*

**s***b*(

*x*,

*y*), master curve) the internal and bed friction angles,

*ϕ*,

*δ*, and the initial values for

*h*,

*u*,

*v*; the output is defined by the three conservative functions

*h*(

*x*,

*y*,

*t*),

*hu*(

*x*,

*y*,

*t*) and

*hv*(

*x*,

*y*,

*t*) or equivalently

*h*(

*x*,

*y*,

*t*),

*u*(

*x*,

*y*,

*t*) and

*v*(

*x*,

*y*,

*t*). Discontinuities may arise in these variables as well as their gradients.

### (b) Numerical schemes and performances

In order to test the above model equations against avalanche events either in nature or in the laboratory, a numerical integration scheme must be constructed. Early attempts (Savage & Hutter 1989, 1991; Hutter & Koch 1991; Greve & Hutter 1993; Hutter & Greve 1993; Greve *et al*. 1994; Koch *et al*. 1994; Hutter *et al*. 1995; Gray *et al*. 1999; Wieland *et al*. 1999) used Lagrangian finite-difference schemes with central difference approximation and leap frog temporal integration steps. These approaches made addition of explicit numerical diffusion necessary, but it was held minimal and made operative where gradients of the avalanche thickness and velocities became large. These methods were not able to capture shocks and may have smoothed these out. Shock capturing finite-difference techniques were introduced later. Tai (2000) and Tai *et al*. (2002) used them in a two-dimensional Eulerian shock capturing scheme and a one-dimensional front-tracking method, respectively. Wang *et al*. (2004) employed a high resolution approach, namely the non-oscillatory central (NOC) scheme, and compared different cell reconstruction techniques—four different second-order total variation diminishing (TVD) limiters and a third-order essentially non-oscillatory (ENO) cell reconstruction scheme. Of the numerical methods under consideration, the NOC scheme with the Minmod TVD limiter showed the best performance in chute flows down a channel merging into a horizontal plane.4 These Eulerian schemes, while superior to the above-mentioned central difference schemes, are comparably accurate to other shock capturing methods. Similarly, Denlinger & Iverson (2004) and Iverson *et al*. (2004) use shock capturing techniques to solve their extended model equations on a rectangular Cartesian coordinate system. Koschdon & Schäfer (2003) use an arbitrary Lagrangian–Eulerian finite-volume method, where unstructured boundary-fitted moving grids are employed to follow the free boundary. The underlying flow solver consists of a Godunov-type approach in the space-time domain, and the fluxes are calculated using Riemann solvers. Vollmöller (2004), on the other hand, uses a wave propagation method in the context of unstructured finite volumes and on the basis of Godunov-type schemes with spatially discretized flux functions.

Comparisons of numerical solutions with laboratory chute flows have been conducted in Hutter & Koch (1991), Savage & Hutter (1991), Greve & Hutter (1993), Greve *et al*. (1994), Koch *et al*. (1994), Hutter *et al*. (1995), Hutter (1996), Gray *et al*. (1999), Tai *et al*. (1999, 2001, 2002), Wieland *et al*. (1999) and Tai (2000), and will not be repeated here. However, the results of the motion of a finite mass of sand down a plane Plexiglas chute merging into a horizontal run-out are summarized in figure 2. It compares the geometry of the deposit as obtained experimentally and computationally by Pudasaini (2003) using the integration technique in Tai (2000) and Wang *et al*. (2004). The chute (4 m long and 1.6 m wide) is made of Plexiglas, see inset in figure 2. The upper plane part is inclined at 45°; it merges into a short cylindrical element with rear boundary at 1.56 m and front boundary at 1.93 m from the top of the chute and is followed by the horizontal run-out plane. The sand is initially held in a scull cap that is quickly tilted upward to release the material, which here is quartz of 4 mm nominal diameter, and internal and bed friction angles are *ϕ*=35° and *δ*=23°, respectively. The cap is taken from a hemisphere of radius 0.195 m and has a base radius of 0.17 m, for a detailed description see Pudasaini *et al*. (submitted).

One of the most important aspects in avalanche dynamics is the determination of the run-out area and the height profile of the avalanche in its deposit. Figure 2 shows thickness contours of the deposit at rest of the laboratory avalanche just described, in panel (*a*) as obtained by computation, in panel (*b*) as determined from the laboratory avalanche by determining the depths at about 200 points with the aid of a penetrometer. The level lines are done by using the raw data, not using any smoothing. This explains the relatively rough shapes of the experimental level lines. Considering that the deep blue outermost ring in panel (*b*) has almost zero depth, the two deposits are very similar; the computed pile is slightly less wide, and its rear margin has moved somewhat farther than its experimental counterpart, but the front margin is in both graphs at the same position. Thus, the lateral and the longitudinal run-out distances, the over all run-out zone as well as the height profiles are well predicted by the theory.

The true test of a shock capturing integration method is a situation where abrupt changes of flow heights and/or velocities occur. Experiments were conducted in the Laboratory of the Department of Mechanics, Darmstadt University of Technology for uniform granular flows down an inclined plane that is diverted by a pyramid, a circular cylinder or a wall. Shock waves, dead zones and/or particle free regions were formed that were computationally reproduced by Gray *et al*. (2003) using the NOC scheme introduced by Nessyahu & Tadmor (1990) and extended to multidimensions by Jiang & Tadmor (1997) and Lie & Noelle (2003). The results of Gray *et al*. have been reproduced by using the integration routine of Wang *et al*. (2004).

A robust test of an integration code for the SH equations is the flow of a finite-granular mass down an inclined plane impinging on obstructions. Consider a hemi-ellipsoidal shell holding the material together, which is suddenly released so that the bulk material commences to slide down an inclined flat plane at 40° into a horizontal run-out plane connected by a smooth transition. The computational domain is the rectangle *x*∈[0, 30] (along the downslope direction) and *y*∈[−10, 10] (in the cross-slope direction) in dimensionless length units, where the inclined section lies in the interval *x*∈[0, 20] and the horizontal region lies where *x*≥24 with a smooth change in the topography in the transition zone, *x*∈[20, 24]. As an example, an obstructing equilateral tetrahedron with dimensionless bottom-side length of *l*=4 and height *H*=5 placed symmetrically in the inclined portion of the chute, one wedge pointing uphill and the side opposite to it arranged horizontally at *x*=17. A direct impression of the granular flow can be obtained from the corresponding three-dimensional evolutions of the flow geometry, which are shown in figure 3*a*. When the granular material hits against the obstacle, the flow is diverted. Immediately below the obstacle, and on either side of it, two knolls form and extend downwards; these can be explained as weak shock waves caused by the obstacle. Furthermore, behind the wedge the reduction of the in-flowing mass causes a dent to form. Basically, the material accelerates until it reaches the horizontal run-out zone where the front comes to rest. When this occurs, the tail is still in the accelerating phase. In this stage the granular body contracts, generally accompanied with a shock wave comprising steep surface gradients, where the velocity changes from a supercritical to a subcritical stage. This occurs just before the end of the transition zone. With the approaching mass from the tail, the shock wave propagates backwards (from *t*=14 to *t*=21) until it reaches the upslope margin of the deposit (*t*=21). At *t*=21 the granular mass comes to rest. The final deposition (for *t*>21) has two local maxima with (nearly) the same height due to the symmetry of the flow with a valley in-between. When the obstacle is sufficiently high, as for the case displayed in figure 3*a*, all the granular material approaching the obstacle is diverted to its sides and flows around the obstacle downwards. The shocks that are formed on either side of the obstacle rapidly spread in the transverse direction and merge into two expansion fans, between which a so-called granular vacuum is formed, which can be seen as the protected zone in a practical application of the avalanche protection.

The behaviour of the granular flows past a blunt-nosed cuboid wall obstacle with dimensionless length *l*=4, thickness *b*=0.8, and height *H*=5, located transversely in the middle of the inclined plane and perpendicularly to it with the centre in (*x*, *y*)=(15, 0), is illustrated in figure 3*b*. The results are generally fairly similar to those past a wedge as shown above in figure 3*a*. For the blunt-nosed obstacle, the granular material is partly deposited in front of the obstacle. While this deposition in front of the blunt-nosed obstacle is built, an upstream propagating shock wave is formed as the velocity changes from supercritical to subcritical (from *t*=7 to *t*=14 in figure 3*b*). This phenomenon is fairly similar to the deposition process of the granular material on the horizontal run-out zone as it can be observed from *t*=14 to *t*=21 in figure 3*a*,*b*. The topography of the deposited mass depends strongly on the flow around and over the obstruction. For a sufficiently high wall, as displayed here in figure 3*b* with *H*=5, the flow is around the sides of the wall and the deposited mass consists of two separate heaps with lobes connected to the sides of the wall at much higher elevation, which are unrealistic and not reproduced when the integration is performed for the flow over and around a tetrahedral wedge (figure 3*a*). Admittedly, when impinging a wall, the SH equations cannot be a valid set of equations to properly predict the flow in the immediate vicinity of the wall. The formation of the dead zone behind the wall, however, changes the geometry of the flow in the right direction. In any case, a numerical scheme that does not show these lobes would certainly be more trustworthy; in spite of this, corroboration of the solution by experiment is still necessary.

It is apparent that numerical integration of the strongly convective SH equations requires a robust numerical programme that must be tested in benchmark problems involving shocks, dead zones, and granular free regions. This is the reason why exact solutions to special problems are useful. Of course, less sophisticated discretizations may prove adequate when shocks or steep gradients of the fields do not arise. In any case, only a trustworthy numerical programme can be used when computational results are compared with experiments.

### (c) Comparison with observations

It is not necessary here to repeat the many comparisons that were conducted with numerical output and experimental findings. In the laboratory the moving granular masses were photographed with fast speed cameras operating between 4 and 15 frames per second. This simple technique allowed comparison of the avalanche circumference as the moving mass evolved through time from initiation to run-out. With less accuracy the photos also permitted estimation of the velocity distribution. However, the latter became only experimentally available, once the particle image velocimetry technique was introduced. Excellent comparisons were obtained in confined flows down chutes consisting of plane segments (Hutter *et al*. 1995), an exponentially (Hutter & Koch 1991) or concave-convexly (Greve & Hutter 1993) curved profile in which a single hump was split into two separate depositions above and below the topographic bump. Further comparisons were performed for transversely unconfined flows down an inclined plane merging into a horizontal plane (Greve *et al*. 1994; Koch *et al*. 1994; Vollmöller 2004; figure 2 shows an example of this sort), or a parabolic channel merging into a horizontal plane with curved thalweg (Gray *et al*. 1999; Wieland *et al*. 1999). Steep shock-like depth changes arose in the transition zones from supercritical (dilating) to subcritical (compacting) flow conditions that could well be reproduced by the shock-capturing integration techniques, but shock forming flows down inclined planes, diverted by wedges, walls and cylinders (Tai *et al*. 1999) could only adequately be reproduced by the shock capturing techniques (Gray *et al*. 2003; Wang *et al*. 2004). While these results are pleasing and provide support for the SH model, there is still the need for comparison of experiments and theory when dead zones are formed. Such studies are presently underway.

Whereas laboratory experiments can be studied under isolated, well controlled conditions, this is not so for natural avalanche events, even if these are artificially released. Comparison of the SH equations with such events are very rare and generally less convincing than for granular avalanches in the laboratory. Problems arise with the estimation of *δ* and *ϕ*, the discretization of the topography, the estimation of the moving mass, the neglection of entraining mass from below and possible deposition, which may be present in some instances to a large extent. The only example we know is reported in Zwinger (2000) and Zwinger *et al*. (2003).

In this regard it seems important that input parameters of the SH model are randomly varied within practically reasonable ranges and that probability distributions are determined for the output quantities.

In summary, it appears that computational routines are now available that reliably allow construction of trustworthy solutions of the SH equations. The computational output, when compared with laboratory experiments proved the equations to be a reliable model for most tested configurations. Further scrutiny does, however, seem to be necessary for geometrically more complex situations for which the validity of the model seems to be limited. Predictions or post-computations of real avalanches should always be performed, with a statistical input such that input parameters are varied in their ranges of expectation and probability distributions are determined for the output quantities.

### (d) Effect of the momentum correction factors *α*_{xx}, *α*_{xy}, *α*_{yy}

#### (i) Estimation of momentum correcting factors

It is our conjecture that for the ranges of shearing and sliding that are typical for avalanches the corrections implied by the incorporation of *α*-values different from unity are not significant. To estimate numerical values for *α*_{xx}, consider gravity-driven shear flow of a power law fluid down an inclined plane, infinitely long, steady and of constant height, allowing for a plane-parallel velocity that consists of a sliding velocity *u*_{s} and a shearing, that is given by , *n*≥1, where *A* and *n* are constant and *τ* is the shear stress. If *u*_{m} denotes the (maximum) surface velocity then(3.6)describes the velocity profile for *z*∈[0, *h*]. If we let *u*_{s}=*Χu*_{m}, *Χ*∈[0,1], then (3.6) becomes(3.7)and describes a range of velocity profiles from ‘only shearing and no sliding’ (*Χ*=0) to ‘only sliding and no shearing’ (*Χ*=1) and all intermediate stages (0<*Χ*<1) for any power law rheology, i.e. Newtonian (*n*=1), Bagnoldian (*n*=1/2) and even the non-realistic case (*n*=0).

It is easy to compute *α*_{xx}≔*α* from (3.6) or (3.7), the result being(3.8)for which values are given in table 1. Note that for *Χ*=1, *α*=1, irrespective of the value of *n*<∞, and *α*(*Χ*, ∞)=1 for all *Χ*, which is often associated with ideal plastic behaviour.

It is difficult to find estimates for *α*_{xy} and *α*_{yy}. For this reason, analyses are conducted with *α*_{xy}=*α*_{yy}=*α*_{xx}=*α*. We shall also make this choice, but say more on this below.

#### (ii) Computed avalanche flow for *α*_{ij}≠1

Computations were performed for *α*_{ij}=*α*, for all *i*, *j*, and values *α*=1, 1.05, 1.1, 1.2, on a bottom profile comprising an inclined plane merging through a cylindrical transition zone into a horizontal plane. The granular material is suddenly released so that the bulk material commences to slide down an inclined flat plane at 40° into a horizontal run-out plane connected by a smooth transition. The computational domain is the rectangle *x*∈[0,30] (along the downslope direction) and *y*∈[−7,7] (in the cross-slope direction) in dimensionless length units, where the inclined section lies in the interval *x*∈[0,17.5] and the horizontal region lies where *x*≥21.5 with a smooth change in the topography in the transition zone, *x*∈[17.5,21.5]. The corresponding thickness contours of the avalanche and the thickness profiles along the central line *y*=0 at three different dimensionless times *t*=6, 12, 24 are displayed in figures 4 and 5, respectively. With the deviation of the correction factor *α* from unity, the flow becomes slower, but its spreading increases. For *α*=1, a final deposition can be reached at *t*=18, but for larger values of *α*, it will take longer (e.g. *t*=24 for *α*=1.2). Comparison of depositions with different values of *α* shows that, although the shapes of the final depositions are fairly similar, the deposition with a large value of *α* extends over a larger area, the position of its maximum height is somewhat behind that for *α* close to unity and the maximum height is decreasing with the increase of *α*.

The above calculations were performed for momentum correcting factors having all the same values. We have repeated computations also for *α*_{xx}=*α*, but *α*_{xy}=*α*_{yy}=1, and found only very small differences between the results of figure 4 and those corresponding to figure 4 for the new *α*-values. Moreover, figure 4 indicates that for 1<*α*<1.1 consideration of shearing within the granular pile may safely be ignored. We have also indicated this range in table 1 by writing those *α*-values in boldface for which a momentum correction due to shearing might be necessary. Realistically, since granular shear flows behave much more like a Bagnold-fluid rather than a Newtonian fluid this is only the case when the no-slip boundary condition applies or . Since the no-slip boundary condition hardly ever occurs, these are very rare cases. Moreover, even in these cases ‘errors’ implied by the use of the SH equations with *α*=1, are very small.

We conclude, a depth-averaged model based on uniform velocity profiles provides a fair description of the dynamics of flow avalanches even on rather rough beds.

## 4. Outlook

This brief review with some additional new results and discussions of the depth-integrated SH equations allows the following inferences to be drawn:

For laboratory avalanches starting on steep slopes, the model seems to reproduce the motion of a finite mass of granular material down inclines from initiation to run-out pretty well, including shock capturing features. Additional studies are necessary for those cases in which obstructions are hit in the vicinity of which the shallowness assumption is formally violated.

Comparison of the model equations with field events—either snow or rock avalanches—are scarce and insufficiently conclusive. Post-calculations of the Madlein avalanche in Austria seems to indicate that the Coulomb basal friction law is insufficient and requires complementation by a viscous contribution. However, to make this statement conclusive additional comparisons are needed.

The model equations have not yet been tested for stresses exerted on walls of obstructing objects. Are computed wall pressures and shear tractions representative for the corresponding tractions in the experiment? Practically, these quantities are important ones, and suitability of the model equations for these would make the model much more valuable. This request holds equally true for all other avalanche models known to us.

There is a need for a wider application of the model to situations with natural topographies. Such topographies should also be reproduced in physical models at smaller scale to test the robustness of the equations as well as the reproducibility of avalanche motion performed with these. A first attempt is given by Iverson

*et al*. (2004).Avalanches starting from gentle slopes of inclination angles less than or equal to 30° and on bumpy beds behave differently, and they often come to a premature stop. They are generally less frequent, since in reality bumpy beds are not the rule, and they are small and consequently less dangerous.

Finally, it ought to be mentioned that this review and all the work reported herein does not touch the practically significant entrainment mechanism from the ground. If snow avalanches move over a layer of deposited snow, they often entrain snow from this layer and grow in mass. This often happens and significantly influences the dynamics of an avalanche. To our knowledge, there is only one paper attacking this problem (Sovilla *et al*. 2001). Of course, in real avalanche events the measurement of the entrained mass is difficult and in most situations not controllable. Needless to say, this adds further to the uncertainty in comparisons of real avalanche events with computational output of any model.

## Acknowledgments

This work was done while the first author was at the Isaac Newton Institute of the University of Cambridge as a member of the programme ‘Granular and particle laden flow’. Financial support is gratefully acknowledged.

## Footnotes

One contribution of 11 to a Theme ‘Geophysical granular and particle-laden flows’.

↵The more common denotation is ‘viscous’ contribution.

↵The centre of mass of a finite mass is moving steadily; the avalanche can still deform.

↵In these experiments

*δ*and*ϕ*were similar to those also relevant in dense flow avalanches.↵All these schemes reproduce the analytical one-dimensional parabolic cap similarity well (Tai 2000); however, in regions of transitions from supercritical to subcritical flow large oscillations may falsify the numerical solution. In this sense the Minmod TVD limiter has the best performance; Wang

*et al*. (2004).- © 2005 The Royal Society