Introduction
Since quantitative measurements of the mass balance of snow avalanches have been reported (Reference Issler, Gauer, Schaer and KellerIssler and others, 1996; Reference Sovilla, Sommavilla and TomaselliSovilla and others, 2001; Reference Vallet, Gruber and DufourVallet and others, 2001; Reference SovillaSovilla, 2004), modellers have renewed their efforts to describe the important but poorly understood process of entrainment. Different erosion mechanisms have been hypothesized (Reference Gauer and IsslerGauer and Issler, 2004), and for some of these (frontal ploughing, ripping of chunks corresponding to an entire snowpack layer and continuous abrasion along the flow bottom) there seems to be conclusive experimental evidence from the deposit structure or data from profiling radar (Reference Gubler and HillerGubler and Hiller, 1984; Reference Issler, Hutter and KirchnerIssler, 2003; Reference SovillaSovilla, 2004). To date, the potential of laboratory chutes for detailed studies of erosion and entrainment mechanisms with diverse bed and flow materials under controlled conditions has barely been tapped (Reference Barbolini, Biancardi, Cappabianca, Natale and PagliardiBarbolini and others, 2005).
Due to the lack of conclusive experimental knowledge, mathematical models of bed entrainment have mostly been heuristic (see Reference Eglit and DemidovEglit and Demidov, 2005, for a comparison of a large variety of entrainment models developed for different phenomena), typically of the form
where a and b are empirical exponents, σ b is the bed shear stress and τ c and u c are the threshold, or critical, bed shear strength and velocity, respectively. More physical approaches rooted in fluid dynamics were proposed by the Moscow school, either treating the flow/bed interface as a shock front (Reference Briukhanov and OuraBriukhanov and others, 1967; Reference Grigorian and OstroumovGrigorian and Ostroumov, 1977) or postulating an analogy with shear-induced mixing at the interface between different fluids (Reference Eglit and ShahinpoorEglit, 1983). In the case of (frictional, dry) granular flows on a bed of the same material, Reference Boutreux, Raphaël and de GennesBoutreux and others (1998) clarified the role played by the slope angle. Their work was extended and applied to snow avalanches by Reference Naaim, Naaim-Bouvet, Faug and BouchetNaaim and others (2004),who may have been the first to recognize the close connection between the shear rate at the bed/flow interface (governed by the rheology of the flowing material) and the entrainment rate.
D. Issler and T. Jóhannesson (unpublished information) further explored the connections between flow rheology, bed properties, erosion of bed particles and their entrainment into the flow, assuming constant flow depth and a perfectly brittle bed material that is characterized by its shear strength, τ c, and is entrained along the flow bottom. With these simplifications, they found analytical solutions for quasistationary flow of a Bingham fluid and realistic erosion rates for parameter values typical of snow avalanches. However, extension to non-stationary situations requires some approximations that are difficult to control a priori, and analytic solutions for the velocity and shear-stress profiles under entrainment may not exist for more complicated rheologies.
Our aim here is therefore to create a tool with which to find suitable approximations to the velocity profiles of entraining flows with nonlinear rheology, or to directly obtain parameterizations of the erosion rate as a function of the depth-averaged flow velocity, 
and the flow depth, h. We formulate and implement a one-dimensional numerical model for the time evolution of the velocity and shear-stress profiles, the erosion rate and the flow depth that can, in principle, be applied to any rheology of the form
where 
and 
are the specific shear and normal stress, respectively, and ρ is the density. The measure of shear, 
reduces to the shear-rate gradient, ∂zux , in our idealized situation.
Evolution Equation For the Velocity Profile
For the sake of simplicity, we make several assumptions, which can easily be relaxed or modified later: (1)We treat the avalanche as infinitely long in the x-direction along a slope inclined at a constant angle, θ. All variables thus depend only on the slope-perpendicular coordinate, z, and the time, t (Fig. 1). (2) the snow cover and the avalanche are assumed to have the same density, ρ. (3) Entrainment stops as soon as the entire snow cover is eroded away. (4) the combined depth of avalanche and snow cover remains constant, i.e. there is no mass redistribution along the body of the avalanche. As a consequence, a particle in the avalanche always moves parallel to the ground and the interface between flow and bed is also parallel to the ground.
Fig. 1. Schematic t–z diagram, showing the erosion front advancing in time and the choice of coordinate system. Note that the x-direction is not shown. The slope of the curve, b(t), is the (non-material) erosion speed, w e(t). There is no mass flux in the z-direction even though mass flows across the moving interface.
Mass balance is trivially satisfied in this set-up, and the momentum balance in the z-direction yields 
if we assume the surface at z = 0 to be stress-free and (g
∥, g⊥
) ≡ (g sin θ, g cos θ) is the gravitational acceleration vector. With u(z, t ) the velocity in the x-direction, we write the balance of x-momentum as
The shear stress is to be expressed as a function of the shear rate, ∂zu (cf. Equation (2)), so Equation (3) in general is a nonlinear partial differential equation of diffusion type for u. Note that 
< 0 for u > 0 in our coordinate system.
We assume no slip at the erosion front, i.e. u(b, t) = 0. To completely specify this problem with a time-dependent computational domain, the erosion speed, w e = ḃb(t ), must be determined. For bed materials that are perfectly brittle (implying that they fail instantaneously and without significant strain once the stress reaches a critical value, −τ c) Issler and Jóhannesson (unpublished information) argued that there is a fundamental feedback mechanism that keeps the bed shear stress, τ b, at the value −τ c as long as entrainment is possible. If |τ b | were smaller than τ c, erosion would cease, so the shear rate and |τ b | would increase until erosion resumes. If |τ b | > τ c, however, the erosion rate would grow and the acceleration of more eroded bed material would reduce |τ b |. We assume that snow can indeed be modelled as a perfectly brittle material to a sufficient degree of accuracy, given the rapid loading induced by the arrival of the avalanche and the substantial shear stresses exerted by the flow, which are of the order of hρg ∥ ≈ 1–3 kPa.
Equation (2) and the boundary condition 
imply that the shear rate at the erosion front is constant in time:
The progression of the erosion front is driven by the need to maintain the stress boundary condition while the recently eroded material is being accelerated. At time t + dt, the velocity of a particle that was eroded at time t is, according to Equation (3) and to first order in dt ,
while the erosion front is now located at
The shear rate at the front thus evaluates to
where 
As this must be equal to 
we find the propagation speed of the erosion front in the z-direction:
Note that 
is time-dependent if the shear stress depends on the normal stress (e.g. due to Coulomb friction). As an example, apply the NIS (Norem–Irgens–Schieldrop) model (Reference Norem, Irgens and SchieldropNorem and others, 1989),
to Equation (8); here μ is the dry-friction coefficient, 
the specific yield strength, 
the specific effective normal stress transmitted by frictional grain contacts, and ν, ν'
are the coefficients of the dispersive stresses (m2). The erosion speed becomes
Numerical Implementation and Simulation Results
Much of the simplicity of Equation (3) is due to the absence of advection in the z-direction – a feature that the numerical scheme should reflect. At the same time, there is a moving (non-material) boundary at z = b(t ) where increased precision is required because the gradient of the shear rate must be computed at the erosion front (cf. Equation (8)). In order to accommodate both requirements, we combine a coarse regular fixed grid of mesh length Δz c with a fine grid of N + 1 equidistant nodes a distance Δz f = Δz c /N apart in a finite-difference scheme. At the beginning of each time-step, the fine mesh is translated so that its end coincides
with the present location of the erosion front, and u, 
and 
are interpolated to the new node positions (see Fig. 2). The field values at the coarse-grid node within the fine grid are also interpolated; as the erosion front progresses, successive coarse-grid nodes are activated.
Fig. 2. Schematic representation of the coarse grid and the superposed fine grid, which is shifted to the new position of the erosion front at the beginning of each time-step.
We use staggered grids, with the velocity, u, located at the nodes and the stresses, 
at the cell centres. The z-derivatives are approximated by central differences, and quadratic extrapolation is used for 
so that the scheme is second-order accurate in space. Time-stepping is by simple forward differencing and thus first-order accurate. The time-step is determined from two criteria: (1) the eroded depth per time-step must be less than the cell width of the fine grid; (2) disturbances of the velocity profile, which are damped in the exact solution due to the diffusive nature of the differential equation, must not be amplified in the numerical approximation. The second criterion often imposes very short time-steps, but computation times are nevertheless short (a few seconds) on a current personal computer.
If the erosion rate is a little too large due to small numerical errors, the specific shear stress will drop slightly below 
and erosion ceases at once. In order to prevent spurious fluctuations of this kind, the erosion rate is set to 98–99% of its calculated value.
Another numerical problem, visible in Figure 3b and e, arises from the displacement of the fine grid relative to the coarse one: as the tail node of the fine grid is shifted past a node of the coarse grid (which then gets activated), the velocity and shear stress gradients change a little, but they do so abruptly. The disturbance propagates to the erosion front and causes the erosion rate to be overestimated briefly. A more precise interpolation scheme is called for, but the error is negligible for all practical matters.
Fig. 3. Results of numerical simulations for the case of Newtonian (left column) and Bagnoldian (right column) fluids in laminar flow. The initial conditions and material properties are summarized in Table 1. The rows show the evolution of, respectively, velocity and acceleration, erosion speed and flow height, and velocity and shear-stress profiles. Note that the Bagnoldian fluid meets the erosion criterion only after ∼6 s when it has almost attained its terminal velocity without erosion of ∼17ms − 1. The profiles differ strongly from the corresponding equilibrium profiles of stationary non-entraining flows with the same rheology.
The code was validated by verifying that non-entraining flows with various rheologies converge to the known analytical solutions for the terminal velocity and stationary velocity profile. In addition, the simulated time evolution of the velocity profile of a non-entraining laminar Newtonian flow matches the theoretical solution obtained by Fourier expansion.
Two exploratory simulations with the new model (Table 1; Fig. 3) assume Newtonian and Bagnoldian rheologies, respectively. The latter is obtained from Equations (9) and (10) as a special case by letting 
μ → 0. Both simulations deliberately allow unrealistic flow height growth in order to highlight the asymptotic behaviour of the model. If the finite avalanche length and mass redistribution along the flow direction were accounted for, realistic erosion depths of 0.5 –2 m and lower velocities and flow heights would result. While the Newtonian flow was started with an initial velocity and eroded a relatively weak snow cover from the beginning, the snow cover strength for the Bagnoldian case was only slightly below the equilibrium bed shear stress of the non-entraining flow. Erosion therefore began only when the flow had nearly reached terminal velocity, the acceleration had become small and the specific bed shear stress approached the equilibrium value, h
0
g
∥. However, after the threshold was passed, the growing flow height made more gravitational force available for entraining snow and accelerating the flow.
Table 1. Initial conditions and material properties of the simulations shown in Figure 3
The long-term behaviours of both flows exhibit striking similarities, in that both the acceleration and the erosion speed tend towards constant values and the velocity profiles are almost linear except near the surface, where they resemble somewhat the equilibrium profiles of non-entraining flows of the same rheology. Accordingly, the shear-stress profiles are only weakly curved near the surface, but ∂z τ differs strongly from the stationary value, ρg ∥ = 1500 Pam − 1. The observed behaviour is a consequence of the crucial postulate that the bed shear stress equal the bed shear strength. This is made clearer in the following section.
Asymptotic Solutions For Slab Models With Entrainment
The results obtained through the numerical solutions suggest there exists an asymptotic solution to the equation of motion of an entraining mass point on an inclined plane with constant entrainment rate and linearly growing flow depth
and velocity under the shear-stress boundary condition, 
We therefore set
in the equation of motion (overlined quantities are depth-averaged)
and obtain
One immediately sees that 
and the bed shear stress, 
becomes insignificant compared to the gravitational traction, 
This conclusion is borne out by the numerical simulations with Newtonian and Bagnoldian rheologies (Fig. 3).
The most notable features of this solution are (1) its independence of the rheology, (2) its apparent stability and (3) the fact that asymptotically one half of the momentum input from gravity goes into accelerating the bulk of the flow and the other half into bringing the eroded mass up to the average flow velocity.
As mentioned above, some (temporary) model assumptions are unrealistic for snow avalanches, and hence so is this asymptotic solution with ever-increasing velocity and flow height. Interestingly, however, Reference Gauer, Kronholm, Lied, Kristensen and BakkehøiGauer and others (2010) pointed out that available experimental data on maximum velocity of avalanches, u
max, plotted against the drop height, H, indicate the relationship u
max
∝ (gH)1/2, which is obtained if the effective friction coefficient is approximately constant, but not if it grows significantly with velocity. The boundary condition on the shear stress we have adopted in this work, 
gives an effective friction coefficient that is independent of the velocity, but diminishes as the inverse of the growing flow depth. However, entrainment of eroded snow at an ever-increasing speed gives rise to a resistive acceleration, 
that mimics dry friction as long as the slope angle is constant and there is enough erodible snow. The dynamics of entrainment thus offers a possible explanation for the observed growth of u
max with H that can be tested further against experimental data.
Conclusions and Outlook
The numerical model described here has proved to be a very useful tool to continue the work of Issler and Jóhannesson (unpublished information) beyond the restrictive assumptions needed to obtain analytical solutions. In particular, the asymptotic solution hinted at by the numerical simulations highlights the decisive dynamical role played by the postulated feedback mechanism that locks the bed shear stress at the value of the shear strength of the snow cover, and by the assumptions concerning the redistribution of entrained snow. We anticipate that future detailed studies of the shape evolution of the velocity and shear-stress profiles will give useful guidance in the development of approximate entrainment models based on the self-consistent approach of Issler and Jóhannesson (unpublished information). It is also conceivable that further asymptotic solutions can be found with different assumptions regarding redistribution of entrained snow.
The numerical model also has potential for future development in its own right. Some of the next steps we plan to implement are: (1) improved interpolation methods to reduce disturbances as the fine grid leaves a node of the coarse grid behind; (2) extension to a slab model on an arbitrary path profile, with the slab length growing as snow is entrained; (3) snow strength varying along the path and with depth in the snow cover; (4) the ideas leading to Equation (8) can and should also be incorporated in more advanced numerical models in two dimensions (x and z). However, a more refined formulation of the momentum balance will be needed because the assumptions of vanishing uz and independence of x are no longer valid.
Last but not least, there is urgent need for innovative experimental studies of the fundamental entrainment mechanisms. Particularly, our central assumption of perfectly brittle failure needs to be confronted with controlled experiments under realistic conditions. If it holds to some degree in snow, the conditions for its validity in other eroding gravity mass flows are a question of great practical interest. While the experimental challenge is substantial, the pioneering study of Reference Barbolini, Biancardi, Cappabianca, Natale and PagliardiBarbolini and others (2005) indicates a promising approach.
Acknowledgements
D.I. gratefully acknowledges a fellowship from the NILS Mobility Project (ABEL Extraordinary Chair action) and the hospitality extended by the Laboratorio de Geotecnia of CEDEX in Madrid, Spain. Completion of this work was financed by the EU Framework Programme 7 project SafeLand and by the Norwegian Water Resources and Energy Directorate through the research project Snow Avalanches. We thank two anonymous referees for careful reading and constructive criticism. This is contribution No. 349 from the International Centre for Geohazards.
