Introduction
The complexity of snow avalanche physics is related to the variability and changing nature of snow (Schweizer and others, Reference Schweizer, Bruce Jamieson and Schneebeli2003; Ancey, Reference Ancey2006; Castebrunet and others, Reference Castebrunet, Eckert and Giraud2012; Steinkogler and others, Reference Steinkogler, Sovilla and Lehning2014). Evidence obtained at full-scale experimental slopes (Sovilla and others, Reference Sovilla, Schaer, Kern and Bartelt2008; Vriend and others, Reference Vriend2013; Prokop and others, Reference Prokop2015; Faug and others, Reference Faug, Turnbull and Gauer2018) shows a myriad of avalanche propagation and stopping regimes (Köhler and others, Reference Köhler, McElwaine, Sovilla, Ash and Brennan2016), and numerical propagation models can reproduce these observations with increasing realism (Bartelt and others, Reference Bartelt, Buser, Vera Valero and Bühler2016; Gaume and others, Reference Gaume, Gast, Teran, van Herwijnen and Jiang2018). However, knowledge concerning the mechanical behaviour of snow during motion and associated processes (granulation, erosion/deposition, etc.) remains incomplete (Steinkogler and others, Reference Steinkogler, Gaume, Löwe, Sovilla and Lehning2015; Truong and others, Reference Truong, Keylock, Eckert, Bellot and Naaïm2018). From a macroscopic point of view, experimental approaches (Casassa and others, Reference Casassa, Narita and Maeno1989; Rognon and others, Reference Rognon2008; Kern and others, Reference Kern, Bartelt, Sovilla and Buser2009) and the proposal of Voellmy (Reference Voellmy1964) suggest rheological behaviours which remain ad hoc. This renders on-site calibration on the basis of local data unavoidable (Ancey and Meunier, Reference Ancey and Meunier2004; Salm, Reference Salm2004; Eckert and others, Reference Eckert2012) to, e.g. predict high-return-period avalanches in land use planning and assess the related risk (Keylock and others, Reference Keylock, McClung and Magnússon1999; Meunier and others, Reference Meunier, Ancey and Taillandier2004; Favier and others, Reference Favier, Eckert, Bertrand and Naaim2014b,Reference Favier, Bertrand, Eckert and Naaima). 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, Reference Borstad and McClung2009; Fischer, Reference Fischer2013).
After deterministic inversion methods had shown their limits (Dent and Lang, Reference Dent and Lang1980; Dent and others, Reference Dent1998; Ancey and others, Reference Ancey, Meunier and Richard2003), 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, Reference Oakley and O'Hagan2002; Carmassi and others, Reference Carmassi, Barbillon, Keller, Parent and Chiodetti2018), 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, Reference Ancey2005; Straub and Grêt-Regamey, Reference Straub and Grêt-Regamey2006; Eckert and others, Reference Eckert, Parent and Richard2007, Reference Eckert, Parent, Naaim and Richard2008, Reference Eckert, Parent, Faug and Naaim2009, Reference Eckert, Naaim and Parent2010; Schläppy and others, Reference Schläppy2014). Specifically, Gauer (Reference Gauer, Medina-Cetina, Lied and Kristensen2009) used a Bayesian framework to calibrate the friction parameters of three avalanche sliding block models and Fischer and others (Reference Fischer, Fromm, Gauer and Sovilla2014) proposed a method to evaluate simulations compared to Doppler radar observations. Ultimately, Naaim and others (Reference Naaim, Durand, Eckert and Chambon2013) 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, Reference Fischer, Kofler, Fellin, Granig and Kleemayr2015). 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, Reference Kuczera1983; McInerney and others, Reference McInerney, Thyer, Kavetski, Lerat and Kuczera2017; Schaefli and Kavetski, Reference Schaefli and Kavetski2017; Sun and others, Reference Sun, Yuan and Liu2017). 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, Reference Kuczera and Parent1998; Evin and others, Reference Evin, Thyer, Kavetski, McInerney and Kuczera2014), 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, Reference Thibert2015). 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, Reference Turnbull and Bartelt2003; Gauer, Reference Gauer2014; Dreier and others, Reference Dreier, Bühler, Ginzler and Bartelt2016), 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, Reference Naaim, Naaim-Bouvet, Faug and Bouchet2004; Christen and others, Reference Christen, Kowalski and Bartelt2010; Bartelt and others, Reference Bartelt, Bühler, Buser, Christen and Meier2012). 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, Reference Ancey and Meunier2004; Gauer, Reference Gauer2014). However, for hazard mitigation, simple models with few parameters remain useful (Salm, Reference Salm2004) 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:
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, Reference Voellmy1964). This means that the friction force is:
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, Reference Ancey, Meunier and Richard2003). 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, Reference Ancey and Meunier2004; Eckert and others, Reference Eckert, Parent and Richard2007).
Statistical model formulation
Let us denote f our avalanche propagation model. The model f predicts the avalanche speed (m · s−1) 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 1, …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:
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 σ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 ε 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:
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, Reference Fisher1922). 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, Reference Bayes1763; Bernardo and Smith, Reference Bernardo and Smith2009). Hence, Bayesian statistics is now widely accepted as a reasonable option in environmental sciences (Berliner, Reference Berliner2003; Clark, Reference Clark2005) and we use this framework in what follows.
For simplicity, let us denote γ, the set of parameters related to the errors, it means γ = {σ2} for ${\cal M}_{0}$
and γ = {ϕ, σ2} 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:
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:
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):
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 (Reference Barbolini and Savi2001), Borstad and McClung (Reference Borstad and McClung2009), Buhler and others (Reference Buhler, Argue, Jamieson and Jones2018). 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, Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller1953; Torre and others, Reference Torre, Boreux and Parent2001). 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 (Reference Robert2015). 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 ψ (0) and a proposal distribution q. Each iteration k of the algorithm consists in:
1. Generating ψ′ ~ q(.|ψ (k−1)).
2. Calculating $u \sim {\cal U}\lpar 0\comma\; 1\rpar$
.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 (Reference Robert2015), 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, Reference Gelman and Rubin1992).
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, Reference Kass and Raftery1995):
where
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 log10B 10 is proposed in the work of Kass and Raftery (Reference Kass and Raftery1995), where log10 denotes the common logarithm. Notably, a log10B 10 value between 1 and 2 indicates strong evidence in favour of rejecting ${\cal M}_0$
compared to ${\cal M}_1$
, and log10B 10 > 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, Reference Thibert2015). 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−3, ranging between 270 and 225 kg · m−3 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, Reference Fierz2009).
Avalanche front positions were determined using a high rate photogrammetric system that was specifically developed (Soruco and others, Reference Soruco, Thibert, Vincent, Blanc and Héno2011; Thibert and others, Reference Thibert2015). 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, Reference Faig, Shih and Liu1990). 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−6 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, Reference Prokop2008; Prokop and others, Reference Prokop2015), 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, Reference Pulfer, Naaim, Thibert and Soruco2013):
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:
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, Reference Thibert2015). 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 (Reference Bühler, Christen, Kowalski and Bartelt2011) 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 μ, ξ, σ2 and ϕ parameters of ${\cal 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 $\sigma ^2_{\rm {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.
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, Reference Gilks, Richardson and Spiegelhalter1995; Eckert and others, Reference Eckert, Naaim and Parent2010). 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 (Reference Salm, Burkard and Gubler1990). 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, Reference Kuczera and Parent1998; Torre and others, Reference Torre, Boreux and Parent2001; Robert and Casella, Reference Robert and Casella2009; Eckert and others, Reference Eckert, Naaim and Parent2010). The R Core Team (2019) package coda created by Plummer and others (Reference Plummer, Best, Cowles and Vines2006) 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 (H0: 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 (Reference Robert2015) reach fast convergence in a random walk. Finally, the Gelman–Rubin convergence criterion was computed (Gelman and Rubin, Reference Gelman and Rubin1992). 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, Reference Robert and Casella2009). 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 (Reference Dent and Lang1980). 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 σ2 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 m2 · s−1 with model ${\cal M}_1$
and 410 m2 · s−1 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 σ2 from model ${\cal M}_{0}$
to model ${\cal M}_{1}$
(Table 4). Indeed, the errors variance σ2 is higher if the autocorrelation ϕ is not included. Specifically, the standard deviation of the errors decreases from 3.6 to 3.1 m · s−1 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 (H0: 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 H0: 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:
and, the second one:
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:
By using a Bayesian approach, we obtained a sample of parameters (μ, ξ, σ 2) 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, Reference Ancey and Meunier2004; Gauer, Reference Gauer2014) 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:
and the second:
Obtained values are Δ1 = 1.31 m · s−1 and Δ2 = 1.72 m · 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 m and x = 250 m (plateau phase), around 15 m · s−1 with ${\cal M}_{0}$
. The underestimation with ${\cal M}_{0}$
with regards to ${\cal M}_{1}$
can be estimated by the ratio between Δ2 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 log10 B10 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 (Reference Kass and Raftery1995, 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, log10B 10 > 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 σ 2 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 σ 2 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 log10B 10 obtained for the 100 synthetic avalanches. The value of 2 in green corresponds to the reference value of Kass and Raftery (Reference Kass and Raftery1995) 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.
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.
Appendix A: Likelihood of an AR1 process
Since |ϕ| < 1, the ε t error has a representation in terms of the innovations η t (also known as MA(∞) representation, see Gouriéroux and Monfort, Reference Gouriéroux and Monfort1995 for more details):
where we recall that ϕ and η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 ε 0 with zero mean and variance equal to σ 2/1 − ϕ 2, independent from (η1, …, ηt), and we get that:
The marginal distribution of ε t is Gaussian because it is the sum of independent and identically normally distributed η t. In particular the first error term ε 1 is Gaussian:
Considering that ε 2 = ϕε 1 + η 2, the conditional distribution of ε 2 given ε 1 is $\epsilon _{2}\vert \epsilon _{1} \sim {\cal N}\lpar \phi \epsilon _{1}\comma\; \sigma ^2\rpar$
. In a more general way:
Then, the joint distribution of errors is:
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:
where i ∈ {0, 1}.
If g is a density function defined on an ensemble Ω:
From Monte Carlo simulations, ${1}/{{p}\lpar v_{{\rm {obs}}}\vert {\cal M}_{i}\comma\; x\rpar }$
can thus be estimated as:
where $\lcub \psi _i^{\lpar j\rpar }\semicolon\; j=1\comma\; \ldots \comma\; N \rcub$
is a sample from the posterior distribution of ${\cal M}_{i}$
(see Kass and Raftery, Reference Kass and Raftery1993, 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 (μ, ξ, σ 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 ${\cal M}_0$
model, marginal vague priors could be tested for all the three parameters instead of informative ones. For ${\cal M}_1$
model, however, the informative prior should be always kept for ξ, 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 ξ and ϕ parameters (close to 0.5, see Fig. 6). For the ϕ parameter, a uniform distribution ${\cal U}\lpar -\!1\comma\; 1\rpar$
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.
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): ${\cal M}_0$
model; panels (d–g): ${\cal M}_1$
model.
Table 6. Marginal prior distributions used for the sensitivity analysis. Priors 3 and 5 are used with ${\cal M}_0$
model only, and prior 6 with ${\cal M}_1$
model only. Γ denotes the Gamma distribution and InvΓ denotes the inverse Gamma distribution. Prior1 is the one used in the paper core. ‘-’ denotes vague marginal priors
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 ${\cal M}_1$
and ${\cal M}_0$
posterior samples and, in nearly all cases, there is strong to decisive evidence to reject ${\cal M}_0$
compared to ${\cal M}_1$
model, as demonstrated in the paper core with informative priors. Indeed, the obtained log10 (B 10) 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 spurious results (Kuczera and Parent, Reference Kuczera and Parent1998) 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).
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 (m3, t6), and for control point m5 at the location of Col du Galibier road. The right image (c) has been taken after the avalanche stops. Temporary control points t1 to t6 are setup after the avalanche in the runout area to improve image orientation in this area according to the avalanche deposit.
Open access
