## Abstract

The idea that snow avalanches might behave as granular flows, and thus be described as Coulomb fluid flows, came up very early in the scientific study of avalanches, but it is not until recently that field evidence has been provided that demonstrates the reliability of this idea. This paper aims to specify the bulk frictional behaviour of snow avalanches by seeking a universal friction law. Since the bulk friction coefficient cannot be measured directly in the field, the friction coefficient must be calibrated by adjusting the model outputs to closely match the recorded data. Field data are readily available but are of poor quality and accuracy. We used Bayesian inference techniques to specify the model uncertainty relative to data uncertainty and to robustly and efficiently solve the inverse problem. A sample of 173 events taken from seven paths in the French Alps was used. The first analysis showed that the friction coefficient behaved as a random variable with a smooth and bell-shaped empirical distribution function. Evidence was provided that the friction coefficient varied with the avalanche volume, but any attempt to adjust a one-to-one relationship relating friction to volume produced residual errors that could be as large as three times the maximum uncertainty of field data. A tentative universal friction law is proposed: the friction coefficient is a random variable, the distribution of which can be approximated by a normal distribution with a volume-dependent mean.

## 1. Introduction

This paper is concerned with the probabilistic calibration of an avalanche dynamics model. The objective of probabilistic calibration is to predetermine the dynamic features (e.g. the velocity, impact pressure and geometry) of extreme avalanches for various land-use planning and engineering applications, including avalanche zoning and mitigation. Although snow avalanches are complex phenomena involving varied processes from snowfall to snow friction, there is growing evidence that their chief dynamic features can be described by using simple deterministic models (Hutter 1996).

In the earliest developments of the scientific study of avalanches in the 1920s, the French forest engineer Mougin (1922) suggested using a sliding block to model the avalanche motion. Assuming that the block experiences a Coulomb frictional force, he was able to compute the maximum velocity reached by the avalanche and the point of furthest reach. The idea was further used in a few engineering applications, but the question of how to adjust the friction coefficient remained open at that time. Since then, different models of the frictional force have been proposed but most of them include a velocity-dependent contribution plus a constant term that can be interpreted as a Coulomb frictional contribution or a yield stress. Despite various attempts to find physical justifications for their expressions, these models remain speculative and suffer from lack of testing with field data.

The interest in the Coulomb model for describing avalanches and debris flows has been renewed recently (Iverson & Denlinger 2001; Tai *et al.* 2001). Dent (1993) argued that in snow avalanches, shear should occur in a thin layer at the avalanche base and that this behaviour is close to that observed in a granular flow down a smooth surface. Savage & Hutter (1989, 1991), Hutter *et al*. (1995) and Pudasaini & Hutter (2003) thoroughly studied the motion of a finite mass of grains along curved and smooth chutes in the laboratory and showed that including a basal Coulomb friction as bottom shear stress in depth-averaged equations makes it possible to describe the behaviour of granular avalanches for a wide range of flow conditions (see Hutter *et al*. 2003; Pudasaini & Hutter 2003 for a recent review on this theory).

Ancey & Meunier (2004) performed a back analysis on 15 well-documented avalanches by inferring the bulk frictional force from avalanche velocity. They found that for 11 events, the Coulomb model provided a fairly good description of the velocity variations along the path profile. Gubler (1993) took measurements on real avalanches by using a Doppler radar. He found that the velocity profile inside the observed avalanches exhibited a plug flow (constant-velocity zone) and a sheared zone at the bottom, clearly revealing that there was shear localization at the bottom. These experimental results provided evidence that the bulk rheological behaviour of flowing avalanches can be described by using the Coulomb model, which makes it possible to use granular-flow theories such as the Savage & Hutter (1989) theory in engineering applications. This would be a real breakthrough given the speculative character of current avalanche dynamics models. Indeed, the ad hoc character of the rheological equations used in these models is considered the main deficiency of the mechanical approach.

For the Savage–Hutter theory to be of practical interest in engineering, the Coulomb friction coefficient must be specified. There is no means of measuring the friction coefficient for real avalanches, so field data and calibration techniques must be used to determine it. As for any model calibration, questions emerge regarding (i) the nature of the parameter to be fitted, (ii) the proper selection of the calibration technique and (iii) the universality of the adjustment, as follows.

The fact that a variable cannot be measured directly and physically leads to questioning whether the variable reflects a physical process or merely an abstract conceptual picture of reality. This is crucial since the variation range of the parameter must be specified or checked. For instance, a friction coefficient

*μ*=2 is questionable if the coefficient has a physical meaning because it outweighs the upper bound of possible values, but it may be acceptable from a conceptual viewpoint.Determining the friction coefficient from field data basically involves solving an

