Bayesian calibration of an avalanche model from autocorrelated measurements along the flow: application to velocities extracted from photogrammetric images

Abstract

Physically-based avalanche propagation models must still be locally calibrated to provide robust predictions, e.g. in long-term forecasting and subsequent risk assessment. Friction parameters cannot be measured directly and need to be estimated from observations. Rich and diverse data are now increasingly available from test-sites, but for measurements made along flow propagation, potential autocorrelation should be explicitly accounted for. To this aim, this work proposes a comprehensive Bayesian calibration and statistical model selection framework. As a proof of concept, the framework was applied to an avalanche sliding block model with the standard Voellmy friction law and high rate photogrammetric images. An avalanche released at the Lautaret test-site and a synthetic data set based on the avalanche are used to test the approach and to illustrate its benefits. Results demonstrate (1) the efficiency of the proposed calibration scheme, and (2) that including autocorrelation in the statistical modelling definitely improves the accuracy of both parameter estimation and velocity predictions. Our approach could be extended without loss of generality to the calibration of any avalanche dynamics model from any type of measurement stemming from the same avalanche flow.

Introduction

The complexity of snow avalanche physics is related to the variability and changing nature of snow (Schweizer and others, 2003; Ancey, 2006; Castebrunet and others, 2012; Steinkogler and others, 2014). Evidence obtained at full-scale experimental slopes (Sovilla and others, 2008; Vriend and others, 2013; Prokop and others, 2015; Faug and others, 2018) shows a myriad of avalanche propagation and stopping regimes (Köhler and others, 2016), and numerical propagation models can reproduce these observations with increasing realism (Bartelt and others, 2016; Gaume and others, 2018). However, knowledge concerning the mechanical behaviour of snow during motion and associated processes (granulation, erosion/deposition, etc.) remains incomplete (Steinkogler and others, 2015; Truong and others, 2018). From a macroscopic point of view, experimental approaches (Casassa and others, 1989; Rognon and others, 2008; Kern and others, 2009) and the proposal of Voellmy (1964) suggest rheological behaviours which remain ad hoc. This renders on-site calibration on the basis of local data unavoidable (Ancey and Meunier, 2004; Salm, 2004; Eckert and others, 2012) to, e.g. predict high-return-period avalanches in land use planning and assess the related risk (Keylock and others, 1999; Meunier and others, 2004; Favier and others, 2014b,a). This is all the more true given that studies have shown that avalanche propagation models are highly sensitive to their friction parameter values (see e.g. Borstad and McClung, 2009; Fischer, 2013).

