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

María Belén Heredia, Nicolas Eckert, Clémentine Prieur, Emmanuel Thibert

► To cite this version:

María Belén Heredia, Nicolas Eckert, Clémentine Prieur, Emmanuel Thibert. Bayesian calibration of an avalanche model from autocorrelated measurements along the flow: application to velocities extracted from photogrammetric images. Journal of Glaciology, 2020, 66 (257), pp.373-385. 10.1017/jog.2020.11. hal-03201750

HAL Id: hal-03201750

https://hal.science/hal-03201750

Submitted on 19 Apr 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Paper

**Cite this article:** Heredia MB, Eckert N, Prieur C, Thibert E (2020). Bayesian calibration of an avalanche model from autocorrelated measurements along the flow: application to velocities extracted from photogrammetric images. *Journal of Glaciology* 66(257), 373–385.
https://doi.org/10.1017/jog.2020.11

Received: 8 August 2019
Revised: 21 January 2020
Accepted: 21 January 2020
First published online: 18 March 2020

**Keywords:**
Avalanches; glaciological instruments and methods; glaciological natural hazards

**Author for correspondence:**
María Belén Heredia,
E-mail: maria-belen.heredia@inrae.fr

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

María Belén Heredia¹, Nicolas Eckert¹, Clémentine Prieur²
and Emmanuel Thibert¹