*inverse problem*. An inverse problem groups different mathematical issues to invert a relationship in the form where*Y*is a model output or parameter that can be measured,*X*is a model input or an internal variable that is not known and is a functional. Different approaches to solving this problem can be applied, as follows.Deterministic methods try to find a solution in the form , but difficult theoretical and numerical issues are usually encountered: existence, uniqueness and stability of the solution are rarely ensured in physical problems. These problems are usually exacerbated when numerical schemes induce discretization errors or when observed data

*Y*are noisy.Stochastic methods are a second family of techniques which attempt to find

*X*by a trial-and-error procedure: guessed values (*X*_{n}) of*X*are randomly generated and then tested and selected so that*F*(*X*_{n}) converges to*Y*. Stability is enforced since*X*is determined without making use of , but other difficulties such as slow convergence can arise. Stochastic methods are particularly useful when data are noisy and uncertain, which is the case in the available avalanche data.

We stress that a model is statistically admissible if its residual error, namely the difference between the observed and computed values, lies below the datum uncertainty. If a model parameter has a physical meaning, similar experiments or events must lead to the same value (to the model precision) when measuring and calibrating it. Alternatively, there is no clear reason to expect universality if a variable is only conceptual. For instance, if a model is calibrated for a given site, it is unlikely that it holds true for another site without further calibration.

This paper begins by presenting the equations of motion used to model avalanche motion and the probabilistic methods needed for determining the friction coefficient (§2); this paper extensively uses Bayes' theorem and inference techniques (Brooks 1998; D'Agostini 2003). The raw data and the preliminary treatment (path profile regularization) is presented in §3 and we test different assumptions regarding the friction coefficient in §4. We show that this coefficient behaves as a continuous random variable dependent on avalanche volume. Further analysis reveals that the friction coefficient may depend on other parameters, but, in the absence of data corroboration, a probabilistic description of friction depending on volume can be proposed.

## 2. Methods

### (a) Equations of motion

In most current numerical simulations, snow avalanches are modelled by using the shallow-flow approximation (Brugnot & Pochat 1981; Eglit 1983; Bartelt *et al.* 1999; Barbolini *et al.* 2000). Although these recent developments are undoubtedly a promising approach in modelling snow avalanches, their level of sophistication contrasts with the crude current knowledge of snow rheology and avalanche physics as well as the poor quality of available field data. Since this investigation is a preliminary step towards the determination of the frictional coefficient, we restrict ourselves to a simplified model that requires less computational effort and time: the sliding-block model. Despite its simplistic character, this model successfully describes the main characteristics of flowing granular mass, as shown by laboratory experiments of Iverson *et al*. (in press) and by the back-analysis from field data produced by Ancey & Meunier (2004).

Here, an avalanche will be considered a slender sliding rigid body of volume *V* and mass *m*. The following assumptions are made.

The variations in body shape are ignored as a first approximation (A1).

The body is assumed to move along a curvilinear two-dimensional profile (A2) of which an equation in a Cartesian frame takes the form

*y*=*f*(*x*), where*y*is the elevation and*x*is an arbitrary distance measured along a horizontal axis.There is no significant lateral spreading of the mass (A3) and there is no sudden kink in the avalanche trajectory (A4).

The longitudinal profile is a smooth and gently varying curve (A5). Basically, this means that the curvature radius

*R*is at least as large as the typical length-scale*L*of the flow or larger:*R*∼*L*or*R*>*L*.The sliding body experiences a Coulomb frictional force, that is, the normal and tangential components of the frictional force

*F*_{n}and*F*_{t}are linearly dependent:*F*_{t}=*μF*_{n}, where*μ*is the Coulomb friction coefficient.The avalanching material is of finite and constant mass.