After deterministic inversion methods had shown their limits (Dent and Lang, 1980; Dent and others, 1998; Ancey and others, 2003), and following progress made in many fields where accurate numerical model calibration is now recognized as a crucial issue (e.g. Oakley and O'Hagan, 2002; Carmassi and others, 2018), the Bayesian framework has become an appealing avenue in snow science over the past years, especially in the frequent case of small data samples (e.g. Ancey, 2005; Straub and Grêt-Regamey, 2006; Eckert and others, 2007, 2008, 2009, 2010; Schläppy and others, 2014). Specifically, Gauer (2009) used a Bayesian framework to calibrate the friction parameters of three avalanche sliding block models and Fischer and others (2014) proposed a method to evaluate simulations compared to Doppler radar observations. Ultimately, Naaim and others (2013) could establish empirical links between friction parameters and physical properties of snow.

However, most of these existing approaches remain limited to rather coarse field data (e.g. samples of runout distances supplemented by input conditions), and when more comprehensive data sets have been considered, improper likelihood formulations have been used more often than not (Fischer and others, 2015). For instance, little attention has been given so far to the specific difficulty induced by potential autocorrelation between the data used for calibration. Whereas this is not a matter, for example, a sample of runout corresponding to distinct avalanche events, the assumption of independent observations is much more questionable in the current context of the increasingly diverse and rich measurements made on test-sites within the same avalanche. Many environmental applications have indeed demonstrated that neglecting potential autocorrelation between different measurements used in a calibration scheme can lead to biases in parameter estimation and/or to lower predictive performances (see e.g. Kuczera, 1983; McInerney and others, 2017; Schaefli and Kavetski, 2017; Sun and others, 2017). If this happens, both the physical interpretation of parameter estimates and the operational use of model predictions can be questioned.

On this basis, in this paper, we propose a Bayesian approach to calibrate the friction parameters of an avalanche propagation model using data of high temporal resolution (of the order of 1 s). Inspired by studies done for the calibration of hydrological models (e.g. Kuczera and Parent, 1998; Evin and others, 2014), the main novelty of our work is to explicitly account for potential autocorrelation between measurements made along the avalanche flow within the calibration framework. In what follows, we demonstrate that doing so results in a model better supported by the data that improves the accuracy of friction parameter estimates and of velocity predictions. Application is made on a well-documented avalanche event from the Lautaret full-scale test-site (Thibert and others, 2015). The used velocity data were obtained from high rate positioning from photogrammetric images. Such data were already applied by various authors to avalanche simulations performed with reference friction parameter values (e.g. Turnbull and Bartelt, 2003; Gauer, 2014; Dreier and others, 2016), but without including them within an explicit calibration scheme. That taking into account potential autocorrelation between measurements within the calibration is unavoidable to get unbiased estimates is further demonstrated with synthetic data analogue to the case-study.

In this work, as a proof of concept, we use the sliding block avalanche model, also known as 1-D Voellmy model instead of state of the art depth-averaged models (Naaim and others, 2004; Christen and others, 2010; Bartelt and others, 2012). We are aware of its limitations, notably that it makes it impossible to depict, e.g. flow depth variations in space and time. Also, some authors have shown that it may underestimate avalanche velocities (see e.g. Ancey and Meunier, 2004; Gauer, 2014). However, for hazard mitigation, simple models with few parameters remain useful (Salm, 2004) and have the advantage to allow fast computation in comparison to more complex ones. A very simple avalanche propagation model is therefore a good choice for developing a calibration approach which could be in the future applied to any other, more advanced, avalanche propagation model as soon as autocorrelation in measurements series is suspected.

Avalanche model calibration principle

Sliding block propagation model

Our model considers the avalanche as a rigid body sliding over a bidimensional curvilinear profile starting from the top of the path. The mass m and body shape variations of the avalanche are neglected. Under these assumptions, the motion equation of the avalanche mass centre is:

(1)$${{\rm{d}} u\over {\rm d} t} = g \sin \theta - {{F}\over m}\comma\; $$

where $u=\Vert \vec {u}\Vert$, du/dt is the acceleration, g is the gravity constant, θ is the local slope angle and ${F}=\Vert \vec {F}\Vert$ is the frictional force. In this study, we consider the classical Voellmy friction law (Voellmy, 1964). This means that the friction force is:

(2)$${F}=\mu m g \cos \theta + {mg\over \xi h} u^{2}\comma\; $$

depending on two friction parameters Θ = {μ, ξ}. Often it is assumed that μ evolves with the physical properties of snow whereas ξ may correspond to the geometry of the avalanche path and to terrain roughness (Ancey and others, 2003). However, whether this interpretation is sound or not is not our debate here. We do not make any further assumption regarding the linkages between (μ, ξ) and snow and topographical variables and simply search for the best couple on the basis of the data. The propagation model also depends on three forcing variables x = {T, h, x _start} where T defines the topography of the terrain, h the mean flow depth of the avalanche and x _start the release abscissa (the protected runout length) of the avalanche mass centre(Ancey and Meunier, 2004; Eckert and others, 2007).

Statistical model formulation

Let us denote f our avalanche propagation model. The model f predicts the avalanche speed (m · s⁻¹) and position along the slope (m) at a time t (s). The model depends on parameters $\Theta \in {\opf R}^{p}$ and forcing variables $x \in {\opf R}^{d}$. The observed velocity, collected on the field, at time t is noted by v _t and we denote v _obs = {v ₁, …v _n} the set of observations where n denotes the number of observations.

The aim of model calibration is (1) to find the optimal combination of parameters Θ that minimizes the discrepancy between the observations v _obs and the model simulation f(Θ, x) and (2) to rigorously quantify the associated uncertainty. To this end, we use the generic statistical model:

(3)$${\cal M}\colon v_{t} = f_{t}\lpar x\comma\; \Theta\rpar + \epsilon_{t}\comma\; \quad \forall t \in \lcub 1\comma\; \ldots n\rcub \comma\; $$

where f _t(Θ, x) denotes the simulation of the avalanche velocity at time t and ε _t is the model error. With this classical additive formulation, propagation model errors and observation errors are modelled altogether in the residuals ε _t.

In nearly all existing avalanche model calibration approaches, model errors are implicitly or explicitly assumed as independent and identically normally distributed (iid). In other words, $\epsilon _{t} \mathop{_{\sim}^{\rm iid}} {\cal N}\lpar 0\comma\; \sigma ^2\rpar$, $\forall t \in \lcub 1\comma\; \ldots n\rcub$, where ${\cal N}$ denotes the Normal distribution and σ² is the common variance of the errors. However, this may be a too strong assumption for different measurements made along the same avalanche. For instance, the errors ε _t are likely to present a non-negligible correlation between two consecutive time steps (in other words, ε _t and ε _t−1 may be correlated). To include the errors' autocorrelation in the calibration we propose to model them as an autoregressive (AR) process. Specifically, we consider an AR model of order 1 only (AR1). AR models of higher order could be used but this would imply the estimation of additional parameters, a non-trivial task in our case of a limited velocity sample. In addition, results obtained for the application indeed suggest that an AR1 is sufficient to accurately account for the data variability (see Application, results and discussion section). Thus, the errors ε _t can be expressed as:

(4)$$\epsilon_{t}=\phi \epsilon_{t-1} + \eta_{t}\comma\; \quad \eta_{t} \mathop{_{\sim}^{\rm iid}} {\cal N}\lpar 0\comma\; \sigma^2\rpar \ \forall t \in \lcub 1\comma\; \ldots n\rcub \comma\; $$

where $\phi \in {\opf R}$ and η_t are the coefficient and innovations of the AR1 process, respectively. If |ϕ| < 1, $\lpar \epsilon _t\rpar _{t \in {\opf N}}$ is defined as the unique stationary solution of Eqn (4). In our study, |ϕ| < 1 (see Application, results and discussion section) that guarantees the stationarity of $\lpar \epsilon _t\rpar _{t \in {\opf Z}}$. Note also that with this model the innovations η_t are normally independently distributed but not the errors ε _t.

Hereafter, the model with the assumption of normally distributed and independent errors is denoted ${\cal M}_{0}$ and the model with AR1 errors is denoted ${\cal M}_{1}$. Depending on the ${\cal M}_{i}$ model with i ∈ {0, 1}, the whole set of parameters to be estimated is different. Table 1 summarizes the parametrization of the competing statistical models considered. Hence, model ${\cal M}_{1}$ with ϕ = 0 corresponds to ${\cal M}_{0}$.

Table 1. Considered statistical models and corresponding parameters

Bayesian framework

The probability of the data can be maximized with respect to the model parameters (Fisher, 1922). Bayesian statistic is an useful framework to estimate model parameters from scarce data. Within this approach, the quantification of uncertainty in parameter estimation is straightforward. In fact, the advantage of the Bayesian approach is that the uncertainty on parameters is assessed through credibility intervals by contrast to traditional methods (confidence intervals) (Bayes, 1763; Bernardo and Smith, 2009). Hence, Bayesian statistics is now widely accepted as a reasonable option in environmental sciences (Berliner, 2003; Clark, 2005) and we use this framework in what follows.

For simplicity, let us denote γ, the set of parameters related to the errors, it means γ = {σ²} for ${\cal M}_{0}$ and γ = {ϕ, σ²} for ${\cal M}_{1}$, and x = {T, h, x _start} the forcing variables. Under a Bayesian framework, the joint posterior distribution of the parameters is the following one:

(5)$$\pi\lpar \Theta\comma\; \gamma \vert v_{{\rm {obs}}}\comma\; x\rpar \propto {\cal L}\lpar v_{{\rm {obs}}}\vert \Theta\comma\; \gamma\comma\; x\rpar \pi\lpar \Theta\comma\; \gamma\rpar \comma\; $$

where ${\cal L}\lpar v_{{\rm {obs}}}\vert \Theta \comma\; \gamma \comma\; x\rpar$ is the likelihood of the observations and π(Θ, γ) is the joint prior distribution. Thus, the likelihood ${\cal L}\lpar v_{{\rm {obs}}}\vert \Theta \comma\; \gamma \comma\; x\rpar$ is required. Under the ${\cal M}_{0}$ model, the observations follow a Gaussian distribution and the likelihood ${\cal L}\lpar v_{{\rm {obs}}}\vert \Theta \comma\; \gamma \comma\; x\comma\; {\cal M}_{0}\rpar$ writes as:

(6)$${1\over \lpar 2\pi \sigma^2\rpar ^{n/2}} \exp \left[-{1\over 2\sigma^2} \mathop{\sum}\limits_{t=1}^{n}\lpar v_{t}-f_{t}\lpar \Theta\comma\; x\rpar \rpar ^2 \right].$$

Under the ${\cal M}_1$ model, the likelihood ${\cal L}\lpar v_{{\rm {obs}}}\vert \Theta \comma\; \gamma \comma\; x\comma\; {\cal M}_{1}\rpar$ writes as (for more detail, see Appendix A):

(7)$$\eqalignno{& \sqrt{{1-\phi^2\over 2 \pi \sigma^2}} \exp \left[-{1-\phi^2\over 2\sigma^2} \lpar v_{1}-f_{1}\lpar \Theta \comma\; x\rpar \rpar ^{2} \right]\cr & \quad \times {1\over \lpar 2\pi \sigma^2\rpar ^{\lpar n-1\rpar /2}} \exp \left[-{1\over 2\sigma^2} \mathop{\sum}\limits_{t=2}^{n} \eta_{t}^2 \right].}$$

Note that, in this study, we do not calibrate the quantities x _start, h, T because they were inferred from the data to put the effort on the friction law calibration. The sensitivity to these quantities in avalanche models has been studied in the works of e.g. Barbolini and Savi (2001), Borstad and McClung (2009), Buhler and others (2018). The authors found that changes in avalanche volume have a larger effect on both runout distance and avalanche velocity and that the friction coefficient μ (in a Coulomb model) is of high importance. However, this question should be analysed deeper in a formal statistical framework and it is out of the scope of this work.

Metropolis–Hastings algorithm

The main practical difficulty in Bayesian inference is how to compute the normalizing constant in Bayes theorem. We overcome it by implementing a sequential Metropolis–Hastings algorithm, hereafter denoted MH (Metropolis and others, 1953; Torre and others, 2001). The MH algorithm proposes a generic way to construct a stationary and ergodic Markov chain that converges, under mild conditions, to the posterior distribution π(Θ, γ|v _obs, x). The Markov chain returned by the algorithm can be considered as a sample from π(Θ, γ|v _obs, x).

The following description of the MH algorithm is obtained from Robert (2015). Let us denote, for simplicity, ψ = {Θ, γ} the set of the error and model parameters. For the application of the MH algorithm, it is needed an initial value ψ ⁽⁰⁾ and a proposal distribution q. Each iteration k of the algorithm consists in:

1. 1. Generating ψ′ ~ q(.|ψ ^(k−1)).

2. 2. Calculating $u \sim {\cal U}\lpar 0\comma\; 1\rpar$.

3. 3. Taking

   (8)$$\psi^{\lpar k\rpar }= \left\{\matrix{ \!\!\!\!\!\!\!\! \psi' & {\rm{if}}\ {\it u} \leq \alpha\comma\; \cr \psi^{\lpar k-1\rpar } & {\rm{otherwise}} }\right. $$
   with
   (9)$$\alpha = \min \left({{\cal L}\lpar v_{{\rm {obs}}}\vert \psi'\comma\; x\rpar \pi\lpar \psi'\rpar \over {\cal L}\lpar v_{{\rm {obs}}}\vert \psi^{\lpar k-1\rpar }\comma\; x\rpar \pi\lpar \psi^{\lpar k-1\rpar }\rpar } {q\lpar \psi^{\lpar k-1\rpar }\vert \psi'\rpar \over q\lpar {\psi'\vert \psi^{\lpar k-1\rpar }}\rpar }\comma\; 1 \right)\comma\; $$

where we recall that ${\cal L}\lpar v_{{\rm {obs}}}\vert \psi \comma\; x\rpar$ and π(ψ) = π(Θ, γ) stand for the likelihood of velocity observations and the joint prior distribution, respectively. As mentioned in Robert (2015), the performance of the algorithm depends on the choice of q. For example, a standard choice for the proposal distribution q is a multinormal distribution centred in ψ ^(k−1) and with a given covariance Σ_q, which defines a random walk. For our application, Σ_q was tuned according to the optimal acceptance rates that grant fast convergence of the MH algorithm (Gelman and Rubin, 1992).

Model selection using Bayes Factor

To compare the accuracy of the two competing models ${\cal M}_0$ and ${\cal M}_1$, we use the Bayes factor. The Bayes factor is a criterion of the evidence provided by the data (in our case v _obs) to reject a model (${\cal M}_{0}$) compared to another one (${\cal M}_{1}$). The Bayes factor is defined as the ratio of marginal probabilities (see more details in Kass and Raftery, 1995):

(10)$$B_{10} = {{\,p}\lpar v_{{\rm {obs}}}\vert {\cal M}_{1}\comma\; x\rpar \over {\,p}\lpar v_{{\rm {obs}}}\vert {\cal M}_{0}\comma\; x\rpar }\comma\; $$

where

(11)$${\,p}\lpar v_{{\rm {obs}}}\vert {\cal M}_{i}\comma\; x\rpar = \int {\cal L}\lpar v_{{\rm {obs}}} \vert \psi_{i}\comma\; {\cal M}_{i}\comma\; x\rpar \pi \lpar \psi_{i}\vert {\cal M}_{i}\comma\; x\rpar \, {\rm d} \psi_{i}\comma\; \quad i \in \lcub 0\comma\; 1\rcub. $$

The conditioning by ${\cal M}_{i}$ highlights that likelihood, prior and posterior distributions depend on the considered model. The Bayes factor can be estimated by applying importance sampling and using the MH sample drawn from the posterior density $\pi \lpar \psi _{i}\vert v_{{\rm {obs}}}\comma\; x\comma\; {\cal M}_{i}\rpar$ (see more details in Appendix B). An interpretation of the numerical value obtained to determine if there is evidence provided by the data to reject the model ${\cal M}_0$ compared to ${\cal M}_1$ according to the value of log₁₀B ₁₀ is proposed in the work of Kass and Raftery (1995), where log₁₀ denotes the common logarithm. Notably, a log₁₀B ₁₀ value between 1 and 2 indicates strong evidence in favour of rejecting ${\cal M}_0$ compared to ${\cal M}_1$, and log₁₀B ₁₀ > 2 a decisive evidence in favour of ${\cal M}_1$.

Application data

Avalanche data from the Lautaret test-site

The case study is an avalanche released at the Lautaret full-scale test-site (Thibert and others, 2015). This test-site, located in the French Alps, holds a succession of avalanche paths. Here, the path referred as ‘path number 2’ was used to artificially trigger an avalanche on 13 February 2013 (Fig. 1). This path is 450 m long, dropping from 2400 to 2100 m a.s.l. (Fig. 2). Its upper part, where the acceleration of the flow occurs, is steep with an average inclination of 37° and a maximal slope around 45° in the starting area. This part of the path is steep-sided and around 10 m-wide, so that the flow is channelized. The lower part of the path is mainly the run-out zone, a large and naturally open slope. At the transition between the two parts, a road (Col du Galibier road, open in summer only) crosses the avalanche path, which generates a local slope rupture in the avalanche track.

Fig. 1. Three-dimensional topography of path number 2 of the Lautaret test-site, some front positions (black lines) of the avalanche released on 13 February 2013 as determined from photogrammetric measurements and changes in snow depth before/after the avalanche as inferred from terrestrial laser scanning.

Fig. 2. Two-dimensional topography of the considered avalanche path (path number 2 of the Lautaret test-site), and position of the avalanche mass centre each second for the studied avalanche released on 13 February 2013 (red circles).

The avalanche released on 13 February 2013 was composed of a dense part with a limited saltation layer on top. Properties of the snow involved in the flow were characterized with a density, temperature and hardness profile of the released snow layer. Snow grain types and dimensions were also characterized. The avalanche was released at 11.58 hours when air temperature was about −10°C. A 0.25 m thick layer of fragmented and decomposing snow particles was released. The mean density was 250 kg · m⁻³, ranging between 270 and 225 kg · m⁻³ from the bottom to the surface of the released snow cover. In the released layer, particles diameter is less than 0.5 mm. Snow temperature was between −4.7 and −5°C. Hardness was ‘fist’ (hand index) and measured as 20 N in Ram Resistance Equivalents (Fierz and others, 2009).

Avalanche front positions were determined using a high rate photogrammetric system that was specifically developed (Soruco and others, 2011; Thibert and others, 2015). We used a low-cost non-metric imagery system (numerical reflex Nikon D2Xs cameras in DX format, CMOS sensor with Nikon 85 mm f/1.4 AF fixed focal lenses). An advanced ad hoc tuning was performed to account for the radial distortion of the lenses, the decentration of the principal point (principal point shift) and the exact focal length of the lenses required for the correct scaling of the images (Faig and others, 1990). The resulting error in positioning avalanche fronts was estimated to be less than 25 cm after image orientation on ground control points and a comparison of direct photogrammetric measurements to laser scanning on a test area. Synchronization between the two cameras was achieved within a precision of 6 · 10⁻⁶ s, therefore the error associated with the time sequence is also negligible. Eventually, terrestrial laser scanning was used to retrieve terrain elevation before and after avalanche triggering (Prokop, 2008; Prokop and others, 2015), so as to quantify the snow mass transfer by the avalanche (Fig. 1). In Appendix D, Figure 10 shows some examples of images used in the photogrammetric process.

The sliding block model is a one-dimensional model representing the successive positions of the centre of mass of the avalanche, thus we used the following procedure to determine the velocity vector v _obs. At each time t, we computed the centre of mass position of the avalanche using a set of front points ${\cal D}_{t} \in {\opf R}^3$ from the delineated front. ${\cal D}_{t}$ excludes the lateral front points because they are less active in the avalanche flow (Pulfer and others, 2013):

(12)$$c_{t} = {\sum_{\,p_{t} \in {\cal D}_{t}} p_{t}\over n_{t}}\comma\; $$

where n _t is the number of elements in ${\cal D}_{t}$. Note that, from the data available, it is not possible to calculate with certainty the position of the mass centre of each avalanche front. In this study, we estimated each of them using eqn (12), arguably a rough estimation but sufficient for obtaining an approximation of the avalanche velocity consistent with the modelling framework we are using. Figure 1 shows the front of the avalanche and the centre of mass at each 5 s. Then, the velocity at time t is estimated as the mean velocity of the centre of mass between two successive images:

(13)$$v_{t}={c_{t}-c_{t-1}\over \Delta t}.$$

In our application, Δt = 1s and the data set is composed of 21 observations. We consider the first 21 s before the avalanche splits in sub-avalanches to ensure the validity of the application of a sliding block model. The mean flow depth h value was calculated as the mean of the difference between the snow depth of the digital elevation model taken after the release of the avalanche and the snow depth registered at the mass centre locations during the avalanche motion. The value calculated was h = 2.19 m with a standard deviation of 0.78 m. This estimation is rough but it is a reasonable value given the characteristics of the avalanche studied.

The topography T was constructed as follows: the site under study has a simple geometry. It is a rather straight avalanche path starting with a cornice and of limited and rather constant width. The chosen 2-D topography is the main flow path starting from the cornice in the middle of the path (see more details in Thibert and others, 2015). From a fine-scale Digital Elevation Model, it was extracted as a grid of horizontal resolution of 1.4 m. This is largely small enough to make the impact of the numerical approximation on the computed velocities negligible. Indeed, Bühler and others (2011) showed that a spatial resolution of 25 m is sufficient for avalanche modelling.

Figure 2 shows the topography T and the points where the fronts were recorded.

Synthetic data generation

In order to further evaluate the accuracy of parameter estimation with the two models, we used synthetic data. This standard approach in statistical developments allows validation because ‘truth’ (i.e. the parameter values used for the synthetic data generation) is known. Any parameter values could have been used at this stage. To resemble the most possible to the measured avalanche velocities, we chose the combination of parameter estimates corresponding to the measured case study. Also, to highlight that taking into account autocorrelation between measurements within the calibration is crucial to get unbiased estimates when such autocorrelation actually exists, our synthetic data were generated with model ${\cal M}_1$.

Specifically, synthetic velocity time series were generated from the maximum a posteriori (MAP) resulting from the application of our calibration approach to the avalanche described above with model ${\cal M}_1$. The MAP is the mode of the posterior distribution, namely the most plausible value given the data (Table 2). In detail, after the MAP estimators of the μ, ξ, σ² and ϕ parameters of ${\cal M}_1$ model were determined for the avalanche, we proceed as follows:

1. 1. An avalanche model simulation was conducted using the μ _MAP and ξ _MAP parameters. The result of this simulation is a velocity time series denoted by f(Θ_MAP, x).

2. 2. An AR(1) error of parameters $\sigma ^2_{\rm {MAP}}$ and ϕ _MAP was simulated.

3. 3. The AR(1) error simulated was added to f(Θ_MAP, x). If negative velocity was obtained, the AR(1) error was resampled.

Table 2. Parameter values used for the generation of the synthetic data set. They correspond to the maximum (MAP) obtained with model ${\cal M}_1$ for the avalanche of the 13 February 2013

We generated 100 samples of synthetic data following this procedure and, for each of these, the 21 velocity values which correspond to the same location of the measurement were kept (to stay with the same data size as for the avalanche). Figure 3 shows some examples of the velocity samples generated versus the avalanche data.

Fig. 3. Samples of the generated synthetic longitudinal velocity profiles (colour lines) and observations corresponding to the studied avalanche (dots).

Application, results and discussion

The avalanche

Prior distributions

Our Bayesian framework was used to calibrate the parameters of the ${\cal M}_{0}$ and ${\cal M}_{1}$ models (see Table 1). Marginal prior distributions used are described in Table 3. Parameters were assumed marginally independent a priori. This assumption is not strong because the dependence of the parameters is reflected in their joint posterior distribution (Gilks and others, 1995; Eckert and others, 2010). Informative marginal priors were used for all parameters but one, ϕ, for which a vague (poorly informative) prior was used instead, namely an uniform distribution ${\cal U}\lpar -\!1\comma\; 1\rpar$. This latter choice was done because prior knowledge for ϕ was unavailable. Instead, informative marginal prior distributions for other parameters could be determined on the basis of expert knowledge and well known reference values from Salm and others (1990). In addition, a comprehensive prior sensitivity analysis was conducted (see Appendix C). It demonstrates that our results are highly robust to the prior choice so that this question is no longer further considered in what follows.

Table 3. Parameters marginal prior distributions. The prior distribution of ϕ applies only to ${\cal M}_1$ model. Γ represents the Gamma distribution and ${\cal U}$ the uniform distribution

Convergence of the MH algorithm

For both models, we generated 50 000 MH samples and the last 25 000 were kept. The first 25 000 iterations were discarded as a burn-in period. This step avoids dependence on initial values. To further assess the convergence of the MH algorithm, we computed standard diagnoses (Kuczera and Parent, 1998; Torre and others, 2001; Robert and Casella, 2009; Eckert and others, 2010). The R Core Team (2019) package coda created by Plummer and others (2006) was used to calculate some of these. Specifically, three parallel chains were generated starting from different initial points and it was graphically checked that they were mixing well enough (not displayed). Also, the Kolmogorov–Smirnov test (H₀: the two samples were drawn from the same continuous distribution) was applied to the samples from the different chains and results showed that these were indeed drawn from the same distribution. In addition, for both models and all parameters, it was checked that acceptance rates (defined as the number of candidates accepted over the total number of iterations of the MH algorithm) were close enough to 0.25, the optimal value according to Robert (2015) reach fast convergence in a random walk. Finally, the Gelman–Rubin convergence criterion was computed (Gelman and Rubin, 1992). It is based on the difference between a weighted estimator of the variance and the variance of the estimators from the different chains (Robert and Casella, 2009). In our case, it equals 1 for all parameters, a perfect result. This all demonstrates that convergence was clearly reached for both models applied to the case study, leading to meaningful numerical approximations of the joint posterior distributions and, hence, reliable posterior estimates of the parameters of both models and of the Bayes factor.

Posterior distributions and posterior estimates

Figure 4 shows the parameter prior and posterior densities according to the MH samples. In Figures 4a–c (respectively, Figs 4d–g), the ${\cal M}_{0}$ densities are shown (respectively, ${\cal M}_{1}$) and Table 4 sums up the corresponding descriptive statistics. Eventually, posterior correlations between model parameters are shown in Figures 5 and 6 for ${\cal M}_0$ and ${\cal M}_1$, respectively. For the two models, the marginal posterior distributions have lower variance than their priors, meaning that the observations have conveyed information into the analysis (see Fig. 4), an expected result.

Fig. 4. Prior and posterior densities of model parameters. Panels (a–c) model ${\cal M}_0$ and panels (d–g) model ${\cal M}_1$.

Fig. 5. Joint posterior distribution of parameters of ${\cal M}_0$ model highlighting inter-parameter correlations.

Fig. 6. Joint posterior distribution of parameters of ${\cal M}_1$ model highlighting inter-parameter correlations.

Table 4. Posterior distribution characteristics: posterior mean, posterior standard deviation (sd.), median (50% percentile) and 95% credibility interval (2.5% and 97.5% percentiles).

The variances of the posterior distribution of the friction parameters μ and ξ in the ${\cal M}_0$ model are lower than in ${\cal M}_1$ model. This result (arguably the sole undesirable one with ${\cal M}_1$ compared to ${\cal M}_0$) could be a consequence of the interactions between the friction coefficients and error parameters in the more parametrized ${\cal M}_1$ model (four free parameters instead of three with ${\cal M}_0$). Indeed, there is strong correlation between the parameters ξ and ϕ (0.49) for the model ${\cal M}_1$, which may preclude reaching sharp estimates of friction parameters with model ${\cal M}_1$. Conversely, even if there is a high correlation between the parameters μ and ξ with both models, switching from ${\cal M}_0$ to ${\cal M}_1$ reduces it from 0.85 to 0.50. The high correlation between the two parameters of the Voellmy friction law is known since the first calibration approaches of, e.g. Dent and Lang (1980). This usually limits robust interpretation of obtained estimates, so that reducing it thanks to ${\cal M}_1$ should be seen as advantageous. This correlation reduction may be a collateral effect of having one more free parameter, allowing a bit more flexibility to fit the data. Similar effect is reflected in the lower correlation between σ² and the friction parameters with ${\cal M}_1$ than with ${\cal M}_0$. Note by the way that fairly assessing such correlations is a real strong point of our formal Bayesian calibration approach.

Remarkably, posterior estimates of friction parameters are very contrasted under both models. Even if posterior densities are not significantly different at the 95% credibility level (notably because of the higher a posteriori variance with model ${\cal M}_1$), differences in posterior estimates reaches 13% for μ (relative difference between 0.34 with ${\cal M}_1$ instead of 0.3 with ${\cal M}_0$), and 70% for ξ (relative difference between 696 m² · s⁻¹ with model ${\cal M}_1$ and 410 m² · s⁻¹ with model ${\cal M}_0$). This indicates that using either ${\cal M}_0$ or ${\cal M}_1$ makes a huge difference if one aims at interpreting the value of posterior estimates, for, e.g. relating avalanche friction characteristics to snow and topographical conditions.

Another remarkable result is the diminution of the error variance σ² from model ${\cal M}_{0}$ to model ${\cal M}_{1}$ (Table 4). Indeed, the errors variance σ² is higher if the autocorrelation ϕ is not included. Specifically, the standard deviation of the errors decreases from 3.6 to 3.1 m · s⁻¹ indicating that with model ${\cal M}_1$ velocity predictions may be seen as 13% more accurate (relative difference between both estimates), another desirable property for practical use in snow science. Eventually, the autocorrelation parameter ϕ is largely positive, with a posterior mean of 0.71. Also, its posterior $95\percnt$ credibility interval whose lower bound is 0.33 firmly excludes the zero value corresponding to ${\cal M}_0$. This pleads for a high and significant autocorrelation between velocities along flow propagation.

Model posterior estimates

We calculated for model ${\cal M}_{0}$ the posterior estimate of Voellmy model simulations f(Θ, x) from the posterior estimates of parameters μ and ξ (Table 4). Resulting posterior estimates of model errors ε _t were analysed by applying the Ljung–Box test (H₀: the data are independently distributed). We found a p-value lower than 0.05 indicating that there is a significant autocorrelation of model errors. In other words, cor(ε _t, ε _t−1) for all t > 1 is significant. Hence, one of the assumptions underlying statistical model ${\cal M}_{0}$ is not fulfilled by the data. On the other hand, applying standard statistical tests to posterior estimates of model errors shows no evidence of non-Gaussian errors (Shapiro test H₀: the sample is normally distributed, p-value > 0.05) or heteroscedasticity. The assumptions of Gaussian errors of common variance underlying statistical model ${\cal M}_{0}$ are thus fulfilled.

We proceeded similarly for statistical model ${\cal M}_{1}$, applying the same statistical tests to the posterior estimates of the innovations η_t of the AR(1) model. Remember that, in the ${\cal M}_{1}$ model, the independent assumption is on the innovations η _t and not on the errors ε _t. Test results show that innovations are indeed independent (p-value of the Ljung–Box > 0.05), and normally distributed (p-value of Shapiro test > 0.05) with common variance. These results show that ${\cal M}_1$ represents correctly the autocorrelation of the errors, and, more widely, that contrary to model ${\cal M}_{1}$ all of its underlying assumption are fulfilled, which promotes its use from a strict statistical point of view.

Predictive velocity distributions

To analyse how the statistical model choice affects velocity estimation, we propagated parameter uncertainty up to model predictions. Two sets of posterior predictive simulations were performed, the first integrating over the posterior distribution of friction parameters only, leading the posterior predictive distribution ${p}\lpar v_{x}\vert v_{{\rm {obs}}}\comma\; x\comma\; {\cal M}_i\rpar$ of the avalanche propagation model, and the second integrating over the distribution of both friction parameters and model error parameters leading the full posterior predictive distribution of avalanche velocities for the case study ${p}'\lpar v_{x}\vert v_{{\rm {obs}}}\comma\; x\comma\; {\cal M}_i\rpar$. More precisely, the first writes:

(14)$${\,p}\lpar v_{x}\vert v_{{\rm {obs}}}\comma\; x\comma\; {\cal M}_i\rpar =\int f\lpar \Theta\comma\; x\rpar {\,p}\lpar \Theta\vert v_{{\rm {obs}}}\comma\; x\comma\; {\cal M}_i\rpar \, {\rm d} \Theta\comma\; $$

and, the second one:

(15)$${\,p}'\lpar v_{x}\vert v_{{\rm {obs}}}\comma\; x\comma\; {\cal M}_i\rpar = \int \left(f\lpar \Theta\comma\; x\rpar +\epsilon \right){\,p}\lpar \Theta\comma\; \gamma\vert v_{{\rm {obs}}}\comma\; x\comma\; {\cal M}_i\rpar \, {\rm d} \Theta \, {\rm d}\gamma.$$

Eventually, to quantify the added value of parameter inference, the prior distribution of the velocity was also calculated by integrating over the prior distribution on the friction and model error parameters:

(16)$$\pi\lpar v_{x}\vert x\comma\; {\cal M}_i\rpar =\int \left(f\lpar \Theta\comma\; x\rpar + \epsilon \right)\pi\lpar \Theta\comma\; \gamma\vert {\cal M}_i\rpar \, {\rm d} \Theta \, {\rm d}\gamma.$$

By using a Bayesian approach, we obtained a sample of parameters (μ, ξ, σ ²) for ${\cal M}_0$ (and also of ϕ for ${\cal M}_1$). Then, we computed avalanche model simulations using these two samples to obtain several curves of velocity. Finally, we drawn the percentile of each curve in Figure 7. Also, this figure shows the resulting 90% credibility intervals. For both models, comparison between prior and posterior credible intervals shows that predicted velocities are logically shifted towards observations. Uncertainty reduction (change in the width of the 90% credible intervals) is stronger with ${\cal M}_{0}$, but also exists with ${\cal M}_{1}$. This simply reflects the larger a posterior variance of model ${\cal M}_{1}$ parameters highlighted above.

Fig. 7. Predictive velocity distributions versus data: (a) model ${\cal M}_0$ and (b) model ${\cal M}_1$. 90% posterior credibility intervals (CI) are computed according to Eqns (14) and (15) for the Voellmy propagation model (orange dotted lines), and the complete propagation and error model (blue dotted lines), respectively. Black plain line denotes the posterior median of the complete model. 90% prior credibility intervals are computed according to Eqn (16) (green dotted lines). Observations used for calibration appear as red points.

In detail, almost all the elements of v _obs, except two with model ${\cal M}_{0}$ and one with model ${\cal M}_{1}$, are inside the 90% posterior credibility intervals drawn on ${p}'\lpar v_{x}\vert v_{{\rm {obs}}}\comma\; x\comma\; {\cal M}_i\rpar$ (blue dotted lines Fig. 7), so that from this perspective the added value of ${\cal M}_{1}$ is not obvious. However, why ${\cal M}_{1}$ model is advantageous clearly appears if one focuses on the deterministic Voellmy model component (dotted orange lines in Fig. 7): predictions come closer to observations and notably the velocity underestimation generally attributed to the Voellmy sliding block model (see e.g. Ancey and Meunier, 2004; Gauer, 2014) is reduced. Indeed, only 33% of the observations are inside the 90% posterior credibility intervals of the sliding block model simulations ${p}\lpar v_{x}\vert v_{{\rm {obs}}}\comma\; x\comma\; {\cal M}_0\rpar$, compared to 62% with ${\cal M}_{1}$ model.

Another quantification is provided by the systematic difference between both models that we evaluated through the mean difference and the mean quadratic difference between the medians of predicted velocities under ${\cal M}_{0}$ (Fig. 7a, black line) and ${\cal M}_{1}$ (Fig. 7b, black line). If we denote $q50_{v_x\comma {\cal M}_i}$ the median (50% percentile) of the posterior predictive distribution of velocities under model ${\cal M}_i$, the first writes:

$$\Delta_1={1\over 240}\int_{0}^{240}\lpar q50_{v_x\comma {\cal M}_1}-q50_{v_x\comma {\cal M}_0}\rpar \, {\rm d}x$$

and the second:

$$\Delta_2=\sqrt{{1\over 240}\int_{0}^{240}\lpar q50_{v_x\comma {\cal M}_1}-q50_{v_x\comma {\cal M}_0}\rpar ^2 \, {\rm d}x}$$

Obtained values are Δ₁ = 1.31 m · s⁻¹ and Δ₂ = 1.72 m · s⁻¹. The positive value of the mean difference and the fact that it is not that much lower than the mean quadratic difference clearly shows that the inclusion of autocorrelation into the modelling truly leads to systematically higher velocity predictions. Predicted velocities are rather constant between x = 50 m and x = 250 m (plateau phase), around 15 m · s⁻¹ with ${\cal M}_{0}$. The underestimation with ${\cal M}_{0}$ with regards to ${\cal M}_{1}$ can be estimated by the ratio between Δ₂ and this velocity to $11.5\percnt$.

All in all, the analysis of predictive velocity distributions further confirms that model ${\cal M}_{1}$ should be preferred to get sharper, less underestimated, point estimates of avalanche velocities and to assess the related uncertainty fairly (credibility intervals larger but more likely to realistically describe the range of possible results).

Model selection

Finally, the logarithm of the Bayes Factor log₁₀ B₁₀ was calculated to compare the ${\cal M}_{0}$ model to the ${\cal M}_{1}$ model and we obtained a value of 2.02. According to the interpretation suggested in Kass and Raftery (1995, p. 777) evidence, given by the data, to reject the ${\cal M}_{0}$ model compared to the ${\cal M}_{1}$ model is thus decisive. Indeed, log₁₀B ₁₀ > 2 simply implies that given the data model ${\cal M}_{1}$ is more than 100 times more likely than the model ${\cal M}_{0}$.

Synthetic data

We then applied our Bayesian framework to calibrate the parameters of the 100 synthetic velocity profiles generated using the MAP estimators of model ${\cal M}_1$. We used the same marginal prior distributions as before (see Table 3). For each synthetic avalanche, we generated 20 000 MH samples and the last 10 000 were kept. Convergence was verified with the same diagnoses as for the measured data. We could then determine the ability of the ${\cal M}_0$ and ${\cal M}_1$ models to retrieve the true model parameter values used for the synthetic data generation. For this, we calculated the 90% coverage rates for the different parameters, this means, the number of 90% credibility intervals recovering the true value of Table 2. According to the results of Table 5, ${\cal M}_1$ model has a much better ability to determine the true model parameter values compared to the ${\cal M}_0$ model. This especially applies to σ ² for which ${\cal M}_0$ model is fully unable to identify the true value and, to a lower extent, to μ, for which the true value is in the posterior 90% credibility interval less than two times over three. Conversely, 90% coverage rates with ${\cal M}_1$ model are rather fair, varying between 82 and 98% for all the four parameters.

Table 5. The 90% coverage rates for the synthetic data sample. For each parameter and both models, the 90% coverage rate corresponds to the number of times over the sample of 100 synthetic avalanches for which the 90% posterior credibility interval includes the true value used for data generation. In other words, 90% is the perfect validation score. Parameter ϕ applies to model ${\cal M}_{1}$ only

To further compare the precision of the estimates led by the two models, the MAP estimators corresponding to the full synthetic sample are presented in Figures 8a–d. In mean, both models result in an underestimation of the μ parameter, but the underestimation with ${\cal M}_{1}$ is much smaller than with ${\cal M}_{0}$. In addition, the true value is well within the range of variability of the different estimates corresponding to the 100 synthetic avalanches with ${\cal M}_{1}$ whereas it is clearly outside with ${\cal M}_{0}$. Similarly, parameter σ ² is much better estimated by ${\cal M}_{1}$ model. Parameter ξ is slightly underestimated by ${\cal M}_{0}$ model and slightly overestimated by ${\cal M}_{1}$ model, but both models perform reasonably. Finally, ${\cal M}_{1}$ model estimates correctly the autocorrelation parameter ϕ. Overall, this analysis confirms that only ${\cal M}_{1}$ model is able to retrieve the true parameter values as soon as autocorrelation between measurements actually exists. In other words, not accounting for autocorrelation within model calibration when autocorrelation actually exists carries high risk to lead to biased estimates. This provides another very strong argument in favour of ${\cal M}_{1}$ since, as evidenced by the case study, significant autocorrelation may indeed exist between measurements made along the same avalanche flow.

Fig. 8. (a–d) Boxplots of the MAP estimators under models ${\cal M}_{0}$ (red colour) and ${\cal M}_{1}$ (blue colour) corresponding to the 100 synthetic avalanches. The true parameter values used for synthetic data generation are shown with a green colour line. (e) Boxplot of log₁₀B ₁₀ obtained for the 100 synthetic avalanches. The value of 2 in green corresponds to the reference value of Kass and Raftery (1995) above which evidence in favour of ${\cal M}_{1}$ is decisive.

Eventually, the common logarithm of the Bayes factor calculated for each synthetic avalanche (see Fig. 8e) indicates for all synthetic velocity series decisive evidence to reject the ${\cal M}_{0}$ model compared to the ${\cal M}_{1}$ model. This result is all the more logical given that the Bayes factor asymptotically selects the true model when it is included within a sample of competing models. However, this can be seen as a last strong point to advise using model ${\cal M}_{1}$ as soon as autocorrelation is suspected.

Conclusion and outlook

In this work, a Bayesian calibration of an avalanche flow dynamics model from data of high temporal resolution [1s] was developed. The objective of this work was to show how potential autocorrelation between measurements made along an avalanche flow can be considered within the calibration of an avalanche model, and to demonstrate that this improves the accuracy of friction parameters estimation and velocity predictions.

Two statistical models representing the discrepancies between observations and simulations were proposed: the first one, ${\cal M}_0$, classically considered the errors as independent and identically normally distributed, and the second one, ${\cal M}_1$, modelled the errors as an autoregressive process of order 1. The latter accounts for potential autocorrelation between measures made along the same avalanche flow, a question which has been poorly addressed in the snow science literature so far. Our objective was to determine the accuracy of parameter estimation when the autocorrelation is included or not in the modelling. More generally, we wanted to develop a framework able to link in a more mathematically consistent way, the rich and diverse data now increasingly available from avalanche test-sites with numerical propagation models.

Application was made on a well-documented avalanche event and on synthetic data. The avalanche was released at the Lautaret full-scale test-site on 13 February 2013 and the corresponding velocity time series was obtained from high rate positioning photogrammetry. A synthetic data set of 100 velocities time series mimicking the avalanche was created in order to further test our approach and to illustrate its benefits. Results for the avalanche showed strong and significant autocorrelation between predicted velocities. It also appeared that the velocity time series was correctly modelled by the statistical model ${\cal M}_1$ only (Ljung–Box test p-values > 0.05) so that there was decisive evidence to reject ${\cal M}_0$ compared to ${\cal M}_1$. In addition, resulting posterior estimates for friction parameters μ and ξ and velocities along the path were shown to be different and arguably more realistic (for the estimated velocities) with ${\cal M}_1$ than those obtained with ${\cal M}_0$. Eventually, our synthetic data confirmed that in the presence of autocorrelation between measurements, not accounting for it may lead to biased estimates. This all demonstrates the necessity and usefulness to explicitly account for autocorrelation within the calibration, both to realistically predict avalanche characteristics and to get friction parameter estimates that can be related to snow and topographical characteristics in a meaningful way, two points of crucial importance in snow science. Hence, since avalanche velocities are time series, as a good practice of modelling, potential autocorrelation of the errors should always be envisaged, since we clearly demonstrate that neglecting this effect may lead to undesirable consequences.

To develop and illustrate our approach, we chose the simple sliding block model with the Voellmy friction law because it is faster to run than more computationally intensive state of the art avalanche models. As it was mentioned in the Introduction, we used this simple model as a proof of concept but our approach could be generalized to other friction laws (e.g. the Coulomb one which has the advantage of having a single parameter (μ)), and to other ways of describing avalanche flows (e.g. depth-averaged equations or a sliding block with a PCM formulation, which would avoid the estimation of the flow depth) after suitable adaptation of the algorithm. Also, on the existing basis, it should be possible to include within the calibration other quantities which can be measured simultaneously as, for example, the snow depth in the release area, the runout distance, the deposited volume, etc. This would better constrain the calibration, leading to potentially less uncertain posterior estimates. In such a case, an appropriate likelihood would need to be proposed and a higher number of observations would be required to calibrate the new likelihood parameters. To overcome this shortcoming using several avalanches altogether could be an option which was out of the scope of this work.

One of the main drawbacks of a Bayesian approach is the high number of simulations that are needed (at least 1000 simulations) to reach convergence in a Markov Chain Monte Carlo setting. Also, a compromise between the number of parameters that can be estimated and the size of the data must be found. Hence, whether or not our approach can be practically implemented with the most computationally intensive avalanche model and how many unknown parameters can be inferred as a function of the available data sets remain the questions to be investigated.

Finally, it is worth to mention that the avalanche size and velocities we found are typical of a medium size avalanche (class 3 on the CAA international avalanche scale). Specifically, friction parameter values correspond to the range of values that can be found in the literature. In future work, it could therefore be informative to apply our approach to a much large sample of avalanches to exploit its potential for inferring relevant physics (e.g. to find the relationship between friction parameters and snow conditions). Ultimately, this all should (1) lead to a more accurate evaluation of the highest-return-period events required for avalanche risk assessment, and (2) simultaneously improve our knowledge of the relevant physics by providing sharper quantification of underlying processes.