¹Univ. Grenoble Alpes, INRAE, UR ETGR, Grenoble, France and ²Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP (Institute of Engineering Univ. Grenoble Alpes), LJK, 38000 Grenoble, France

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 propa- gation, potential autocorrelation should be explicitly accounted for. To this aim, this work pro- poses 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 sam-
ples (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 observa-
tions. 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

© The Author(s) 2020. This is an Open Access
article, distributed under the terms of the
Creative Commons Attribution licence (http://
creativecommons.org/licenses/by/4.0/), which
permits unrestricted re-use, distribution, and
reproduction in any medium, provided the
original work is properly cited.
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; Schaeffli 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:

$$ \frac{du}{dt} = g \sin \theta - \frac{F}{m}, \qquad (1) $$

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

$$ F = \mu mg \cos \theta + \frac{mg}{\xi h} u^2, \qquad (2) $$

depending on two friction parameters $\Theta = \{\mu, \xi\}$. Often it is assumed that $\mu$ evolves with the physical properties of snow whereas $\xi$ 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 $(\mu, \xi)$ 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_{\text{start}}\}$ where $T$ defines the topography of the terrain, $h$ the mean flow depth of the avalanche and $x_{\text{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 \cdot s^{-1}$) and position along the slope ($m$) at a time $t$ ($s$). The model depends on parameters $\Theta \in \mathbb{R}^p$ and forcing variables $x \in \mathbb{R}^d$. The observed velocity, collected on the field, at time $t$ is noted by $v_t$ and we denote $v_{\text{obs}} = \{v_1, \dots, 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 $\Theta$ that minimizes the discrepancy between the observations $v_{\text{obs}}$ and the model simulation $f(\Theta, x)$ and (2) to rigorously quantify the associated uncertainty. To this end, we use the generic statistical model:

$$ \mathcal{M}: v_t = f_t(x, \Theta) + \epsilon_t, \quad \forall t \in \{1, \dots, n\}, \qquad (3) $$

where $f_t(\Theta, x)$ denotes the simulation of the avalanche velocity at time $t$ and $\epsilon_t$ is the model error. With this classical additive formulation, propagation model errors and observation errors are modelled altogether in the residuals $\epsilon_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$ is the $\mathcal{N}(0, \sigma^2)$, $\forall t \in \{1, \dots, n\}$, where $\mathcal{N}$ denotes the Normal distribution and $\sigma^2$ 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 $\epsilon_t$ are likely to present a non-negligible correlation between two consecutive time steps (in other words, $\epsilon_t$ and $\epsilon_{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 $\epsilon_t$ can be expressed as:

$$ \epsilon_t = \phi \epsilon_{t-1} + \eta_t, \quad \eta_t \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \sigma^2) \quad \forall t \in \{1, \dots, n\}, \qquad (4) $$

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

Hereafter, the model with the assumption of normally distributed and independent errors is denoted $\mathcal{M}_0$ and the model with AR1 errors is denoted $\mathcal{M}_1$. Depending on the $\mathcal{M}_i$ model with
<table><caption>**Table 1.** Considered statistical models and corresponding parameters</caption><thead><tr><th>Model</th><th>Errors</th><th>Parameters</th></tr></thead><tbody><tr><td>M<sub>0</sub></td><td>Independent and identically distributed N(0, σ<sup>2</sup>)</td><td>Θ, σ<sup>2</sup></td></tr><tr><td>M<sub>1</sub></td><td>AR1</td><td>Θ, φ, σ<sup>2</sup></td></tr></tbody></table>

$i \in \{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 $M_1$ with $\phi = 0$ corresponds to $M_0$.

## 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 $\gamma$, the set of parameters related to the errors, it means $\gamma = \{\sigma^2\}$ for $M_0$ and $\gamma = \{\phi, \sigma^2\}$ for $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:

$$ \pi(\Theta, \gamma | v_{\text{obs}}, x) \propto \mathcal{L}(v_{\text{obs}} | \Theta, \gamma, x) \pi(\Theta, \gamma), \quad (5) $$

where $\mathcal{L}(v_{\text{obs}}|\Theta, \gamma, x)$ is the likelihood of the observations and $\pi(\Theta, \gamma)$ is the joint prior distribution. Thus, the likelihood $\mathcal{L}(v_{\text{obs}}|\Theta, \gamma, x)$ is required. Under the $M_0$ model, the observations follow a Gaussian distribution and the likelihood $\mathcal{L}(v_{\text{obs}}|\Theta, \gamma, x, M_0)$ writes as:

$$ \frac{1}{(2\pi\sigma^2)^{n/2}} \exp \left[ -\frac{1}{2\sigma^2} \sum_{t=1}^{n} (v_t - f_t(\Theta, x))^2 \right]. \quad (6) $$

Under the $M_1$ model, the likelihood $\mathcal{L}(v_{\text{obs}}|\Theta, \gamma, x, M_1)$ writes as (for more detail, see Appendix A):

$$ \sqrt{\frac{1 - \phi^2}{2\pi\sigma^2}} \exp\left[-\frac{1 - \phi^2}{2\sigma^2}(v_1 - f_1(\Theta, x))^2\right] \times \frac{1}{(2\pi\sigma^2)^{(n-1)/2}} \exp\left[-\frac{1}{2\sigma^2}\sum_{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 $\mu$ (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 $\pi(\Theta, \gamma|v_{\text{obs}}, x)$. The Markov chain returned by the algorithm can be considered as a sample from $\pi(\Theta, \gamma|v_{\text{obs}}, x)$.

The following description of the MH algorithm is obtained from Robert (2015). Let us denote, for simplicity, $\psi = \{\Theta, \gamma\}$ the set of the error and model parameters. For the application of the MH algorithm, it is needed an initial value $\psi^{(0)}$ and a proposal distribution $q$. Each iteration $k$ of the algorithm consists in:

1. Generating $\psi' \sim q(\cdot|\psi^{(k-1)})$.

2. Calculating $u \sim U(0, 1)$.

3. Taking

$$ \psi^{(k)} = \begin{cases} \psi' & \text{if } u \le \alpha, \\ \psi^{(k-1)} & \text{otherwise} \end{cases} \quad (8) $$

with

$$ \alpha = \min \left( \frac{\mathcal{L}(v_{\text{obs}} | \psi', x) \pi(\psi')}{\mathcal{L}(v_{\text{obs}} | \psi^{(k-1)}, x) \pi(\psi^{(k-1)})} \frac{q(\psi^{(k-1)} | \psi')}{q(\psi' | \psi^{(k-1)})}, 1 \right), \quad (9) $$

where we recall that $\mathcal{L}(v_{\text{obs}}|\psi, x)$ and $\pi(\psi) = \pi(\Theta, \gamma)$ 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 $\psi^{(k-1)}$ and with a given covariance $\Sigma_q$, which defines a random walk. For our application, $\Sigma_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 $M_0$ and $M_1$, we use the Bayes factor. The Bayes factor is a criterion of the evidence provided by the data (in our case $v_{\text{obs}}$) to reject a model ($M_0$) compared to another one ($M_1$). The Bayes factor is defined as the ratio of marginal probabilities (see more details in Kass and Raftery, 1995):

$$ B_{10} = \frac{p(v_{\text{obs}} | M_1, x)}{p(v_{\text{obs}} | M_0, x)}, \quad (10) $$

where

$$ p(v_{\text{obs}} | M_i, x) = \int \mathcal{L}(v_{\text{obs}} | \psi_i, M_i, x) \pi(\psi_i | M_i, x) \, d\psi_i, \quad (11) \\ i \in \{0, 1\}. $$

The conditioning by $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(\psi_i|v_{\text{obs}}, x, M_i)$ (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 $M_0$ compared to $M_1$ according to the value of $\log_{10} B_{10}$ is proposed in the work of Kass and Raftery (1995), where $\log_{10}$ denotes the common logarithm. Notably, a $\log_{10} B_{10}$ value between 1 and 2 indicates strong
<figure><img src="image_2.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>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).</figcaption></figure>

evidence in favour of rejecting $M_0$ compared to $M_1$, and $\log_{10} B_{10} > 2$ a decisive evidence in favour of $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.

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_{\text{obs}}$. At each time $t$, we computed the centre of mass position of the avalanche using a set of front points $D_t \in \mathbb{R}^3$ from the delineated front. $D_t$ excludes the lateral front points because they are less active in the avalanche flow (Pulfer and others, 2013):

$$c_t = \frac{\sum_{p_t \in D_t} p_t}{n_t}, \qquad (12)$$

where $n_t$ is the number of elements in $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:

$$v_t = \frac{c_t - c_{t-1}}{\Delta t}. \qquad (13)$$

In our application, $\Delta t = 1$ s 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
<table><caption>**Table 2.** Parameter values used for the generation of the synthetic data set. They correspond to the maximum (MAP) obtained with model $M_1$ for the avalanche of the 13 February 2013</caption><thead><tr><th>Parameter</th><th>MAP estimate</th></tr></thead><tbody><tr><td>μ<sub>MAP</sub></td><td>0.32</td></tr><tr><td>ξ<sub>MAP</sub>[m·s<sup>-2</sup>]</td><td>495.73</td></tr><tr><td>σ<sup>2</sup><sub>MAP</sub>[m<sup>2</sup>·s<sup>-2</sup>]</td><td>8.07</td></tr><tr><td>φ<sub>MAP</sub></td><td>0.65</td></tr></tbody></table>

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 $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 $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 $μ$, $ξ$, $σ^2$ and $φ$ parameters of $M_1$ model were determined for the avalanche, we proceed as follows:

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. An AR(1) error of parameters $σ^2_{MAP}$ and $φ_{MAP}$ was simulated.

3. The AR(1) error simulated was added to $f(Θ_{MAP}, x)$. If negative velocity was obtained, the AR(1) error was resampled.

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.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Samples of the generated synthetic longitudinal velocity profiles (colour lines) and observations corresponding to the studied avalanche (dots).</figcaption></figure>

<table><caption>**Table 3.** Parameters marginal prior distributions. The prior distribution of $φ$ applies only to $M_1$ model. $Γ$ represents the Gamma distribution and $U$ the uniform distribution</caption><tr><td>μ ~ Γ(9, 0.03)</td><td>σ² ~ InvGamma(3, 20)</td></tr><tr><td>ξ ~ Γ(25, 40)</td><td>φ ~ U(-1, 1)</td></tr></table>

## Application, results and discussion

### The avalanche

#### Prior distributions

Our Bayesian framework was used to calibrate the parameters of the $M_0$ and $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 $U(-1, 1)$. 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.

#### 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_0$: 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)
<figure><img src="image_5.png" /><figcaption>Fig. 4. Prior and posterior densities of model parameters. Panels (a-c) model $M_0$ and panels (d-g) model $M_1$.</figcaption></figure>

<table><caption>**Table 4.** Posterior distribution characteristics: posterior mean, posterior standard deviation (sd.), median (50% percentile) and 95% credibility interval (2.5% and 97.5% percentiles).</caption><thead><tr><th rowspan="2">Model</th><th colspan="3">M<sub>0</sub></th><th colspan="4">M<sub>1</sub></th></tr><tr><th>μ</th><th>ξ[m·s<sup>-2</sup>]</th><th>σ<sup>2</sup>[m<sup>2</sup>·s<sup>-2</sup>]</th><th>μ</th><th>ξ[m·s<sup>-2</sup>]</th><th>σ<sup>2</sup>[m<sup>2</sup>·s<sup>-2</sup>]</th><th>φ</th></tr></thead><tbody><tr><td>Mean</td><td>0.3</td><td>409.74</td><td>12.71</td><td>0.34</td><td>696.12</td><td>9.79</td><td>0.71</td></tr><tr><td>sd.</td><td>0.07</td><td>101.12</td><td>3.53</td><td>0.08</td><td>305.71</td><td>2.62</td><td>0.17</td></tr><tr><td>q<sub>0.025</sub></td><td>0.16</td><td>264.30</td><td>7.66</td><td>0.18</td><td>306.97</td><td>5.93</td><td>0.33</td></tr><tr><td>q<sub>0.50</sub></td><td>0.3</td><td>393.34</td><td>12.17</td><td>0.34</td><td>627.36</td><td>9.33</td><td>0.73</td></tr><tr><td>q<sub>0.975</sub></td><td>0.43</td><td>649.64</td><td>20.98</td><td>0.50</td><td>1475.67</td><td>16.15</td><td>0.96</td></tr></tbody></table>

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 $M_0$ densities are shown (respectively, $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 $M_0$ and $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.

The variances of the posterior distribution of the friction parameters $\mu$ and $\xi$ in the $M_0$ model are lower than in $M_1$ model. This result (arguably the sole undesirable one with $M_1$ compared to $M_0$) could be a consequence of the interactions between the friction coefficients and error parameters in the more parameterized $M_1$ model (four free parameters instead of three with $M_0$). Indeed, there is strong correlation between the parameters $\xi$ and $\phi$ (0.49) for the model $M_1$, which may preclude reaching sharp estimates of friction parameters with model $M_1$. Conversely, even if there is a high correlation between the

parameters $\mu$ and $\xi$ with both models, switching from $M_0$ to $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 $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 $\sigma^2$ and the friction parameters with $M_1$ than with $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 $M_1$), differences in posterior estimates reaches 13% for $\mu$ (relative difference between 0.34 with $M_1$ instead of 0.3 with $M_0$), and 70% for $\xi$ (relative difference between 696 m$^2$·s$^{-1}$ with model $M_1$ and 410 m$^2$·s$^{-1}$ with model $M_0$). This indicates that using either $M_0$ or $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 $\sigma^2$ from model $M_0$ to model $M_1$ (Table 4). Indeed, the errors variance $\sigma^2$ is higher if the autocorrelation $\phi$ is not included. Specifically, the standard deviation of the errors decreases from 3.6 to 3.1 m·s$^{-1}$ indicating that with model $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 $\phi$ is largely positive, with a posterior mean of 0.71. Also, its posterior 95% credibility interval whose lower bound is 0.33 firmly excludes the zero value corresponding to $M_0$. This pleads for a high and significant autocorrelation between velocities along flow propagation.

## Model posterior estimates

We calculated for model $M_0$ the posterior estimate of Voellmy model simulations $f(\Theta, x)$ from the posterior estimates of parameters $\mu$ and $\xi$ (Table 4). Resulting posterior estimates of model errors $\epsilon_t$ were analysed by applying the Ljung-Box test ($H_0$: 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, $\text{cor}(\epsilon_t, \epsilon_{t-1})$ for all $t$
<figure><img src="image_6.png" /><figcaption>Fig. 5. Joint posterior distribution of parameters of M₀ model highlighting inter-parameter correlations.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Joint posterior distribution of parameters of M₁ model highlighting inter-parameter correlations.</figcaption></figure>

> 1 is significant. Hence, one of the assumptions underlying statistical model M₀ 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 M₀ are thus fulfilled.

We proceeded similarly for statistical model M₁, applying the same statistical tests to the posterior estimates of the innovations ηᵢ of the AR(1) model. Remember that, in the M₁ model, the independent assumption is on the innovations ηᵢ and not on the errors εᵢ. 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 M₁ represents correctly the autocorrelation of the errors, and, more widely, that contrary to model M₁ 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(vₓ|vobs, x, Mᵢ) 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'(vₓ|vobs, x, Mᵢ). More precisely, the first writes:

$$p(v_x|v_{\text{obs}}, x, \mathcal{M}_i) = \int f(\Theta, x) p(\Theta|v_{\text{obs}}, x, \mathcal{M}_i) d\Theta, \quad (14)$$

and, the second one:

$$p'(v_x|v_{\text{obs}}, x, \mathcal{M}_i) = \int (f(\Theta, x) + \epsilon) p(\Theta, \gamma|v_{\text{obs}}, x, \mathcal{M}_i) d\Theta d\gamma. \quad (15)$$

Eventually, to quantify the added value of parameter inference, the prior distribution of the velocity was also calculated by
<figure><img src="image_8.png" /><figcaption>Fig. 7. Predictive velocity distributions versus data: (a) model $M_0$ and (b) model $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.</figcaption></figure>

integrating over the prior distribution on the friction and model error parameters:

$$ \pi(v_x|x, M_i) = \int (f(\Theta, x) + \epsilon) \pi(\Theta, \gamma|M_i) d\Theta d\gamma. \quad (16) $$

By using a Bayesian approach, we obtained a sample of parameters $(\mu, \xi, \sigma^2)$ for $M_0$ (and also of $\phi$ for $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 $M_0$, but also exists with $M_1$. This simply reflects the larger a posterior variance of model $M_1$ parameters highlighted above.

In detail, almost all the elements of $v_{obs}$, except two with model $M_0$ and one with model $M_1$, are inside the 90% posterior credibility intervals drawn on $p'(v_x|v_{obs}, x, M_i)$ (blue dotted lines Fig. 7), so that from this perspective the added value of $M_1$ is not obvious. However, why $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(v_x|v_{obs}, x, M_0)$, compared to 62% with $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 $M_0$ (Fig. 7a, black line) and $M_1$ (Fig. 7b, black line). If we denote $q50_{v_x,M_i}$ the median (50% percentile) of the posterior predictive distribution of velocities under model $M_i$, the first writes:

$$ \Delta_1 = \frac{1}{240} \int_0^{240} (q50_{v_x,M_1} - q50_{v_x,M_0}) dx $$

and the second:

$$ \Delta_2 = \sqrt{\frac{1}{240} \int_0^{240} (q50_{v_x,M_1} - q50_{v_x,M_0})^2 dx} $$

Obtained values are $\Delta_1 = 1.31 \text{ m} \cdot \text{s}^{-1}$ and $\Delta_2 = 1.72 \text{ m} \cdot \text{s}^{-1}$. 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 \text{ m}$ and $x=250 \text{ m}$ (plateau phase), around $15 \text{ m} \cdot \text{s}^{-1}$ with $M_0$. The underestimation with $M_0$ with regards to $M_1$ can be estimated by the ratio between $\Delta_2$ and this velocity to 11.5%.

All in all, the analysis of predictive velocity distributions further confirms that model $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_{10} B_{10}$ was calculated to compare the $M_0$ model to the $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 $M_0$ model compared to the $M_1$ model is thus decisive. Indeed, $log_{10} B_{10} > 2$ simply implies that given the data model $M_1$ is more than 100 times more likely than the model $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 $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 $M_0$ and $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,
<table><caption>**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 $\phi$ applies to model $M_1$ only</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="3">M<sub>0</sub></th><th colspan="4">M<sub>1</sub></th></tr><tr><th>μ</th><th>ξ</th><th>σ<sup>2</sup></th><th>μ</th><th>ξ</th><th>σ<sup>2</sup></th><th>φ</th></tr></thead><tbody><tr><td>Coverage rate</td><td>64</td><td>84</td><td>0</td><td>98</td><td>82</td><td>98</td><td>93</td></tr></tbody></table>

the number of 90% credibility intervals recovering the true value of Table 2. According to the results of Table 5, $M_1$ model has a much better ability to determine the true model parameter values compared to the $M_0$ model. This especially applies to $\sigma^2$ for which $M_0$ model is fully unable to identify the true value and, to a lower extent, to $\mu$, for which the true value is in the posterior 90% credibility interval less than two times over three. Conversely, 90% coverage rates with $M_1$ model are rather fair, varying between 82 and 98% for all the four parameters.

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 $\mu$ parameter, but the underestimation with $M_1$ is much smaller than with $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 $M_1$ whereas it is clearly outside with $M_0$. Similarly, parameter $\sigma^2$ is much better estimated by $M_1$ model. Parameter $\xi$ is slightly underestimated by $M_0$ model and slightly overestimated by $M_1$ model, but both models perform reasonably. Finally, $M_1$ model estimates correctly the autocorrelation parameter $\phi$. Overall, this analysis confirms that only $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 $M_1$ since, as evidenced by the case study, significant autocorrelation may indeed exist between measurements made along the same avalanche flow.

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 $M_0$ model compared to the $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 $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, $M_0$, classically considered the errors as independent and identically normally distributed, and the second one, $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 $M_1$ only (Ljung-Box test p-values > 0.05) so that there was decisive evidence to reject $M_0$ compared to $M_1$. In addition, resulting posterior estimates for friction parameters $\mu$ and $\xi$ and velocities along the path were shown to be different and arguably more realistic (for the estimated velocities) with $M_1$ than those obtained with $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

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** (a-d) Boxplots of the MAP estimators under models $M_0$ (red colour) and $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_{10} B_{10}$ 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 $M_1$ is decisive.</figcaption></figure>
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.

**Acknowledgments.** M.B. Heredia holds a PhD grant from Labex OSUG@2020. Within the CDP-Trajectories context, this work is supported by the French National Research Agency in the framework of the *Investissements d'avenir program* (ANR-15-IDEX-02). We thank Eric Parent for his valuable advice regarding the calculation of the Bayes factor. We are grateful to the two anonymous reviewers and the scientific editors Hester Jiskoot and Alec Van Herwijnen for their comments which helped us to improve our manuscript.

## References

*Ancey C* (2005) Monte Carlo calibration of avalanches described as Coulomb fluid flows. *Philosophical Transactions: Mathematical, Physical and Engineering Sciences* **363**(1832), 1529–1550.

*Ancey C* (2006) *Dynamique des avalanches*. Lausanne, Switzerland: Presses Polytechniques et Universitaires Romandes.

*Ancey C* and *Meunier M* (2004) Estimating bulk rheological properties of flowing snow avalanches from field data. *Journal of Geophysical Research: Earth Surface* **109**(F01004). doi: 10.1029/2003JF000036

*Ancey C*, *Meunier M* and *Richard D* (2003) Inverse problem in avalanche dynamics models. *Water Resources Research* **39**(4), doi: 10.1029/2002WR001749

*Barbolini M* and *Savi F* (2001) Estimate of uncertainties in avalanche hazard mapping. *Annals of Glaciology* **32**, 299–305. doi: 10.3189/172756401781819373

*Bartelt P*, *Bühler Y*, *Buser O*, *Christen M* and *Meier L* (2012) Modeling mass-dependent flow regime transitions to predict the stopping and depositional behavior of snow avalanches. *Journal of Geophysical Research: Earth Surface* **117**(F1). doi: 10.1029/2010JF001957

*Bartelt P*, *Buser O*, *Vera Valero C* and *Bühler Y* (2016) Configurational energy and the formation of mixed flowing/powder snow and ice avalanches. *Annals of Glaciology* **57**(71), 179–188. doi: 10.3189/2016AoG71A464

*Bayes T* (1763) LII. An essay towards solving a problem in the doctrine of chances. By the late Rev. Mr. Bayes, FRS communicated by Mr. Price, in a letter to John Canton, AMFR S. *Philosophical Transactions of the Royal Society of London* **53**, 370–418.

*Berliner LM* (2003) Physical-statistical modeling in geophysics. *Journal of Geophysical Research: Atmospheres* **108**(D24). doi: 10.1029/2002JD002865

*Bernardo JM* and *Smith AF* (2009) *Bayesian theory, volume 405*. John Wiley & Sons, Baffins Lane, Chichester, England. doi: 10.1002/9780470316870

*Borstad CP* and *McClung D* (2009) Sensitivity analyses in snow avalanche dynamics modeling and implications when modeling extreme events. *Canadian Geotechnical Journal* **46**(9), 1024–1033. doi: 10.1139/T09-042

*Buhler R*, *Argue C*, *Jamieson B* and *Jones A* (2018) Sensitivity Analysis of the RAMMS Avalanche Dynamics Model in a Canadian Transitional Snow Climate, iSSW, Innsbruck, Austria.

*Bühler Y*, *Christen M*, *Kowalski J* and *Bartelt P* (2011) Sensitivity of snow avalanche simulations to digital elevation model quality and resolution. *Annals of Glaciology* **52**(58), 72–80. doi: 10.3189/172756411797252121

*Carmassi M*, *Barbillon P*, *Keller M*, *Parent E* and *Chiodetti M* (2018) Bayesian calibration of a numerical code for prediction. *ArXiv e-prints*.

*Casassa G*, *Narita H* and *Maeno N* (1989) Measurements of friction coefficients of snow blocks. *Annals of Glaciology* **13**, 40–44. doi: 10.3189/S0260305500007618

*Castebrunet H*, *Eckert N* and *Giraud G* (2012) Snow and weather climatic control on snow avalanche occurrence fluctuations over 50 yr in the french alps. *Climate of the Past* **8**(2), 855–875. doi: 10.5194/cp-8-855-2012

*Christen M*, *Kowalski J* and *Bartelt P* (2010) Ramms: numerical simulation of dense snow avalanches in three-dimensional terrain. *Cold Regions Science and Technology* **63**(1), 1–14. doi: 10.1016/j.coldregions.2010.04.005

*Clark JS* (2005) Why environmental scientists are becoming Bayesians. *Ecology Letters* **8**(1), 2–14. doi: 10.1111/j.1461-0248.2004.00702.x

*Dent J* and *5 others* (1998) Density, velocity and friction measurements in a dry-snow avalanche. *Annals of Glaciology* **26**, 247–252. doi: 10.3189/1998AoG26-1-247-252

*Dent JD* and *Lang TE* (1980) Modeling of snow flow. *Journal of Glaciology* **26**(94), 131–140. doi: 10.3189/S0022143000010674

*Dreier L*, *Bühler Y*, *Ginzler C* and *Bartelt P* (2016) Comparison of simulated powder snow avalanches with photogrammetric measurements. *Annals of Glaciology* **57**(71), 371–381. doi: 10.3189/2016AoG71A532

*Eckert N* and *6 others* (2012) Quantitative risk and optimal design approaches in the snow avalanche field: review and extensions. *Cold Regions Science and Technology* **79–80**, 1–19. doi: 10.1016/j.coldregions.2012.03.003

*Eckert N*, *Naaim M* and *Parent E* (2010) Long-term avalanche hazard assessment with a Bayesian depth averaged propagation model. *Journal of Glaciology* **56**, 563–586. doi: 10.3189/002214310793146331

*Eckert N*, *Parent E*, *Faug T* and *Naaim M* (2009) Bayesian optimal design of an avalanche dam using a multivariate numerical avalanche model. *Stochastic Environmental Research and Risk Assessment* **23**(8), 1123–1141. doi: 10.1007/s00477-008-0287-6

*Eckert N*, *Parent E*, *Naaim M* and *Richard D* (2008) Bayesian stochastic modelling for avalanche predetermination: from a general system framework to return period computations. *Stochastic Environmental Research and Risk Assessment* **22**(2), 185–206. doi: 10.1007/s00477-007-0107-4
Eckert N, Parent E and Richard D (2007) Revisiting statistical-topographical methods for avalanche predetermination: Bayesian modelling for runout distance predictive distribution. *Cold Regions Science and Technology* **49**(1), 88–107, selected Papers from the General Assembly of the European Geosciences Union (EGU), Vienna, Austria, 25 April 2005. doi: **10.1016/j.col-dregions.2007.01.005**

Evin G, Thyer M, Kavetski D, McInerney D and Kuczera G (2014) Comparison of joint versus postprocessor approaches for hydrological uncertainty estimation accounting for error autocorrelation and heteroscedasticity. *Water Resources Research* **50**(3), 2350–2375. doi: **10.1002/2013WR014185**

Faig W, Shih T and Liu X (1990) A non-metric stereo system and its calibration. *Proceeding of the ACSM-ASPRS Annual Convention, American Society for Photogrammetry and Remote Sensing (ASPRS) and American Congress on Surveying and Mapping (ACSM), Denver, Colorado (USA) vol*, volume 5, pp. 18–25.

Faug T, Turnbull B and Gauer P (2018) Looking beyond the powder/dense flow avalanche dichotomy. *Journal of Geophysical Research: Earth Surface* **123**(6), 1183–1186. doi: **10.1002/2018JF004665**

Favier P, Bertrand D, Eckert N and Naaim M (2014a) A reliability assessment of physical vulnerability of reinforced concrete walls loaded by snow avalanches. *Natural Hazards and Earth System Sciences* **14**(3), 689–704. doi: **10.5194/nhess-14-689-2014**

Favier P, Eckert N, Bertrand D and Naaim M (2014b) Sensitivity of avalanche risk to vulnerability relations. *Cold Regions Science and Technology* **108**, 163–177. doi: **10.1016/j.coldregions.2014.08.009**

Fierz C and 8 others (2009) *The International Classification for Seasonal Snow on the Ground*, IHP-VII Technical Documents in Hydrology N83, IACS Contribution No1. UNESCO/IHP, Paris.

Fischer JT (2013) A novel approach to evaluate and compare computational snow avalanche simulation. *Natural Hazards and Earth System Sciences* **13**(6), 1655–1667. doi: **10.5194/nhess-13-1655-2013**

Fischer JT, Fromm R, Gauer P and Sovilla B (2014) Evaluation of probabilistic snow avalanche simulation ensembles with Doppler radar observations. *Cold Regions Science and Technology* **97**(Supplement C), 151–158. doi: **10.1016/j.coldregions.2013.09.011**

Fischer JT, Kofler A, Fellin W, Granig M and Kleemayr K (2015) Multivariate parameter optimization for computational snow avalanche simulation. *Journal of Glaciology* **61**(229), 875–888. doi: **10.3189/2015JoG14J168**

Fisher RA (1922) On the mathematical foundations of theoretical statistics. *Philosophical Transactions of the Royal Society of London. Series A*, Containing Papers of a Mathematical or Physical Character **222**(594–604), 309–368. doi: **10.1098/rsta.1922.0009**

Gauer P (2014) Comparison of avalanche front velocity measurements and implications for avalanche models. *Cold Regions Science and Technology* **97**, 132–150. doi: **10.1016/j.coldregions.2013.09.010**

Gauer P, Medina-Cetina Z, Lied K and Kristensen K (2009) Optimization and probabilistic calibration of avalanche block models. *Cold Regions Science and Technology* **59**(2), 251–258, international Snow Science Workshop (ISSW) 2008. doi: **10.1016/j.coldregions.2009.02.002**

Gaume J, Gast T, Teran J, van Herwijnen A and Jiang C (2018) Dynamic anticrack propagation in snow. *Nature Communications* **9**, 3047. doi: **10.1038/s41467-018-05181-w**

Gelman A and Rubin DB (1992) Inference from iterative simulation using multiple sequences. *Statistical Science* **7**(4), 457–472.

Gilks WR, Richardson S and Spiegelhalter D (1995) *Markov chain Monte Carlo in practice*. Chapman and Hall/CRC, London. doi: **10.1201/b14835**

Gouriérox C and Monfort A (1995) *Séries temporelles et modèles dynamiques*. Economica, Paris.

Kass RE and Raftery AE (1993) *Bayes Factors and Model Uncertainty* (Technical report No. 254). University of Washington, Department of Statistics, Seattle, Washington.

Kass RE and Raftery AE (1995) Bayes Factors. *Journal of the American Statistical Association* **90**(430), 773–795.

Kern M, Bartelt P, Sovilla B and Buser O (2009) Measured shear rates in large dry and wet snow avalanches. *Journal of Glaciology* **55**(190), 327–338. doi: **10.3189/00221430978860874**

Keylock CJ, McClung DM and Magnússon MM (1999) Avalanche risk mapping by simulation. *Journal of Glaciology* **45**(150), 303–314. doi: **10.3189/S0022143000001805**

Köhler A, McElwaine JN, Sovilla B, Ash M and Brennan P (2016) The dynamics of surges in the 3 February 2015 avalanches in Vallée de la

Sionne. *Journal of Geophysical Research: Earth Surface* **121**(11), 2192–2210. doi: **10.1002/2016JF003887**

Kuczera G (1983) Improved parameter inference in catchment models: 1. Evaluating parameter uncertainty. *Water Resources Research* **19**(5), 1151–1162. doi: **10.1029/WR019i005p01151**

Kuczera G and Parent E (1998) Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm. *Journal of Hydrology* **211**, 69–85. doi: **10.1016/S0022-1694(98)00198-X**

McInerney D, Thyer M, Kavetski D, Lerat J and Kuczera G (2017) Improving probabilistic prediction of daily streamflow by identifying Pareto optimal approaches for modeling heteroscedastic residual errors. *Water Resources Research* **53**(3), 2199–2239. doi: **10.1002/2016WR019168**

Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH and Teller E (1953) Equation of state calculations by fast computing machines. *The Journal of Chemical Physics* **21**(6), 1087–1092. doi: **10.1063/1.1699114**

Meunier M, Ancey C and Taillandier JM (2004) Fitting avalanche-dynamics models with documented events from the Col du Lautaret site (France) using the conceptual approach. *Cold Regions Science and Technology* **39**(1), 55–66. doi: **10.1016/j.coldregions.2004.03.004**

Naaim M, Durand Y, Eckert N and Chambon G (2013) Dense avalanche friction coefficients: influence of physical properties of snow. *Journal of Glaciology* **59**(216), 771–782. doi: **10.3189/2013JoG12J205**

Naaim M, Naaim-Bouvet F, Faug T and Bouchet A (2004) Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects. *Cold Regions Science and Technology* **39**(2), 193–204, snow And Avalanches: Papers Presented At The European Geophysical Union Conference, Nice, April 2003. Dedicated To The Avalanche Dynamics Pioneer Dr. B. Salm. doi: **10.1016/j.coldregions.2004.07.001**

Oakley J and O'Hagan A (2002) Bayesian inference for the uncertainty distribution of computer model outputs. *Biometrika* **89**(4), 769–784. doi: **10.1093/biomet/89.4.769**

Plummer M, Best N, Cowles K and Vines K (2006) CODA: convergence diagnosis and output analysis for MCMC. *R News* **6**(1), 7–11.

Prokop A (2008) Assessing the applicability of terrestrial laser scanning for spatial snow depth measurements. *Cold Regions Science and Technology* **54**(3), 155–163, snow avalanche formation and dynamics, doi: **10.1016/j.coldregions.2008.07.002**.

Prokop A and 6 others (2015) Merging terrestrial laser scanning technology with photogrammetric and total station data for the determination of avalanche modeling parameters. *Cold Regions Science and Technology* **110**, 223–230. doi: **10.1016/j.coldregions.2014.11.009**

Pulfer G, Naaim M, Thibert E and Soruco A (2013) Retrieving avalanche basal friction law from high rate positioning of avalanches. *International Snow Science Workshop (ISSW)*, Irstea, ANENA, Meteo France, Grenoble – Chamonix Mont-Blanc, France, pp. 1418–1424.

R Core Team (2019) *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing.

Robert CP (2015) The Metropolis-Hastings algorithm. *ArXiv e-prints*.

Robert CP and Casella G (2009) *Introducing Monte Carlo Methods with R* (Use R), 1st ed. Berlin, Heidelberg: Springer-Verlag.

Rognon PG and 5 others (2008) Rheology of dense snow flows: inferences from steady state chute-flow experiments. *Journal of Rheology* **52**(3), 729–748. doi: **10.1122/1.2897609**

Salm B (2004) A short and personal history of snow avalanche dynamics. *Cold Regions Science and Technology* **39**(2), 83–92, Snow And Avalanches: Papers Presented At The European Geophysical Union Conference, Nice, April 2003. Dedicated To The Avalanche Dynamics Pioneer Dr. B. Salm. doi: **10.1016/j.coldregions.2004.06.004**

Salm B, Burkard A and Gubler H (1990) *Berendnung von fließlanwinen*, eine anleitung für paktiker mit beispielen (Technical Report 47). Davos: Eidgenössches Institut für Schnee und Lawinenforschung.

Schaeffi B and Kavetski D (2017) Bayesian spectral likelihood for hydrological parameter inference. *Water Resources Research* **53**(8), 6857–6884. doi: **10.1002/2016WR019465**

Schlappy R and 7 others (2014) Validation of extreme snow avalanches and related return periods derived from a statistical-dynamical model using tree-ring techniques. *Cold Regions Science and Technology* **99**, 12–26. doi: **10.1016/j.coldregions.2013.12.001**

Schweizer J, Bruce Jamieson J and Schneebeli M (2003) Snow avalanche formation. *Reviews of Geophysics* **41**(4). doi: **10.1029/2002RG000123**

Soruco A, Thibert E, Vincent C, Blanc R and Héno R (2011) Measurement of avalanche front velocity from high-speed terrestrial digital
photogrammetry. Geophysical Research Abstracts, EGU General Assembly, EGU2011-EGU8177, volume 13.

Sovilla B, Schaeer M, Kern M and Bartelt P (2008) Impact pressures and flow regimes in dense snow avalanches observed at the Vallée de la Sionne test site. Journal of Geophysical Research: Earth Surface 113(F1). doi: 10.1029/2006JF000688

Steinkogler W, Gaume J, Löwe H, Sovilla B and Lehning M (2015) Granulation of snow: from tumbler experiments to discrete element simulations. Journal of Geophysical Research: Earth Surface 120(6), 1107–1126. doi: 10.1002/2014JF003294

Steinkogler W, Sovilla B and Lehning M (2014) Influence of snow cover properties on avalanche dynamics. Cold Regions Science and Technology 97, 121–131. doi: 10.1016/j.coldregions.2013.10.002

Straub D and Grêt-Regamey A (2006) A Bayesian probabilistic framework for avalanche modelling based on observations. Cold Regions Science and Technology 46(3), 192–203. doi: 10.1016/j.coldregions.2006.08.024

Sun R, Yuan H and Liu X (2017) Effect of heteroscedasticity treatment in residual error models on model calibration and prediction uncertainty estimation. Journal of Hydrology 554, 680–692. doi: 10.1016/j.jhydrol.2017.09.041

Thibert E and 17 others (2015) The full-scale avalanche test-site at Lautaret pass (french alps). Cold Regions Science and Technology 115, 30–41. doi: 10.1016/j.coldregions.2015.03.005

Torre F, Boreux JJ and Parent E (2001) The Metropolis-Hastings algorithm, a handy tool for the practice of environmental model estimation: illustration with biochemical oxygen demand data. Cybergeo: European Journal of Geography [En ligne] (187), 9447–9463. doi: 0.4000/cybergeo.4750

Truong H, Keylock C, Eckert N, Bellot H and Naaïm M (2018) Refining the processing of paired time series data to improve velocity estimation in snow flows. Cold Regions Science and Technology 151, 75–88. doi: 10.1016/j.col-d regions.2018.03.004

Turnbull B and Bartelt P (2003) Mass and momentum balance model of a mixed flowing/powder snow avalanche. Surveys in Geophysics 24(5), 465–477.

Voelmy A (1964) On the destructive force of avalanches (translation by R. E. Tate). U.S. Department of Agriculture, Forest Service, Alta Avalanche Study Center National Forest No 2.

Vriend NM and 5 others (2013) High-resolution radar measurements of snow avalanches. Geophysical Research Letters 40(4), 727–731. doi: 10.1002/grl. 50134

## Appendix A: Likelihood of an AR1 process

Since $|\phi| < 1$, the $\epsilon_t$ error has a representation in terms of the innovations $\eta_t$ (also known as MA($\infty$)) representation, see Gouriéroux and Monfort, 1995 for more details):

$$ \epsilon_t = \phi^t \epsilon_0 + \sum_{i=0}^{t-1} \phi^i \eta_{t-i}, \qquad (A1) $$

where we recall that $\phi$ and $\eta_{t-i}$ are the coefficient and the innovations of the AR1 process, respectively (see eqn 4).

From the last equation, the unique stationary solution is obtained for $\epsilon_0$ with zero mean and variance equal to $\sigma^2/1 - \phi^2$, independent from ($\eta_1, \dots, \eta_t$), and we get that:

$$ E(\epsilon_t) = E(\epsilon_0) + \sum_{i=0}^{t-1} \phi^i E(\eta_{t-i}) = 0, \qquad (A2) $$

$$ \mathrm{Var}(\epsilon_t) = \phi^{2t} \mathrm{Var}(\epsilon_0) + \sum_{i=0}^{t-1} \phi^{2i} \mathrm{Var}(\eta_{t-i}) = \frac{\sigma^2}{(1 - \phi^2)}. \qquad (A3) $$

The marginal distribution of $\epsilon_t$ is Gaussian because it is the sum of independent and identically normally distributed $\eta_t$. In particular the first error term $\epsilon_1$ is Gaussian:

$$ \epsilon_1 \sim \mathcal{N}\left(0, \frac{\sigma^2}{(1 - \phi^2)}\right). \qquad (A4) $$

Considering that $\epsilon_2 = \phi\epsilon_1 + \eta_2$, the conditional distribution of $\epsilon_2$ given $\epsilon_1$ is $\epsilon_2|\epsilon_1 \sim \mathcal{N}(\phi\epsilon_1, \sigma^2)$. In a more general way:

$$ \epsilon_t|\epsilon_{t-1} \sim \mathcal{N}(\phi\epsilon_{t-1}, \sigma^2). \qquad (A5) $$

Then, the joint distribution of errors is:

$$ p(\epsilon_1, \dots, \epsilon_n | \phi, \sigma^2) = p(\epsilon_1) \prod_{i=2}^{n} p(\epsilon_i | \epsilon_{i-1}). \qquad (A6) $$

From this equation, the likelihood expression of eqn (7) is obtained.

## Appendix B: Numerical evaluation of the Bayes Factor

Thanks to the Bayes theorem, we can write:

$$ \frac{1}{p(v_{\mathrm{obs}}|M_i, x)} = \frac{\pi(\psi_i|v_{\mathrm{obs}}, M_i, x)}{\mathcal{L}(v_{\mathrm{obs}}|\psi_i, M_i, x)\pi(\psi_i|M_i, x)}, \qquad (A7) $$

where $i \in \{0, 1\}$.

If $g$ is a density function defined on an ensemble $\Omega$:

$$ \frac{1}{p(v_{\mathrm{obs}}|M_i, x)} = \frac{1}{p(v_{\mathrm{obs}}|M_i, x)} \int_{\psi_i \in \Omega} g(\psi_i) d\psi_i \\ = \int_{\psi_i \in \Omega} \frac{g(\psi_i) \pi(\psi_i|v_{\mathrm{obs}}, M_i, x)}{\mathcal{L}(v_{\mathrm{obs}}|\psi_i, M_i, x) \pi(\psi_i|M_i, x)} d\psi_i $$

From Monte Carlo simulations, $1/p(v_{\mathrm{obs}}|M_i, x)$ can thus be estimated as:

$$ \frac{1}{p(v_{\mathrm{obs}}|M_i, x)} = \frac{1}{N} \sum_{j=1}^{N} \frac{g(\psi_i^{(j)})}{\mathcal{L}(v_{\mathrm{obs}}|\psi_i^{(j)}, M_i, x) \pi(\psi_i^{(j)}|M_i, x)}, \qquad (A8) $$

where $\{\psi_i^{(j)}; j = 1, \dots, N\}$ is a sample from the posterior distribution of $M_i$ (see Kass and Raftery, 1993, p. 19, eqn 12).

For this application, the $g$ function was chosen as a multinormal distribution with mean equal to the empirical mean and covariance matrix equal to the estimated covariance matrix of our Metropolis-Hastings sample, respectively.

## Appendix C: Prior sensitivity analysis

To study the robustness of our results, we conducted a prior sensitivity analysis by varying prior information for model parameters ($\mu, \xi, \sigma^2$). We explored the range of priors from the marginal informative priors used in the paper core up to vague (poorly informative) priors as much as possible. In the case of $M_0$ model, marginal vague priors could be tested for all the three parameters instead of informative ones. For $M_1$ model, however, the informative prior should be always kept for $\xi$, because, if for this parameter a vague prior was used, the MH algorithm convergence could not be achieved. This is explained by the high correlation between the $\xi$ and $\phi$ parameters (close to 0.5, see Fig. 6). For the $\phi$ parameter, a uniform distribution $U(-1, 1)$ was used in all the analyses. The different combinations of informative and non-informative priors tested are shown in Table 6. Corresponding posterior distributions are shown in Figure 9.

<table><caption>**Table 6.** Marginal prior distributions used for the sensitivity analysis. Priors 3 and 5 are used with $M_0$ model only, and prior 6 with $M_1$ model only. $\Gamma$ denotes the Gamma distribution and $\mathrm{Inv}\Gamma$ denotes the inverse Gamma distribution. Prior1 is the one used in the paper core. ‘,’ denotes vague marginal priors</caption><thead><tr><th>Parameter</th><th>μ</th><th>ξ</th><th>σ<sup>2</sup></th></tr></thead><tbody><tr><td>Prior1</td><td>Γ(9, 0.03)</td><td>Γ(25, 40)</td><td>InvΓ(3, 20)</td></tr><tr><td>Prior2</td><td>-</td><td>Γ(25, 40)</td><td>InvΓ(3, 20)</td></tr><tr><td>Prior3</td><td>Γ(9, 0.03)</td><td>-</td><td>InvΓ(3, 20)</td></tr><tr><td>Prior4</td><td>Γ(9, 0.03)</td><td>Γ(25, 40)</td><td>-</td></tr><tr><td>Prior5</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Prior6</td><td>-</td><td>InvΓ(3, 20)</td><td>-</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Prior sensitivity analysis. Post *k*, where *k* ranges from 1 to 6, are the posterior distributions obtained with the priors 1–6 of Table 6. Prior1 is the informative prior used in the paper core. Panels (a-c): *M*₀ model; panels (d-g): *M*₁ model.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Images used in the photogrammetric process. (a) and (b) are the left and right images for the avalanche released on 13 February 2013 for *t* = 20 s after triggering. This couple is used for restitution and to map the location of the avalanche head. Markers plot the permanent ground control points used for image orientation. To illustrate the spatial extension, the coordinates and elevations (in metres) are indicated for the upper and lower control points (*m*3, *t*6), and for control point *m*5 at the location of Col du Galibier road. The right image (c) has been taken after the avalanche stops. Temporary control points *t*1 to *t*6 are setup after the avalanche in the runout area to improve image orientation in this area according to the avalanche deposit.</figcaption></figure>

Even if slight differences between the different posteriors can be noted as
one replaces one prior by another, overall, for both models, results remain
quite similar, leading in all cases to posterior estimates very close to the
ones obtained with the informative priors used in the paper core (see
Table 4). Also, the common logarithm of the Bayes Factors was calculated
between all combinations of M₁ and M₀ posterior samples and, in nearly
all cases, there is strong to decisive evidence to reject M₀ compared to M₁
model, as demonstrated in the paper core with informative priors. Indeed,
the obtained log₁₀ (B₁₀) values are between 0.57 and 2.49 with a mean of
1.38 and a standard deviation of 0.61. Hence, all in all, posterior inference

and model selection show little sensitivity to prior specification. Since it is
advised to use priors with a finite domain as much as possible to avoid spuri-
ous results (Kuczera and Parent, 1998) we kept the informative priors in the
paper core.

**Appendix D: Images used in the photogrammetric process**

In this appendix, some examples of images used in the photogrammetric
process are shown (Fig. 10).