The position of the centre of mass is given by its curvilinear abscissa (where *f*_{x} is the *x*-derivative of *f*). Therefore, we have *x*=*ξ* cos *θ*, with *θ* being the mean path inclination computed over the interval [0, *x*]. The ordinate of the centre of mass (relative to the curve *f*) is denoted by *η*. In the natural basis (**e**_{1}, **e**_{2}) associated with the curvilinear coordinates (*ξ*, *η*), the contravariant components of the velocity vector are denoted by (*u*^{(1)}, *u*^{(2)})=(d*ξ*/d*t*, d*η*/d*t*) and its physical components are given by (*u*^{〈1〉}, *u*^{〈2〉})=((1−*η*/*R*)*u*^{(1)}, *u*^{(2)}) (Sedov 1975). The contravariant components of acceleration in the natural basis arewhere are the Christoffel symbols. Because the natural basis is orthogonal, the Christoffel coefficients are zero, except for and , where *C*=1/*R* is the curvature. The velocity in the *ξ*-direction is *u*=*u*^{〈1〉}=(1−*ηC*)d*ξ*/d*t*; assumptions A1 and A2 imply that *η* is fairly constant and the velocity *u*^{〈2〉} in the *η*-direction is close to zero. The downward and normal components of the momentum equation can be expressed in the physical curvilinear basis as(2.1)

(2.2)The first term on the left-hand side of equation (2.1) represents the downward component of the acceleration, while the second term reflects radial effect due to the curvature of the path profile. The first contribution on the right-hand side of equation (2.1) is the driving action of gravity, while the second term stands for the frictional force exerted by the bottom (ground or snowcover) of the avalanche. Using the Coulomb model *F*_{t}=*μF*_{n}, we expand the equations of motion with respect to the small parameter *ϵ*=*η*/*R* and neglect terms of the order of *ϵ* and higher degrees. Therefore, the equations of motion can be simplified into(2.3)This equation holds for the avalanche's centre of mass, and for any other point because of assumption A1. Notably, if *ℓ* denotes the curvilinear distance between the leading edge *ξ*_{f} and the mass centre *ξ*, it is possible to replace *ξ* by *ξ*_{f}−ℓ in equation (2.3), which provides the equation governing the leading-edge variations with time.

### (b) Estimation of the probability density function of *μ*

The equation of motion (2.3) being known, it is possible to compute the run-out distance *x*_{stop} (point of farthest reach) and its elevation *y*_{stop}=*f*(*x*_{stop}) for a given path and friction coefficient. Repeating the procedure for different values of *μ* makes it possible to obtain a numerical estimate of the function *y*_{stop}=*F*(*μ*) or, if this function is monotonically decreasing, a numerical estimate of *μ*(*y*_{stop}). The function *y*_{stop}(*μ*) is monotonically decreasing when the path profile is monotonically decreasing (*f*_{x}<0 along the path).

We solved equation (2.3) for different values of *μ*, usually in the range 0.2–0.7 and with an increment of 0.0025. For each increment of *μ*, we computed the run-out distance. The resulting array of values was then interpolated to yield either an estimate of *y*_{stop}(*μ*) or *μ*(*y*_{stop}).

When a time-series of run-out distance has been recorded for a given path, we can use the function *μ*(*y*_{stop}) to obtain the *μ* values for all the recorded events. If we are interested in determining the probability distribution of the *μ* sample, a common practice is binning, that is, grouping data into categories. The resulting histogram then provides an estimate of the probability density function of *μ*. Here, because of the limited size of the samples, different choices of bins lead to drastic changes in the shape of the histogram. It is thus better to use a more refined approach. We used the approach developed by Holy (1997).

A sample of *μ* values is available and the probability distribution (*Q*) that best approximates the empirical distribution of *μ* is sought. The candidate functions of *Q* belong to the space of continuous, normalized and positive functions, that is, *Q*=*ψ*^{2} and . In addition, *Q* is assumed to be smooth or, from a probabilistic viewpoint, the probability of observing large gradients in *ψ* is very small. Holy (1997) synthesized this information by stating that the *a priori* distribution of *Q*(*P*[*Q*]) iswhere *Z* is a normalization factor, *ℓ* is a free parameter controlling the smoothness of the distribution and *δ* is the Dirac function. Applying Bayes' theorem enables us to consider the data in order to improve the knowledge of the distribution *Q*,where is the *a posteriori* probability of *Q* given the data, is the *likelihood* of the data —knowing the distribution *Q*—and the denominator is a normalization factor, . We then obtainwhere *S* is the functional

The most likely distribution, given the data, is the function that minimizes the functional *S*. It can be shown that the solution has the form (Holy 1997)where *κ*^{2}=2*λ*/*ℓ*^{2}, with *λ* being a Lagrange multiplier that must be adjusted to enforce the normality of *Q* and *a*_{i}(1≤*i*≤*N*) are coefficients. The Lagrange multiplier *λ* and the coefficients *a*_{i} are solutions of the system of *N*+1 equations(2.4)The free parameter *ℓ* controls the penalty supported by large gradients in *Q*. Better agreement is obtained when a ‘small’ value of *ℓ* is chosen, but in this case, the solution may be widely oscillating. By contrast, the best smoothness is produced when a relatively ‘large’ value of *ℓ* is used, but the resulting function may poorly represent the data.

### (c) Seeking the dependence of *μ* on *V*

In their rheological back-analysis of 15 well-documented snow avalanches, Ancey & Meunier (2004) showed that in many cases, the snow avalanche motion can be modelled with the Coulomb sliding-block model. By adjusting the *μ* value for the computed run-out distance to match the recorded value, the velocity variations along the path profile can be fairly well described over a large part of the path. Figure 1 reports the adjusted value of the friction coefficient *μ* as a function of the snow volume involved in the avalanche for the 11 events that were considered by Ancey & Meunier (2004) to reach a Coulomb regime. The plot seems to imply a reduction of the friction coefficient *μ* with growing volume *V*, which is similar to that observed for other geophysical mass flows (Savage 1989; Legros 2002). It is attractive to quantify this decrease by adjusting a curve through the points, but as shown in figure 1 with three curves, the data scattering is too significant and the volume range is too limited (two orders of magnitude) to deduce the most suitable trend.

Additional field data should provide better insight into the possible dependence of *μ* on volume. A major difficulty in considering additional data comes from their large uncertainty.

Snow is a compressible material; typically, the ratio between the released and deposited snow densities can be as high as three.

In most cases, the deposited volume is estimated by the naked eye by evaluating a typical length, depth and width in the field. Depending on the topography, relative errors exceeding 100% may be induced when using this crude estimation technique.

Significant uncertainty is also met when estimating the run-out distance (see §4

*a*).

With such prerequisites, attention must be paid to data uncertainty. A deterministic method such as the least-square method is not sufficient to fully address the problem; hence the need for using more elaborate techniques. The Bayesian approach offers an interesting way of investigating both possible dependences of *μ* on *V* and the induced error in adjusting a curve to field data (D'Agostini 2003). Let us now formulate our problem in the following way.

There are errors introduced by model approximations, as well as errors in estimating the run-out distance and avalanche volume. We assume that both error sources can be addressed using a single error parameter. We then introduce the deviation

*ϵ*between the computed (*y*_{stop}) and recorded (*y*^{obs}) run-out elevations:*ϵ*=*y*^{obs}−*y*_{stop}.We further assume that

*ϵ*is a random realization from a normal distribution of mean 0 and unknown variance . Given the volume, we assume that there is a one-to-one correspondence between the friction coefficient and the run-out elevation:*y*_{stop}=*F*(μ(*V*)).We assume that we have an

*a priori*idea of the*μ*dependence on*V*:*μ*=*G*(*V*;*Θ*), where*Θ*denotes the free-parameter set of the functional*G*(e.g. for a power-law dependence), and we can express*G*as*G*=*αV*^{−β}, with*Θ*={*α*,*β*}.We have a series of

*N*events for which run-out elevations and volumes were observed: . The Bayes rule allows us to update the parameters*Θ*with the data and to quantify the uncertainty on*Θ*,(2.5)where, in the numerator,*P*(*Θ*) and*P*(*σ*) refer to the probabilities or*priors*of*Θ*and*σ*. The quantityis called the ‘likelihood’ and refers to the probability of observing the samplewhen the functional**d***G*, its parameters*Θ*and the s.d.*σ*are known. The denominator is a normalizing constant. Basically, the Bayes rule is an updating process, where our knowledge of*Θ*and*σ*is updated using the available informationto provide the**d***posterior*distribution*P*(*Θ*,*σ*|,**d***G*).In order to obtain an estimate of the best choice for the values of

*Θ*and*σ*, random values are drawn from the posterior distribution*P*(*Θ*,*σ*|,**d***G*). Finally, the best choice of (*Θ*,*σ*) can be made by determining the modes or the means of the posterior distribution. The marginal probability density function of*Θ*can be estimated by integrating the posterior probability with respect to*σ*, leading to an estimation of a confidence interval (credibility in the Bayesian approach) of*Θ*. Note that the strength of the Bayesian approach lies not only in a proper way of choosing the parameters*Θ*, but also in a realistic assessment of the uncertainty on*Θ*and the overall error (combining model and observation errors)*ϵ*.

What we need now is to specify the priors *P*(*Θ*) and *P*(*σ*) and the iterative procedure to draw random samples from the posterior distribution *P*_{p}=*P*(*Θ*, *σ*|*G*, * d*):

Analysing the data in figure 1 suggests possible priors for

*Θ*when a function*G*has been adjusted. A common assumption is to consider that the prior of*Θ*is a multivariate normal distribution of mean*Θ*_{0}and covariance matrix where the vector**Θ**_{0}may stand for the values obtained by adjusting the parameters**Θ**_{0}using the data in figure 1 and the least-square method. Similarly,*Σ*may be interpreted as the asymptotic covariance matrix determined in the regression.Sampling from the posterior distribution

*P*_{p}can be carried out by using the efficient and robust Metropolis–Hastings algorithm, which is based essentially on Markov chain sampling and Monte Carlo simulations (MCMC simulations; Brooks 1998; Robert 2001). We simply expose the principles of this method, which is widely applicable and easy to implement.

The basic idea of MCMC algorithms is to introduce a probability distribution from which sampling is straightforward, instead of directly sampling from the posterior distribution. We refer to this distribution as the *instrumental* distribution *q*. We will generate random samples from *q*, and explore the probability space of the posterior distribution; *q* is then a transition probability, which is used to move from a probability state of *P*_{p} to another one. Iterating the procedure leads to a sample of values, the empirical probability distribution of which is close to *P*_{p}. In practice, the following steps are performed.

Given a current state

*X*_{n}=*x*, draw a candidate value*y*^{*}from the instrumental distribution*q*(*y*|*x*).Define the acceptance rate

*r*as(2.6)Accept the value

*y*_{*}with probability*r*, i.e. draw a random value*u*from the uniform distribution . If*r*<*u*, accept*y*_{*}and set*X*_{n+1}=*y*^{*}; otherwise reject it and set*X*_{n+1}=*X*_{n}.Repeat the procedure.

Central to ensuring the efficiency of MCMC simulations lies in the proper selection of an instrumental distribution. Here, we have adopted the random-walk version of Metropolis–Hastings algorithm (Robert 2001), which involves selecting a symmetric probability distribution *q*=*q*(|*x*−*y*|). This choice leads to simplifying the expression of the acceptance rate *r* in (2.6): *r*=min[1, *P*_{p}(*y*^{*})/*P*_{p}(*x*)]. Convergence of the empirical distribution of (*X*_{n}) towards *P*_{p} is ensured here because of the exponential decrease of the tail of *P*_{p}. A common choice is to take an uncorrelated, multivariate, normal distribution with a tuneable covariance matrix The scale matrix *ρ* must be tuned such that there is a trade-off between the acceptance rate and the ability of the algorithm to fully explore the probability space. An accepted procedure for this algorithm version is to adjust *ρ* such that the acceptance rate *r* falls in the range of 0.25–0.5. Choosing an uncorrelated distribution *q* (i.e. *ρ* is a diagonal matrix) makes it possible to adjust the acceptance rate for each component of the vector (*X*_{n}).

In summary, we are looking for a nonlinear dependence of the friction coefficient *μ* on the avalanche volume *V* by adjusting parametric functions *G*(*V*|*Θ*) to field data. Assuming that we have an *a priori* knowledge (subjective or based on preliminary attempts at fitting *G*) of the distribution probability of *Θ*, Bayes rule (2.5) makes it possible to update this knowledge and to obtain an estimate of the precision *ϵ* of the model.

### (d) Least-square method revisited

Let us assume that we have a set of data * M*=(

*x*

_{i},

*y*

_{i})(1≤

*i*≤

*N*), and we are interested in adjusting a straight line across the corresponding points. This is a classical regression problem that can be solved (among other ways) by using the least-square method. When data are scattered around the straight line, the predictive capacity of the linear model is weak, and determining the probability of finding a point at a given distance from the line is needed. Here, we suggest dealing with this issue by using a Bayesian approach and conjugate distributions of the exponential families (normal random-effects model). We first assume that the deviation from the linear response is normally distributed:

*y*

_{i}=

*a*+

*nx*

_{i}+

*ϵ*

_{i}, where are identically distributed random values with zero mean and variance

*σ*

^{2}. We are trying to adjust (

*a*,

*n*,

*σ*) from data. Using Bayes' rule, one obtains(2.7)where we have assumed that the variables (

*a*,

*n*,

*σ*) are not correlated and

*P*[

*a*],

*P*[

*n*],

*P*[

*σ*] denote the prior probability distribution functions of

*a*,

*n*and

*σ*, respectively. Here, we assume, without loss of generality,We introduced

*τ*=

*σ*

^{−2}to simplify the computations. In the prior distributions, we also introduced the free parameters (

*a*

_{0},

*σ*

_{a},

*n*

_{0},

*σ*

_{n},

*γ*

_{σ},

*δ*

_{σ}). The Gibbs sampler method can be used to find the posterior probability

*P*

_{p}=

*P*[

*a*,

*n*,

*σ*|

*]. This technique first involves determining the full conditional probabilities of the random variables*

**M***a*,

*n*and

*σ*(Brooks 1998; Robert 2001). Here, we deduce from (2.7)

The Gibbs sampler algorithm operates as follows. We want to build a chain (*a*_{k}, *n*_{k}, *σ*_{k}), whose empirical joint distribution function converges towards *P*_{p}. At the first iteration, we initialize the values (*a*_{0}, *n*_{0}, *σ*_{0}) and then proceed by iteration from *k* to *k*+1 (*k*≥0) as follows.

Draw

*a*_{k+1}from*P*[*a*|*n*_{k},*τ*_{k}] with*τ*_{k}=*σ*_{k}^{−2}.Draw

*n*_{k+1}from*P*[*n*|*a*_{k+1},*τ*_{k}].Draw

*τ*_{k+1}from*P*[*τ*|*a*_{k+1},*n*_{k+1}] and setRepeat the procedure.

Since the full conditional probabilities are common distributions, it is computationally straightforward to generate random samples from them.

The best-fit parameters can be deduced by taking the mean or the mode of the generated samples (after removing the first elements, which are influenced by the initial conditions). Finally, we obtain the deterministic linear relationshipor the full probabilistic descriptionNote that the procedure that we have presented is quite general and the assumptions (linearity, prior selection, etc.) can be modified, usually leading to a substantial increase in the computational efforts. Notably, leaving aside linearity in the response and conjugate priors leads to coping with full conditional probabilities from which it may be intricate to generate random samples.

## 3. Data

### (a) Raw data

To supplement the small dataset (11 data) provided by Ancey & Meunier (2004), additional data were sought. Table 1 summarizes the origins of the 173 events and the main features of the seven paths where events were observed. Figure 2 shows the path profile together with the recorded run-out distances. All of the paths are located in the northern French Alps: Tours-en-Savoie and la Gurraz are two villages in Savoie; the other sites are situated in Haute-Savoie. Field data were extracted from a national database (*Enquête Permanente des Avalanches*). The time-series is not complete. For example, during the two world wars, no observation was reported. Avalanche volume was no longer measured after 1975 (even earlier for some paths). Moreover, the datum sample forms censored observations; for an avalanche to be observed or its deposit to be visible from the valley bottom, it must have travelled a sufficiently long distance from the starting zone. This implies that avalanches with small run-out distances (i.e. a high friction coefficient) could not be observed.

The database provides the run-out elevations and the typical dimensions of the snow deposit by reporting the length *L*_{v}, the width *W*_{v} and the deposit thickness *H*_{v}. Here, deposit volume was estimated as follows: *V*=(1/2)*L*_{v}*H*_{v}*W*_{v}; the multiplicative factor of one-half is arbitrary and was introduced to reflect the irregular cross-section of the avalanche deposits, which differs usually from a rectangular section. Uncertainty on the run-out elevation varies with time. As stated in §2*c*, at the beginning of the twentieth century, it probably exceeded ±25 m, whereas it is currently expected to be below ±10 m. Similarly, uncertainty on snow volume is high. Typically, it can be as high as 100%.

### (b) Path profile

For each site, the path profile was computed from a digital map. We needed to represent the path profile as a continuous curve. Here, we used the best polynomial approximations of the discretized path profile; that is, the polynomial curve that passes through the points: denotes the *N*th order polynomial approximation of the path profile, where *f*_{i}(1≤*i*≤*n*) represent the *n* discretization points of the path profile at a given set of path points *ξ*_{i}. Using orthogonality properties of Legendre polynomials, one writes the polynomial as where *P*_{k}(*ξ*)=*γ*_{k}d^{k}(*r*^{2}−1)^{k}/d*ξ*^{k} is the normalized Legendre polynomial of order *k*, defined over the domain [−1,1]; is the normalizing factor; and *r* is the scaled variable defined by *r*(*ξ*)=(2*ξ*−(*ξ*_{1}+*ξ*_{n}))/(*ξ*_{n}−*ξ*_{1}). The coefficients *α*_{k} are determined using the least-square method, that is, by minimizing the functional , where is a discretized regularization operator and the Lagrange multiplier *λ* must be tuned such that there is a good compromise between agreement and smoothness. For most paths, it is possible to obtain accurate approximations of the path profile by setting *N*=8 and *λ*=10^{−4}.

## 4. Results

### (a) Probability distribution function of *μ*

For each event, we computed the friction coefficient *μ* by numerically solving the equation *F*(*μ*)=*y*^{obs}. We then gathered the values related to the same path to compute the empirical probability distribution using the method described in §2*b*. This procedure is purely deterministic and does not take into account the uncertainty of observed run-out distances. To evaluate the sensitivity of our results to measurement errors, we randomly perturbed the recorded values and determined the *μ* values related to these perturbed samples. We found that there was no substantial difference in the curve shape, which led us to conclude that there is a certain structural stability of the *μ* distribution function against error measurement.

Figure 3 shows the empirical distribution functions of the *μ* sample for each path (thin lines) and for the total sample (thick line). The curves are bell-shaped with a peak at *μ*_{c}≈0.5. This bell shape is not surprising because the observations were censored at large *μ* values and there are little data at low *μ* values (related to extreme and rare avalanches). However, note that for other frictional laws, the probability distribution of the friction coefficient can be significantly modified; for instance, Ancey *et al*. (2003) studied the inverse problem for the Voellmy-like model, in which the frictional force is expressed in the form *μmg*+*Au*^{2} (Coulomb frictional contribution+turbulent-like term, with *A* another frictional coefficient) and did not find a smooth distribution for *μ* given that *μ* behaved as a discrete random rather than a continuous variable.

A closer look at figure 3 shows that there is a substantial variation in the peak position for paths taken individually: *μ*_{c} lies in the range of 0.46–0.60. Because it is unlikely that measurement errors can induce such a deviation between the distribution functions, we conclude that no universal probability distribution exists for the friction parameter *μ*. Physically, this shortcoming is expected because, clearly, different parameters such as snow consistency and avalanche volume can significantly influence avalanche mobility.

### (b) Looking for an average dependence of *μ* on *V*

Based on the preliminary analysis (see figure 1) and the results of the previous section, we will now investigate the *V* dependence of *μ* by using the methodology described in §2*c*. We have a set of data , where *y*_{ki} means the *i* th recorded elevation of path *k*(1≤*k*≤7). For each path, we computed a relationship linking the run-out distance and the friction coefficient by solving equation (2.3): *y*_{stop}=*F*^{(k)}(*μ*). We then tried to adjust a parametric function *μ*=*G*(*V*|*Θ*) such that the deviations between the observed and recorded run-out elevations are minimized. We were also interested in determining the uncertainty on the parameters *Θ* as well as the model precision; that is, the s.d. *σ* of the mismatch . We tested a two-parameter function (exponential model *G*_{2}(*V*)=*a*e^{−nV}). A four-parameter model, the power-law model *G*_{4}(*V*)=*a*(*b*+(*v*/*c*)^{n})^{−1}, was also tested, but we found that increasing the number of parameters did not enhance the model precision.

Table 2 provides the estimates of the best-fit parameters of *Θ*=(*a*, *n*) for the exponential law and indicates the model precision *σ* when all the data are considered; individual results for each path are also provided. Figure 4 displays the histograms of *a*, *n* and *σ* built from the values obtained by the MCMC simulations. Apart from *n*, the distribution of which is bimodal and non-symmetric, parameter *a* and precision *σ* can be described by a normal distribution. Note that the model accuracy is low; typically, if we know the avalanche volume, then we can predict the run-out elevation to 180 m. In comparison, given the measurement uncertainty, the worst estimate that we expect is: *ϵ*∼|*δ*_{y}|+|*F*′*δ*_{μ}|, where is the *μ* error related to the volume uncertainty *δV* and *F*′ is the sensitivity of *y*_{stop} to *μ* variations. Typically, for a middle-size avalanche with *V*=10^{5} m^{3}, we have the following estimates: |*δy*|=O(30) m (in the range 20–50 m), |*δv*|=5×10^{4} m^{3} (§3*a*), *F*′=O(2000) m (range 1000–5000 m); we deduce: |*δμ*|=O(2×10^{−2}); and *ϵ*∼70 m. Thus, the model uncertainty substantially exceeds the expected uncertainty due to measurement errors when all the data are considered. Better results are achieved for Haute-Savoie paths (Scey, Entremene, Nants and Grand Chantey), for which *σ* lies in the range of 41–83 m (see table 2).

This shortcoming can be made clearer by plotting *μ* as a function of the avalanche volume (see figure 5). The points were computed by assuming that the data * d* were exact and by inverting

*μ*=

*F*

^{−1}(

*y*

_{stop}). We have superimposed the two parametric curves

*G*

_{2}and

*G*

_{4}on the same plot and have indicated the 95% confidence interval. The confidence interval was built from the MCMC value sample. The points are widely scattered and this scatter is much more substantial than the measurement uncertainty. This explains why the residual error of the model is so significant (nearly three times larger than the uncertainty resulting from the measurement error). This prompts us to think that additional parameters must be taken into account in the

*V*dependence of

*μ*. The data nevertheless define a trend which could be approximated by eye as a plateau phase for volumes in the range of 10

^{3}–10

^{5}m

^{3}(

*μ*values close to 0.5 on average) and a slow decreasing phase for

*V*>10

^{5}m

^{3}.

In figure 5, in spite of the data scatter, the data provided by Ancey & Meunier (2004) lie in the upper part of the plot. Typically, when *V*→0, the *μ* value extrapolated from their data is approximately 0.68 versus 0.5 for the French database. A reasonable explanation is that Ancey & Meunier (2004) reported values related to monitored avalanches triggered under natural conditions or using explosives, whereas the avalanche database includes censored observations of naturally released avalanches, which leads to discarding avalanches characterized by a high *μ* value.

### (c) Normal random-effects model for *μ*(*V*)

The conclusions of the last subsection prompt us to consider friction as a random process. Instead of studying a functional dependence *μ*=*G*(*V*), we are seeking the probability *P*[*μ*|*V*] of observing *μ* provided the avalanche volume *V* is given. To this end, we will use the normal random-effects model presented in §2*d*, i.e. we assume that , with *G*_{2}(*V*)=*a*e^{−nV}. This means that the dependence of *μ* on *V* holds only on average, and the deviation from the mean trend is normally distributed with variance *σ*^{2}. We are looking for *σ*^{2}, *a* and *n*. To use the linear model of §2*d*, we first need to transform the data: (*μ*_{i}, *V*_{i})→(*ν*_{i}=−log *μ*_{i}, *V*_{i}).

The Gibbs sampler technique described in §2*d* was used. The results of typical simulations are plotted in figure 6 in the form of histograms for *a*, *n* and *σ*. We have also reported the residual error histogram . Table 3 collects the parameter estimates for *a*, *n* and *σ* for the two samples: (i) data from the French database and (ii) all the data (combining the 173 events of the French database and the 11 events from figure 1).

Comparison with results reported in figure 4 shows few changes for parameter *a* (value close to 0.50), while a significant alteration in *n* appears (3.22×10^{−6} against 7.0×10^{−7} previously). Note the much smaller value of the s.d. for *a* and *n*, which indicates a greater accuracy in the normal random-effects model adjusted on (*ν*_{i}, *V*_{i}) as compared with the model used in the previous section. The main difference between the two models is that, in the former model, a normal distribution of the *y*_{stop}−*y*^{obs} was assumed, whereas in the latter model, the normal distribution assumption concerns the friction coefficient *μ*.

In the normal random-effects model, the deviation is assumed to be normally distributed. The panel on the bottom-right quarter of figure 6 shows a more complex empirical distribution, exhibiting two peaks asymmetrically positioned on each side of 0 and a flat asymmetric tail extending over the interval [−0.2,0.15]. The normal distribution is only a crude approximation of the empirical distribution.

Figure 7 displays the function *G*_{2} adjusted on either sample as well as the *μ* data. The curves corresponding to quantiles 0.05 and 0.95 have been superimposed. Compared with figure 5, the average trend seems to describe the observed behaviour reasonably well, in spite of the scatter. Although the residual error *ϵ* is not normally distributed, it is seen that almost all the data lie within the quantile bounds, indicating that approximating *P*[*ϵ*] as a normal distribution provides realistic results concerning the extreme values (i.e. values far away from the mean average).

The major drawback in using the model to describe the dependence *μ*(*V*) can be identified on figure 7. For very large avalanche volumes (*V*>5×10^{6} m^{3}), negative values of *μ* can be generated, which is meaningless from a physical viewpoint. Furthermore, for *V*>5×10^{7} m^{3}, the friction coefficient comes very close to zero, which also seems unrealistic. However, at present, there is no historical information in the Alps concerning huge avalanches with such a snow volume. A refinement of the model would be to include a volume-dependence variance, for instance , with *G*_{2}(*V*)=*a*e^{−nV} and *σ*(*V*)=*a*′e^{−n′V}, where *a*′ and *n*′ are two additional parameters to be determined. Another criticism is that both intrinsic friction variability and model error are included in the probability distribution *P*[*μ*|*V*].

## 5. Concluding remarks

Different authors have put forward the idea that snow avalanches and debris flows might behave as granular flows and thus be described as Coulomb fluids. Such an idea has recently received partial support from a back analysis performed by Ancey & Meunier (2004) on well-documented snow avalanches and the intermediate-scale outdoor experiments on debris flows performed by Denlinger & Iverson (2001). This work offers interesting perspectives, especially concerning the use of sophisticated and robust theories such as the Savage–Hutter model to predetermine the behaviour of large avalanches. For this purpose, a better knowledge of the friction coefficient *μ* is required.

The objective of this paper was to calibrate the friction coefficient *μ* for snow avalanches. Because *μ* cannot generally be directly measured for real avalanches, an inverse problem must be solved, which involves determining *μ* for the Coulomb model outputs to match with field data. A major difficulty was that most field data available from different sources (in Europe) are inaccurate, making any attempt at fitting *μ* by using common adjustment methods difficult. Therefore, it was crucial to address the issue of data uncertainty. For the calibration work to be meaningful, it is also necessary to test the universality of the calibrated values; if it turns out that *μ* values significantly differ among avalanche paths, we can assume that the adjusted friction coefficient includes other effects than friction. This then leads to the question of the universality of the calibrated values and the interest of the Coulomb model to predetermine avalanche behaviour for any site.

The first hypothesis tested in this paper was to consider *μ* as an independent random variable. A sample of 173 events was used. It was found that the overall empirical distribution function of this sample was smooth and bell-shaped; its shape showed a slight dependence on the avalanche path. However, the variations in the distribution function of *μ* remained larger than data uncertainty for the Coulomb model. This suggests that *μ* depends on other parameters such as avalanche volume and snow property.

The second hypothesis was the possibility of a one-to-one correspondence between friction and avalanche volume. Adjusting different parametric functions *μ*=*G*(*V*) produced residual errors between the computed and observed run-out elevations that were, on average, three times as large as the uncertainty on field data. Therefore, additional variables should affect this relationship. We then opted for a probabilistic description of *μ* for a given avalanche volume, and assumed that the probability of observing *μ* conditionally to *V* was in the form .

It is too early to offer predictions about the universality of such a description, notably because the data that were used in this investigation were censored (avalanches with short run-out distance discarded) and the sample size was small. However, its relevance for engineering purposes is certain. A number of practical problems involve determining the 100 year old avalanche. If the period of return is defined in terms of the probability of observing an avalanche volume, then it is possible to obtain the bounds within which the run-out distance must lie for any long return-period avalanches. Indeed, the avalanche volume can be estimated, at least roughly, from snowfalls, whose probability distribution can be approximated by a generalized extreme value distribution (Coles 2001). Thus, it is possible to estimate the snow volume involved on average in an extreme avalanche. Then, by drawing *μ* samples from its probability distribution, one deduces the run-out distance, its expected mean value and the confidence interval.

## Acknowledgments

The author thanks the consulting group Toraval (France) for providing avalanche data, digital maps, and further comments on datum quality.

## Footnotes

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

- © 2005 The Royal Society