Hostname: page-component-7c8c6479df-7qhmt Total loading time: 0 Render date: 2024-03-29T12:32:57.329Z Has data issue: false hasContentIssue false

Multivariate parameter optimization for computational snow avalanche simulation

Published online by Cambridge University Press:  10 July 2017

Jan-Thomas Fischer*
Affiliation:
Department of Natural Hazards, Austrian Research Centre for Forests (BFW), Innsbruck, Austria
Andreas Kofler
Affiliation:
Department of Natural Hazards, Austrian Research Centre for Forests (BFW), Innsbruck, Austria
Wolfgang Fellin
Affiliation:
Division of Geotechnical and Tunnel Engineering, Institute of Infrastructure, University of Innsbruck, Innsbruck, Austria
Matthias Granig
Affiliation:
Snow and Avalanche Center, Avalanche and Torrent Control (WLV), Innsbruck, Austria
Karl Kleemayr
Affiliation:
Department of Natural Hazards, Austrian Research Centre for Forests (BFW), Innsbruck, Austria
*
Correspondence: Jan-Thomas Fischer <jt.fischer@uibk.ac.at>
Rights & Permissions [Opens in a new window]

Abstract

Snow avalanche simulation software is a commonly used tool for hazard estimation and mitigation planning. In this study a depth-averaged flow model, combining a simple entrainment and friction relation, is implemented in the software SamosAT. Computational results strongly depend on the simulation input, in particular on the employed model parameters. A long-standing problem is to quantify the influence of these parameters on the simulation results. We present a new multivariate optimization approach for avalanche simulation in three-dimensional terrain. The method takes into account the entire physically relevant range of the two friction parameters (Coulomb friction, turbulent drag) and one entrainment parameter. These three flow model parameters are scrutinized with respect to six optimization variables (runout, matched and exceeded affected area, maximum velocity, average deposition depth and mass growth). The approach is applied to a documented extreme avalanche event, recorded in St Anton, Austria. The final results provide adjusted parameter distributions optimizing the simulation–observation correspondence. At the same time, the degree of parameter–variable correspondence is determined. We show that the specification of optimal values for certain model parameters is near-impossible, if corresponding optimization variables are neglected or unavailable.

Type
Research Article
Creative Commons
Creative Common License - CCCreative Common License - BY
Copyright © International Glaciological Society 2015 This is an Open Access article, distributed under the terms of the Creative Commons Attribution license. (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.
Copyright
Copyright © International Glaciological Society 2015

Introduction

Snow avalanche simulation tools are used for hazard estimations and protection planning (Reference Sampl and ZwingerSampl and Zwinger, 2004; Reference Christen, Bartelt and KowalskiChristen and others, 2010a). Initial conditions and flow model parameters have to be defined carefully to obtain meaningful simulation results. The challenges arising in this process are manifold. Existing avalanche dynamics models contain parameters some of which are more conceptual than physical. A direct calibration, i.e. measuring the specific material parameters in the field directly, is hardly possible. Existing guidelines provide model parameter suggestions or estimates for release conditions for extreme avalanches (Reference Salm, Burkhard and GublerSalm and others, 1990; Reference Maggioni and GruberMaggioni and Gruber, 2003; Reference Gruber and BarteltGruber and Bartelt, 2007). They are used to quantify avalanche danger for hazard scenarios with different return periods. Back calculation of real avalanche events requires modification of these parameter suggestions to include physical processes such as flow regime transitions (Reference Issler and GauerIssler and Gauer, 2008; Reference Bartelt, Bühler, Buser, Christen and MeierBartelt and others, 2012) as well as the effects of snow temperature and entrainment (Reference Naaim, Durand, Eckert and ChambonNaaim and others, 2013; Reference Vera Valero, Wikstroem Jones, Bühler and BartelVera Valero and others, 2015). Only then is it possible to reproduce observed runout, flow velocities, impact pressures and deposition depths. However, at present it is still not clear whether process-oriented models can be used in engineering practice since they demand the specification of detailed snow-cover and release information. Therefore there is still an urgent need to find an optimization method to establish model parameters for hazard scenarios solved with existing simulation software in avalanche practice.

Early optimization approaches date back to Reference Dent and LangDent and Lang (1980, Reference Dent and Lang1983), who utilized velocity and deposition data from single snow-flow test experiments to fit model parameters. This work was followed by, for example, Reference Ancey, Meunier and RichardAncey and others (2003), who introduced a deterministic inversion method, finding a parameter distribution for a simplified flow model, reproducing avalanche runouts. Besides deterministic approaches, there has been a strong development of probabilistic (especially Bayesian) frameworks originating from the hydrological community. These developments were accompanied by debates concerning different (formal and informal) optimization strategies (Reference Beven and BinleyBeven and Binley, 1992; Reference Vrugt, Ter Braak, Gupta and RobinsonVrugt and others, 2008). The main issues concern the arising quantification of uncertainty (Reference Montanari, Shoemaker and Van de GiesenMontanari and others, 2009).

In the field of snow and avalanches, several authors (Reference AnceyAncey, 2005; Reference Eckert, Parent and RichardEckert and others, 2007, Reference Eckert, Parent, Naaim and Richard2008; Reference Gauer, Medina-Cetina, Lied and KristensenGauer and others, 2009) have employed Bayesian techniques to solve the inverse problem for lumped mass propagation models and avalanche runout, by, for example, analyzing occurrence frequencies for multiple avalanche paths. Reference Naaim, Durand, Eckert and ChambonNaaim and others (2013) employed this method to link model parameters to physical properties of snow. However, challenges in the model validation arise, because formal Bayesian approaches explicitly consider model input and parameter uncertainty (Reference Vrugt, Ter Braak, Gupta and RobinsonVrugt and others, 2008). Therefore the individual error sources have to be identified and quantified in order to assign the resulting probability to one of them (Reference Beven, Smith and FreerBeven and others, 2008; Reference Vrugt, Ter Braak, Gupta and RobinsonVrugt and others, 2008). This problem is especially severe with regard to avalanches because of the interaction of input and parameter uncertainty (Reference Eckert, Naaim and ParentEckert and others, 2010). For example, Reference AnceyAncey (2005) showed the dependence between a frictional parameter and input parameters such as avalanche volume and snow properties for a sliding-block model.

Informal statistical approaches, to which the proposed method is similar, do not explicitly consider model uncertainties (Reference McMillan and ClarkMcMillan and Clark, 2009). They are based on a more arbitrary function to quantify the correspondence between simulation results and observation. With this objective measure, adjusted parameter distributions can be computed. At the same time, these resulting distributions represent an estimate of total uncertainty. Despite the differences between informal and formal Bayesian approaches it has been found that they can yield very similar estimates of the total result uncertainty (Reference Vrugt, Ter Braak, Gupta and RobinsonVrugt and others, 2008).

The application of optimization methods with simulation tools operating in three-dimensional (3-D) terrain is limited. Most studies of avalanche simulation tools are based on multi-parameter models, but they have mostly been optimized for single optimization variables, namely the avalanche runout (Reference AnceyAncey, 2005; Reference Gruber and BarteltGruber and Bartelt, 2007; Reference BozhinskiyBozhinskiy, 2008). Providing single constraints is insufficient to obtain a unique multivariate parameter set. Information on flow variables (e.g. velocity) is rarely accessible and is therefore used in few case studies (Reference Sailer, Rammer and SamplSailer and others, 2002; Reference Ancey and MeunierAncey and Meunier, 2004; Reference IsslerIssler and others, 2005; Reference Gauer, Medina-Cetina, Lied and KristensenGauer and others, 2009; Reference Fischer, Fromm, Gauer and SovillaFischer and others, 2014). In cases where no information is available, empirical analysis can provide rough estimates for missing measurements. For example, the analyses of Reference McClung and SchaererMcClung and Schaerer (2006) allow us to estimate the maximum avalanche velocity by scaling it with the total fall height along the avalanche path. The estimate is based on basic energy relations, recently confirmed for a variety of extreme avalanches (Reference GauerGauer, 2014).

The main focus of the presented optimization concept is to provide adjusted parameter distributions employing a systematic, multivariate comparison of simulation results with field observations and their related uncertainties. The proposed framework is sketched in Figure 1. It is tested for a catastrophic avalanche. A simple, three-parameter flow model including entrainment is employed and implemented in the operational snow avalanche simulation software, SamosAT (Snow Avalanche MOdelling and Simulation – Advanced Technology) (Reference Zwinger, Kluwick, Sampl, Hutter and KirchnerZwinger and others, 2003; Reference Sampl and ZwingerSampl and Zwinger, 2004). A large number of simulation runs is performed following an input parameter distribution. The results are analyzed in an avalanche path dependent coordinate system (Reference FischerFischer, 2013). The multivariate parameter optimization is carried out with respect to three varying model parameters and six different optimization variables, enabling the quantification of simulated and observed avalanche characteristics. By introducing a selection rule, parameter combinations with optimal simulation–observation correspondence are identified. The main results are problem-suited parameter distributions. These adjusted distributions provide peak values for the flow model parameters leading to an optimized simulation result and provide a base for future guidelines.

Fig. 1. The proposed framework consists of three parts: observation, simulation and optimization. Monte Carlo simulation runs are performed for a large number of parameter combinations Θ = μ, ξ, eb following the input parameter distributions . The simulation software, SamosAT, is utilized, taking into account boundary and initial conditions, which are obtained from the observations. The set of optimization variables X = {…} (with variables runout r, matched affected area (true) T, exceeded affected area (false) F, maximum velocity u max, average deposition depth d and mass growth G) is defined in terms of the simulation results (flow depths h, velocities u and mass m of the avalanche) as well as observations and their related uncertainty . The simulation–observation correspondence α(Θ) is quantified for each simulation run. Based on a consistently defined correspondence limit α lim, parameter combinations are withdrawn or accepted, yielding the final results: adjusted parameter distributions ΩΘ.

Simulation Concept

The simulation concept comprises the choice of

  1. 1. the simulation software including the physical flow model and its numerical implementation (simulation

  2. approach)

  3. 2. the initial and boundary conditions, including the digital elevation model (DEM), initial distribution of snow and the model parameter settings (simulation input)

  4. 3. the survey and interpretation of simulation results in view of the model evaluation, and the comparison to the field observations (simulation output).

Simulation approach

Simulation software for the dense, most destructive part of snow avalanches is based on two-dimensional depth-averaged, deterministic flow models (Reference Savage and HutterSavage and Hutter, 1989; Reference Naaim, Furdada and MartinezNaaim and others, 2002; Reference Pitman, Nichita, Patra, Bauer, Sheridan and BursikPitman and others, 2003; Reference Sampl and ZwingerSampl and Zwinger, 2004; Reference Christen, Kowalski and BarteltChristen and others, 2010b; Reference Mergili, Schratz, Ostermann and FellinMergili and others, 2012; Reference PudasainiPudasaini, 2012). An adequate mathematical description requires switching between different coordinate systems. The following model equations will be written in a Lagrangian framework using the notation for flow variables depth and velocity and the mountain surface in a natural coordinate system moving with the avalanche. In the Eulerian framework, variables and surface are denoted h, u, z with respect to the coordinate system aligned with the avalanche path. The mountain surface is represented by a DEM with a spatial resolution of 5 m × 5 m, which is assumed to represent the winter, snow-covered surface sufficiently. Reference Sampl and ZwingerSampl and Zwinger (2004) presented a Lagrangian formulation of the mass and momentum balance, describing the spatio-temporal evolution of the primary variables: depth-averaged flow depth and velocity (Reference IversonIverson, 2012). Here the equations are formulated for an incompressible, isotropic material, with a general basal shear stress τb and an entrainment rate , integrated over an infinitesimal control volume , that moves with the avalanche (Reference Zwinger, Kluwick, Sampl, Hutter and KirchnerZwinger and others, 2003). This yields a locally orthogonal coordinate system. i = 1 is in the direction of the surface-parallel velocity vector, i = 2 is surface-parallel and orthogonal to the velocity vector and i = 3 appears naturally normal to the surface :

(1)

(2)

with the components of the gravitational acceleration gi , the surface-parallel velocity components u 1, ū 2, the surface-normal flow depth and the resulting normal stress . The term accounts for the change in the normal acceleration due to surface curvature in the flow direction (Reference Gray, Wieland and HutterGray and others, 1999; Reference Pudasaini and HutterPudasaini and Hutter, 2003; Reference Fischer, Kowalski and PudasainiFischer and others, 2012). The first term on the right-hand side of Eqn (2) accounts for the acceleration due to gravity. The second term arises due to the pressure gradients on the control volume V, with boundary line ∂A with elements dl and the normal vector ni (Reference Zwinger, Kluwick, Sampl, Hutter and KirchnerZwinger and others, 2003). The third term describes the frictional decelerations opposing the direction of movement i = 1, with the Kronecker delta δij . The last term arises due to the modified mass balance (right-hand side of Eqn (1)) and causes a momentum loss if additional mass has to be accelerated by the avalanche.

Equations (1) and (2) are solved with a smoothed particle hydrodynamics (SPH) scheme (Reference MonaghanMonaghan, 1992) for the three variables ū = (ū 1, ū 2) and depth , by discretization of the released avalanche volume in a large number of mass elements. The number of mass elements is calculated in accordance with the claim that the mass per numerical particle is ∼2000 kg (cf. Reference Sampl and ZwingerSampl and Zwinger, 2004). The simulation end time t end was chosen carefully according to the criterion that the pressure isoline of p = p lim showed no significant changes over time, which was sufficiently reached with t end = 350 s for the specific avalanche simulations (Reference FischerFischer, 2013; Reference Teich, Fischer, Feistl, Bebi, Christen and Grêt-RegameyTeich and others, 2014). The total duration of a computation is in the order of several minutes with a standard computer.

Basal shear stress τ b and entrainment rate

Over the years many different (mostly phenomenological) friction and entrainment relations have been implemented in different flow models (Reference HarbitzHarbitz, 1998). The goal of this work is not to discuss or compare different approaches but to show that a systematic, multivariate parameter optimization is possible. Therefore the well-known Voellmy friction relation for the basal shear stress and a simple assumption for the entrainment rate are employed. The basal shear stress τ b combines a Coulomb bottom friction with a velocity-dependent drag term:

(3)

with dimensionless friction parameter μ and turbulent friction coefficient ξ(m s−2) (Reference VoellmyVoellmy, 1955).

The entrainment rate,

(4)

is proportional to the bottom shear stress τ b, which is similar to other definitions found in the literature (Reference Christen, Kowalski and BarteltChristen and others, 2010b), and includes the phenomenological parameter eb (m2 s−2), which allows for interpretation as specific erosion energy. For small parameter values eb → 0 the entrainment rate is very large, ; however, independently of the entrainment rate, the amount of entrained snow is limited by the available snow reservoir q < h snow (cf. Eqn (5) below). Entraining the entire snow reservoir h snow at the flow front corresponds to the process of frontal plowing. Larger eb values allow for a gradual erosion, from the front to the tail of the avalanche (Reference Gauer and IsslerGauer and Issler, 2004). For large parameter values eb → ∞ the entrainment rate diminishes , i.e. no snow is entrained.

Simulation input

To perform snow avalanche simulations, parameter set-up for the employed flow model and initial conditions have to be defined.

Initial conditions

For the presented analysis, initial conditions are assumed to be constant and are estimated through direct measurement or empirical methods. Release areas are either delineated by direct event observation or a set of empirical rules, which are mainly based on slope and planar curvature (Reference Maggioni and GruberMaggioni and Gruber, 2003; Reference Bühler, Kumar, Veitinger, Christen and StoffelBühler and others, 2013). Typically the initial distribution of released snow is either directly measured or estimated by means of an extreme snowfall. The estimated snow depth h 0 is often linked to the sum of new snow over 3 days for a certain return period (Reference Burkard and SalmBurkard and Salm, 1992) measured on flat ground at a reference altitude z 0.

Here we estimate a smooth initial snow-cover distribution h snow, assuming equal precipitation at each location. This approach allows us to determine the initially released snow mass and the potentially erodible snow mass in a consistent manner. Precipitation varies with altitude through the snow depth–altitude gradient Δh. The influence of wind transport is neglected:

(5)

θ is the local slope angle. The snow depth gradient Δh is a regional coefficient and varies for different precipitation characteristics (Reference Burkard and SalmBurkard and Salm, 1992). The smooth snow distribution leads to lower (higher) snow depth on steep (flat) slopes. Once the avalanche is released, the snow reservoir h snow is depleted inside the release area and evolves over time. for the rest of the mountain surface.

Reference event

Documentation of extreme events characterizing processes in terms of release conditions, flow path and runout zone provides important information for performing hazard assessment and model optimization. The Wolfsgruben avalanche path starts in a release area of ∼20 ha, with mean slope angle 36.5° at ∼2244 m a.s.l. It follows a gully, with a width of ∼100 m, and finally reaches the community of St Anton a.A., Austria (at ∼1260 m a.s.l.; Fig. 2). On 13 March 1988 a catastrophic avalanche struck the village, affecting an area of ∼6.5 ha with mean slope angle of 14.5°. The avalanche led to severe loss of life and property. Three houses and nine cars were destroyed, and several other buildings, about 20 cars and existing infrastructure were damaged (back-calculated pressures range between about 7.5 and 17.5 kPa). Several people were killed or injured. The event return period has been estimated to be sufficiently large to serve as a design event (>150 years in Austria; cf. Reference Johannesson, Gauer, Issler and Liedjohannesson and others, 2009).

Fig. 2. Wolfsgruben avalanche (images WLV). (a) Wolfsgruben avalanche path (vertical drop Δz = 984 m): central flowline in black, release area (196 225 m2, with mean slope angle 36.5°) in blue, affected area (64 153 m2, at 14.5°) in orange. (b) Destroyed house in the runout area. Back-calculated pressures of damages on infrastructure range from about 7.5 to 17.5 kPa.

The observations allow us to reconstruct the initial snowcover distribution h snow of the event, with reference snow height of h 0 = 1.61 m at z 0 = 1289 m a.s.l. and a snow-depth–altitude gradient of Δh = 8cm (100 m)−1 for the Arlberg region. Considering Eqn (9) and the given topography this information leads to a total release volume of V release ≈ 354 600 m3. The snow reservoir depth for entrainment ranges between about 0.6 and 1.8 m. Besides basic observations of snow avalanches (e.g. the delineation of release areas, affected path and runout zone), information on physical properties like deposition depth (∼4 m) and density (∼400 kg m−3), as well as flow velocity, is essential for the parameter optimization. This information allows us to define optimization variables (e.g. runout, maximum velocity or mass growth), which are related to the avalanche (Table 1; next section).

Table 1. Observational variables for the Wolfsgruben avalanche

Flow model parameters Θ = {μ, ξ, eb }

For the presented optimization the input parameter distributions, given by , yield simulation samples of size N. For each simulation run, a set of flow model parameter combinations Θ n = {μn , ξn , eb , n } with n = 1, … , N is assigned. The number of simulation runs N in a sample has to be sufficiently large that the obtained results are statistically significant and stable (N ≳ 8000; cf. Analysis subsection below, ‘investigating the statistical significance’). Plausible combinations of input parameters for the deterministic flow model are obtained through a Latin hypercube sampling (Reference SteinStein, 1987). This sampling method provides a probabilistic representation of the input distributions dividing the range of each variable into equally probable intervals. The flow model parameters Θ = {μ, ξ, eb } are assumed independent. Naturally the random samples include a certain degree of correlation. Since independence in the input samples is desirable, a correlation control is applied on the initial parameter sample (Reference Oberguggenberger, King and SchmelzerOberguggenberger and others, 2009; Reference Fischer, Fromm, Gauer and SovillaFischer and others, 2014). With this empirical method the parameter samples are rearranged in order to get closer to parameter independence (i.e. close to a diagonal rank correlation matrix; Reference Iman and ConoverIman and Conover, 1982).

The input parameter distributions are specified by their distribution function and the related parameter range. Since the estimation of the adjusted distribution for a specific event is a goal of the analysis, we assume an equal distribution function for the input samples, which leaves the parameter ranges to be specified. Assigning the parameter ranges too small may exclude possible solutions; defining the ranges too large multiplies the computational efforts. The interval bounds may be constrained by the physically relevant parameter space, values found in the literature or results of experimental work.

For example, considering the limits μ → 0, ξ → ∞ or eb → ∞, the effect of Coulomb friction, turbulent friction or entrainment is negligible. In contrast, the limits μ → ∞, ξ → 0 lead to infinite friction, i.e. the avalanche cannot move; eb → 0 corresponds to entrainment by frontal plowing, which may initiate uncontrolled mass gain, often accompanied and identifiable by unrealistic lateral spreading of the avalanche. In practice, the limits do not coincide with 0 or ∞ but may be determined by scaling analysis of the respective terms in the model equations. Taking into account values of back calculations for similar flow models allows us to further specify the parameter ranges. Reference VoellmyVoellmy (1955) estimated values in the range μ = 0.08–0.15 and ξ = 400–600 m s−2. Reference GublerGubler (1987) used velocity data, finding different values for the Coulomb friction in the path and the runout zone, in the range μ = 0.15–0.5. For practical purposes different authors have recommended different values. For example, Reference Buser and FrutigerBuser and Frutiger (1980), Reference Salm, Burkhard and GublerSalm and others (1990) and Reference Gruber and BarteltGruber and Bartelt (2007) (often referred to as Swiss Guidelines) proposed Coulomb friction values μ = 0.155–0.3 and turbulent friction coefficients ξ = 400–1000 m s−2 depending on different variables such as return period, avalanche size or terrain features. Reference Barbolini, Gruber, Keylock, Naaim and SaviBarbolini and others (2000) found Coulomb friction values between 0.13 and 0.4 and turbulent friction coefficients between 1000 and 4500 m s−2 for different models and sites. Reference Naaim, Faug, Naaim and EckertNaaim and others (2010) fixed ξ to 1500 m s−2, obtaining μ values between 0.15 and 0.35 based on historical data from the avalanche path Taconnaz. In a similar framework, Reference Naaim, Durand, Eckert and ChambonNaaim and others (2013) scanned values for the static Coulomb friction from 0.1 to 0.7 and turbulent friction from 500 to 1500 m s−2. However, with some exceptions (Reference Sovilla, Margreth and BarteltSovilla and others, 2007; Reference BozhinskiyBozhinskiy, 2008; Reference Naaim, Durand, Eckert and ChambonNaaim and others, 2013), snow entrainment is disregarded or not explicitly treated. Parameter specifications derived by experiments fit in similar ranges; for example, experiments performed on snow chutes by Reference Tiefenbacher and KernTiefenbacher and Kern (2004) and Reference Platzer, Bartelt and KernPlatzer and others (2007) allowed estimation of an effective Coulomb friction of μ = 0.22–0.72. However, in the present context, friction coefficients appear more conceptual than physical.

Based on this, we specify the following parameter range specifications for the equally distributed input parameters , ξ =[400, 15 000 m s−2 and eb = [0, 75 000] m2 s−2, with minimal spacing Δμ = 0.01, Δξ = 50 m s−2, Δeb = 250 m2 s−2.

Simulation results

Primary simulation results are the time evolution of flow depth and velocities. However, the most important simulation results for the evaluation are the peak values, i.e. maximum values over time, of flow depth h, velocity |u| and impact pressure p = ρ |u| 2, with ρ = 200 kg m−3, the density of flowing snow. The formulation of impact pressure is chosen according to guidelines used for avalanche simulations but may differ significantly from this general form (Reference Sovilla, Schaer, Kern and BarteltSovilla and others, 2008). Defining runout or impacted area in terms of impact pressures is in accordance with avalanche hazard mitigation guidelines (Reference Johannesson, Gauer, Issler and LiedJohannesson and others, 2009; Reference FischerFischer, 2013; Reference Teich, Fischer, Feistl, Bebi, Christen and Grêt-RegameyTeich and others, 2014). However, they could equivalently be defined in terms of velocities. The simulation results are evaluated in an avalanche path dependent coordinate system, with flow-path coordinate s and lateral coordinate l, according to the main flow path shown in Figure 2a with a domain width of 500 m (Reference FischerFischer, 2013).

Optimization Variables

The optimization variables represent the different categories that are accessed through the observational data and simulation results. To perform an objective analysis, a set of six optimization variables X = {r, T, F, u max, d, G},

  1. 1. runout r

  2. 2. matched affected area (true) T

  3. 3. exceeded affected area (false) F

  4. 4. maximum velocity u max

  5. 5. average deposition depth d

  6. 6. mass growth G,

is defined in terms of both observation and simulation. Observational variables and their associated uncertainty are denoted by (Table 1), simulation variables by plain X. In this work we do not explicitly account for simulation uncertainties introduced by the deterministic flow model. The information about the source of uncertainty is dispensable for the employed optimization concept. Thus simulation uncertainties can implicitly be associated with the uncertainty of the corresponding observational variable .

Runout r

The definition of avalanche runout in the sense of a simulation and observation is not straightforward, especially in 3-D terrain. Here we define a peak impact pressure related runout in an avalanche-dependent coordinate system. Utilizing a pressure-related measure as runout has several advantages for operational simulation tools (Reference FischerFischer, 2013; Reference Teich, Fischer, Feistl, Bebi, Christen and Grêt-RegameyTeich and others, 2014). For each simulation run, runout r refers to the farthest coordinate s, measured as projected distance in the avalanche-path flow direction, where the maximum value of the peak impact pressure in the cross section still exceeds the predefined pressure limit p lim, i.e. max l p(s, l) > p lim. Here we set p lim =1 kPa, which may be adapted for different hazard-mapping guidelines.

To determine the observational runout the affected area is dealt with in exactly the same manner as the simulation pressure results, assuming peak impact pressures larger than the pressure limit, , inside the observed affected area (Fig. 3). The uncertainty associated with the delineation of the affected area is the range of .

Fig. 3. Sketched simulation results (e.g. simulation outline p lim; blue) and affected area (orange), superimposed with avalanche path domain and coordinate system along the central flowline (black). Runout r, , matched area A matched (green) and exceeded area A exceeded (red) and their complements (blue and gray area), which lead to the optimization variables for true and false prediction T, F, maximum velocity u max, average deposition depth d and mass growth G.

Relative matched and exceeded affected area, T, F

Peak impact pressures serve as a basis to define observed and simulated affected areas, assuming that peak pressures observed in the avalanche exceed the pressure limit p lim, i.e. inside the observed affected area. Comparing the simulated and observed affected areas, four different classifications arise. The four cases include all combinations of matching/non-matching or exceeding/non-exceeding observed area with simulated affected area (cf. the four differently colored areas in Fig. 3; Reference Mergili, Marchesini, Rossi, Guzzetti and FellinMergili and others, 2013). Considering the given affected and total area, two of them are independent. The optimization variables T and F for matching and exceeded areas are defined relative to the affected area and are specified as:

true prediction T = A matching/A affected – simulated area with p > p lim matching observed affected area ,

false prediction F = A exceeded /A affected – simulation area with p > p lim exceeding observed affected area .

It is crucial to consider the main flow direction, represented by the path coordinate s (Reference FischerFischer, 2013), since information on affected areas is mainly available in the runout. To properly account for the false prediction, only areas beyond (in flow direction) the path coordinate s, where , are considered, due to lack of observations in the upstream direction.

The observational value for the optimization variable true prediction T = A matching/A affected is given by , and consequently for the false prediction F = A exceeded/A affected by . For the situation sketched in Figure 3 the respective simulation variable values could be approximated to T ≈ 0.7 and F ≈ 0.05. The associated uncertainty is estimated and expressed relative to the affected area .

Maximum velocity, u max

The maximum velocity u max is defined for each simulation run by taking the maximum of the peak velocities over the entire simulation domain:

(6)

The observational value for the optimization variable along an avalanche path with fall height Δz is empirically estimated bt (Reference McClung and SchaererMcClung and Schaerer, 2006). For the investigated Wolfsgruben avalanche path the maximum velocity is , which is a reasonable value for avalanches of that size (Reference Gauer, Kern, Kristensen, Lied, Rammer and SchreiberGauer and others, 2007; Reference Fischer, Fromm, Gauer and SovillaFischer and others, 2014). The associated uncertainty may be determined by regression (Reference GauerGauer, 2014) and is estimated at .

Average deposition depth d

The average deposition depth is defined as observed depth, averaged in the affected area. It can be directly measured in the field. The documentation by avalanche experts for the Wolfsgruben avalanche includes deposition depths between 3 and 5 m with a density of . This allows the estimate .

Densification in snow avalanches, comparing released, flowing and deposited snow, can reach a factor of three (Reference AnceyAncey, 2005). Thus, the average deposition depth needs to be defined in terms of the simulation results. We take the peak flow depth h and define averaged in the affected area.

Mass growth G

The mass growth index G is a dimensionless number, describing the increase of flowing avalanche mass due to entrainment. It is defined as the ratio of deposited to released mass

(7)

and has been measured for a variety of avalanche events. However, measurements of the growth index are associated with large variations and uncertainties (Reference Sovilla, Burlando and BarteltSovilla and others, 2006, Reference Sovilla, Margreth and Bartelt2007), for example due to the densification (Reference AnceyAncey, 2005).

Considering densification, released snow mass, affected area and documented deposition depths leads to an estimate of for the Wolfsgruben avalanche.

In terms of the avalanche simulations (cf. Eqns (1) and (4)) the growth index may be written as

(8)

Optimization Concept

The goal of the optimization concept is to provide an objective function as an intuitive, scalar metric, describing the correspondence between simulation and documentation in different categories. The metric allows us to perform a ranking, determining simulation runs and according parameter sets with the highest correspondence to the observation. This matches a selection rule allowing us to accept or withdraw certain parameter combinations, providing input distributions and yielding problem-suited, adjusted output parameter distributions ΩΘ, representing optimal parameter combinations. A flow chart of the optimization concept is depicted in Figure 1. The resulting parameter distributions include model and parameter uncertainties and are of fundamental interest for engineers and scientists.

Using the target function presented here, the optimization of the model parameters could be performed straightforwardly with a Gauss–Newton algorithm, or – more appropriately for the usually coarse DEM grid – a simplified gradient method (Reference SailerSailer and others, 2008). The numerical gradients (Jacobian) obtained in such methods for the optimal parameter set could be used for first-order sensitivity studies (Reference Fellin and OstermannFellin and Ostermann, 2006). However, such inverse calculations are not unique and several local minima could exist. Depending on the initial guess of the model parameters, one of these local minima is found and could be falsely seen as the global minimum, so that a wrong optimal parameter set is chosen (cf., e.g., Reference AnceyAncey, 2005). Investigating the whole physically relevant parameter space and performing statistics on the best simulation runs as proposed here is computationally more expensive, but much more information is produced: local minima can easily be detected and excluded, a complete sensitivity analysis can be performed (Reference FischerFischer, 2013) and reasonable ranges of the model parameters can be determined.

As a measure of the correspondence between observational and simulation optimization variables, we determine αX (Θ) for each optimization variable X = {r, T, F, u max, d, G}, conditional upon the choice of the parameter set Θ = {μ, ξ, eb } summarized in a target function α(Θ). The function describing the correspondence is chosen to be a normalized, Gaussian function with mean and variance :

(9)

A metric is then defined in the interval [0, 1], where αX (Θ) → 0 indicates negligible correspondence and αX (Θ) → 1 optimal correspondence between observation and simulation with respect to the investigated optimization variable. However, the choice of the Gaussian function is arbitrary and it could be replaced by another function (e.g. Heaviside, triangle).

The final correspondence target function α(Θ) is defined as

(10)

With , such that α(Θ) is also bounded by the interval [0, 1]. A ranking of α allows for determining the parameter set Θ with optimal agreement between observation and simulation. The results of the optimization clearly depend on the non-unique definition of the target function, which is a heuristic construction (e.g. not accounting for dependencies between the different optimization variables). The weighting factors wX allow us to emphasize or reduce the impact of certain optimization variables X. For the presented investigation the weighting factors are kept equal, with the exception of excluding single optimization variables. For example, if no information on the variable is available it may not be included in the optimization process, corresponding to a zero weighting factor wX = 0 and a reduced set of optimization variables. The optimization procedure can then be adapted to cases with more or fewer observational data.

Defining a limit α lim matches a selection rule, where simulation runs with a simulation–observation correspondence larger than the limit α(Θ) ≥ α lim are accepted and the other parameter combinations with α(Θ) < α lim are withdrawn. The correspondence limit α lim may be a user-defined value or, in accordance with other engineering applications, be defined in terms of a design event α design (e.g. assuming an acceptable deviation of 5% for each observational variable ). The resulting parameter distributions (including parameter and model uncertainties) then allow us to determine characteristic values for actions, which are based on certain quantiles in engineering concepts like the Eurocodes. Thus, in terms of design events, the limit α lim = α design is calculated following Eqn (10):

(11)

For the simulation runs with α(Θ) ≥ α lim, a frequency analysis of the model parameters Θ is performed. The frequency distribution yielding the adjusted distribution ΩΘ is analyzed for each parameter of the model parameters Θ = {μ, ξ, eb }. The adjusted distribution ΩΘ allows us to further investigate the model behavior. Of particular interest are statistical features, such as the 25%, 50% (median) and 75% quantiles for each parameter Θ25%, Θ50%, Θ75%, minimum and maximum values Θmin, Θmax and the parameter value that correspond to the highest simulation–observation correspondence, i.e. maxα (Θ). Additionally, violine plots (Reference KampstraKampstra, 2008) are a helpful tool for visual interpretation of the adjusted distributions, showing an approximate form of the frequency distribution (Fig. 4).

Fig. 4. Adjusted parameter frequency ΩΘ for the simulation runs with αα lim. The violine plot above summarizes statistical features of the adjusted distributions ΩΘ such as the minimum value Θmin the 25%, 50% and 75% quantiles Θ25%, Θ50% and Θ75% and the maximum value Θmax. To provide a reference to the input distribution the minimum and maximum values are shown.

Analysis of Adjusted Output Parameter Distributions ΩΘ for Θ = {μ, ξ, eb }

The analysis of the output parameter distributions ΩΘ is carried out, varying the set-up with respect to:

optimization variable weighting factors wX , i.e. changing the size of the optimization variable set X = {r, T, F, u max, d, G},

the simulation–observation correspondence limit α lim, the sample size N.

Full set of optimization variables X = {r, T, F, u max, d, G}, N = 104, αlim = αdesign

Here we consider the full set of optimization variables X = {r, T, F, u max, d, G}, with N = 104 sample size of the input parameter distribution at a correspondence limit α lim = α design. Figure 4 shows the adjusted output parameter distributions ΩΘ with . For Ω μ and a clear peak is found, i.e. Θ25%, Θ50%, Θ75% and are relatively close (Table 2). The distribution Ω ξ indicates the tendency of ξ to only provide the desired simulation–observation correspondence lim above a certain boundary value around 7500 < ξ (m s−2). However, no clear peak tendency for ξ is found, which is underlined by the fact that is outside the 75% quantile.

Table 2. Information on the adjusted output parameter distributions ΩΘ, for different sets of optimization variables X. Listed are minimum and maximum values Θmin, Θmax, 25%, 50% and 75% quantiles for each parameter Θ25%, Θ50%, Θ75% and . The same data are visualized in Figure 6

Range and sensitivity of simulated optimization variables X with respect to the adjusted output parameter distributions ΩΘ

The range of simulated optimization variables (X; Table 3) allows us to estimate the quality of the performed back calculation. Simulated runouts are in the range 2190–2263 m with an average value of , which corresponds to the observed range (Table 1). Similar correspondences are obtained for all optimization variables.

Table 3. Distribution of simulated optimization variables X = {r, T, F, u max, d, G} with adjusted output parameter distributions ΩΘ and α lim = α design. denotes the average, and σX the variance

Helpful information to interpret the resulting parameter distributions ΩΘ is the quantification of the flow model parameter Θ sensitivity with respect to the optimization variables X. A Spearman rank correlation analysis is performed (Reference FischerFischer, 2013). The correlation coefficients range from −1 to +1, −1 indicating negative monotonic correlation (increasing parameter leads to decreasing variable), +1 indicating positive monotonic correlation (increasing parameter leads to increasing variable) and 0 indicating no correlation. Table 4 summarizes the parameter-optimization variable correlation coefficients. Besides the obvious relations, which reflect the meaning of the parameters in the employed flow model (e.g. increasing friction leads to decreasing runout or velocities; entrainment rate determines growth index), the quantification allows for a relative evaluation. The information that turbulent and Coulomb friction equally influence the maximum velocity or that no parameter but eb significantly influences the deposition depth is important for model tweaking and the interpretation of the adjusted parameter distributions, especially with respect to a reduced set of optimization variables.

Table 4. Information on the correlation of the adjusted parameter distributions ΩΘ and the optimization variables X. Shown is the Spearman (rank) correlation coefficient r s of optimization variables X and parameters Θ. Only highly significant (p < 0.01) correlations with ∥r s∥ > 0.25 are shown

Reduced set of optimization variables X, N = 104, α lim = α design

With two examples we highlight the benefits of a multivariate parameter optimization. The multivariate parameter optimization is based on a set of different optimization variables X = {r, T, F, u max, d, G}; reducing this set has a significant effect on the adjusted output parameter distributions ΩΘ and thus the ability to quantify optimal parameter set-ups. To investigate the effect, we manipulate the correspondence target function (Eqn (10)), excluding certain optimization variables X by changing the related weighting factors wX = 0.

Reduced set of optimization variables X = {r, T, F, d, G}, i.e.

Excluding the optimization variable u max by setting has a considerable effect on the adjusted output parameter distributions, in particular on the ξ distribution Ω ξ . Figure 5 shows the adjusted parameter distributions ΩΘ with . For μ and eb clear peaks are observed, i.e. the values of μ 25%, μ 50%, , μ 75% and e b,25%, e b,50%, , e b,75% are very close (Table 2). For μ few additional outliers with small values are observed. Investigating the adjusted distribution Ω ξ , only very small constraints for optimal ξ values can be drawn. Besides the effect that no high simulation–observation correspondence solutions are found for small ξ < 1300 m s−2, ξ values are found in the entire investigated range evenly distributed.

Fig. 5. Adjusted parameter frequency ΩΘ for the simulation runs with αα lim: (left) ; (right) WG = Wd = 0. The violine plot above summarizes statistical features of the adjusted distributions ΩΘ such as the minimum value Θmin the 25%, 50% and 75% quantiles Θ25%, Θ50% and Θ75% and the maximum value Θmax. To provide a reference to the input distribution the minimum and maximum values , are shown.

Reduced set of optimization variables X = {r, T, F, umax}, i.e. wG = wd = 0

Excluding the optimization variable average deposition depth d and growth index G has little effect on the adjusted output distribution Ω μ but significantly influences the output ξ and eb distributions (Fig. 5). The number of simulation runs above the correspondence limit α lim is . While clear constraints are found for μ, the distributions Ω ξ and are spread out in almost the entire investigated range (Table 2). Interestingly, ξ values are spread evenly in the investigated interval with a lower bound at ξ min = 5300 m s−2. This means for ξ > ξ min the chances are equal of finding ξ values that lead to simulation runs with simulation–observation correspondence αα lim, which is also reflected by the fact that the parameter value leading to highest correspondence is found at the upper bound of the investigated range . For only a lower bound is found, eb > eb , min = 6250 m2 s−2, corresponding to an exclusion of simulation runs with frontal plowing as entrainment mechanism.

Comparison of different sets of optimization variables X

Excluding single or multiple optimization variables has a significant effect on the information value of adjusted output parameter distributions ΩΘ. The different examples show that form and range of the adjusted parameter distributions significantly varies. Results are summarized in Table 2 and Figure 6 for a visual interpretation of the violine plots. For the presented cases the distribution Ω μ is hardly influenced by the reduced set of optimization variables. The adjusted distribution is unaffected by missing knowledge about the maximum velocity u max but heavily influenced by the optimization variable G. This observation is not surprising, taking into account the correlation of optimization variables and model parameters. Excluding information on maximum velocity u max significantly decreases the information value of the adjusted distribution Ω ξ . A similar effect, but much less pronounced, is observed for leaving out average deposition depth d and growth index G. However, all adjusted output distributions Ω ξ have one thing in common, a lower ξ value bound for high simulation–observation correspondence.

Fig. 6. Comparison of violine plots for the adjusted output parameter frequency ΩΘ for the three cases: (1) all optimization variables (Fig. 6); (2) ; and (3) wG = wd = 0 (Fig. 5). The violine plots summarize statistical features of the adjusted distributions ΩΘ such as the minimum value Θmin the 25%, 50% and 75% quantiles Θ25%, Θ50% and Θ75% and the maximum value Θmax. To provide a reference to the input distribution , the minimum and maximum values are shown. The same data are summarized in Table 2.

Full set of optimization variables X, N = 104, varying αlim

The choice of the correspondence limit α lim determines the number of simulation runs with αα lim. Figure 7 shows the nonlinear decrease of with increasing α lim. Above a certain limit α lim > 0.75, no simulation runs are found, i.e. . For α lim < 0.2, increases dramatically. This can be interpreted as residual correspondence and is due to the fact that some optimization variables are complementary (e.g. an avalanche reaching the affected area implies αT ,d = 0 but αF = 1). For α lim = 0 it naturally is . Figure 8 shows the evolution of the adjusted output distributions Ω μ in dependence on the α lim. Lower α lim lead to larger bounds of ΩΘ, depending on the parameter Θ, i.e. for Ω ξ . In the presented example ΩΘ converges to for large correspondence limits α lim. This is an important observation since multiple distribution maxima may also exist (cf., e.g., Ω μ for α lim = 0.15).

Fig. 7. Number of simulation runs , that are assigned to the adjusted parameter distributions ΩΘ, i.e. with sufficiently large simulation–observation correspondence αα lim.

Fig. 8. Comparison of violine plots for the adjusted output parameter frequency ΩΘ for varying α lim. The violine plots summarize statistical features of the adjusted distributions ΩΘ such as the minimum value Θmin, the 25%, 50% and 75% quantilesΘ25%, Θ50% and Θ75% and the maximum value Θmax. To provide a reference to the input distribution the minimum and maximum values are shown.

Full set of optimization variables X, αlimdesign, variable sample size N, investigating the statistical significance

To determine the statistical stability and draw conclusions on the significance of our results we investigate the dependence of ΩΘ on N. To do so we employ a bootstrapping technique drawing random samples of size N from the sample with 104 simulation runs, increasing the sample size N = 250, … , 104. A first observation is that is directly proportional to sample size N indicating the equal random distribution of the full sample. The minimum sample size for statistically stable results is reached when the parameter distributions stay constant for increasing sample size N (e.g. when Θ50 converges for increasing sample size N).

In the presented case the changes on ΩΘ appear negligible when the sufficient sample size N ≳ 8000 is reached (Fig. 9). Between 3000 < N < 8000 the results are intermediately stable. However, if N is chosen too low, in the presented case N ≳ 2000 the adjusted output distributions are not stable and one may identify statistically non-significant parameter choices, i.e. is too low.

Fig. 9. Comparison of violine plots for the adjusted output parameter frequency ΩΘ for varying sample size N. The violine plots summarize statistical features of the adjusted distributions ΩΘ such as the minimum value Θmin the 25%, 50% and 75% quantiles Θ25%, Θ50% and Θ75% and the maximum value Θmax. Note the scale of the parameter ranges.

One should note that the minimum sample size is dependent on , which is a function of α lim and directly proportional to the sample size N itself. For example, with α lim ≈ 0.4, N ≈ 4000–5000 appears to be sufficient. For smaller N, which corresponds to the situation where the input parameter distribution is not well chosen, the optimal parameter distribution ΩΘ may vary significantly. However, it is beyond the scope of this work to derive a general rule determining the sufficient sample size.

Conclusions

A new optimization concept for snow avalanche simulation has been presented. The method allows us to optimize multiple model parameters using a multivariate evaluation by comparing simulation results with field data based on objective, well-defined optimization variables. With these variables the simulation–observation correspondence is determined, and adjusted output parameter distributions representing optimal parameter combinations are found.

A large number of simulation runs are performed and for each parameter combination Θ = {μ, ξ, eb } a set of simulated optimization variables X = {r, T, F, u max, d, G} is determined. Parameter combinations with high correspondence are identified by analyzing the correspondence between simulated and observed optimization variables using the scalar metric α(X). A selection rule based on a correspondence limit α lim is introduced and linked to a design event correspondence level α design in a consistent manner. Parameter frequency distributions ΩΘ are derived and analyzed statistically (e.g. to determine quantiles or bounds for the model parameters). Violine plots, which allow an intuitive interpretation of the parameter distributions (e.g. to identify multiple distribution maxima), are utilized for visualization. Additionally, an extensive sensitivity analysis allows for linking model parameters and simulation outcome and determining their predictive importance. For the investigated event, the statistical evaluation of the adjusted output parameter distributions showed a clear peak for the Coulomb friction parameter μ 25% = 0.24 < μ < μ 75% = 0.26 and the erosion energy parameter e b,25% = 8825m2 s−2 < eb < e b,75% = 13 925 m2 s−2. Moreover these peaks coincide with the parameter values of maximum correspondence and , indicating high information reliability. For the turbulent friction parameter ξ a lower bound ξ > ξ min = 7500 m s−2, but no clear peak was identified, i.e. the parameter value of maximum correspondence is slightly larger than the 75% quantile value ξ 75% = 13 925ms−2. This means that the optimal value of is some arbitrary value larger than ξ min. However, for ξ values in this range, the magnitude of the turbulent friction diminishes, making its existence questionable, which is in agreement with other observations (Reference Ancey and MeunierAncey and Meunier, 2004). Compared to prior optimization approaches, the optimal values of the turbulent friction appear to be about one order of magnitude larger (0.5–1.5 × 10−3 m s−2 → 0.5–1.5 × 104 m s−2). However, care has to be taken with this comparison, since in other studies entrainment or curvature effects (that introduce additional velocity-dependent friction and thus may have a significant influence on the effective values of the turbulent friction coefficient; cf. Reference Fischer, Kowalski and PudasainiFischer and others, 2012) were partly disregarded. Single parameters cannot be exchanged between different models due to differences in the model implementation. The fact that lower values for the turbulent friction coefficient ξ (i.e. higher friction) appear as high-correspondence solutions when leaving out the maximum velocity u max as optimization variable is in correspondence with prior studies and indicates that lower ξ values lead to avalanche velocities that are too low (as previously noted by Reference GublerGubler, 1987; Reference Fischer, Fromm, Gauer and SovillaFischer and others, 2014; Reference GauerGauer, 2014).

At this point, the advantage of adjusted parameter distributions over singular parameter sets should be highlighted. Statistical information is a major topic in modern engineering design. For example, the upper limit of a 90% confidence interval is usually used as a characteristic value for an action such as impact pressure in engineering codes. The derived parameter distributions include the total model and parameter uncertainties, which are related to a certain deviation of the design event. Thus predictive ensemble simulations are possible and can be used to determine confidence intervals for impact pressures or runouts, i.e. towards a conceptual approach (Reference Meunier and AnceyMeunier and Ancey, 2004) with operational models in 3-D terrain. Furthermore the combination of parameter distributions rather than singular parameter values is indispensable in order to consider the full spectrum of multi-event or multipath analyses.

With different examples we highlight the importance of the multivariate approach. Decreasing the size of the optimization variable set significantly reduces the information value of the adjusted parameter distributions. To adapt the set of optimization variables may be crucial when parameter distributions for different questions or applications have to be determined (e.g., for dam planning, an accurate flow velocity is more important than an accurate runout). The importance of each optimization variable varies for the different model parameters and is directly linked to the results of the sensitivity analysis. For example, the optimization variable of maximum velocity u max is very important to increase the information value of the turbulent friction parameter distribution Ω ξ but only slightly influences the distribution of the Coulomb friction Ω μ . This is also reflected by the fact that the parameter μ is linked to many other optimization variables. The influence of the correspondence limit α lim is studied and the related statistical stability is investigated.

One should note that the presented outcomes may change for other design event definitions, other avalanche paths and may not be transferable to flow models of the same family, or their implementation. Furthermore model parameters may not be constant but may vary in time and space, i.e. for different flow regimes and along the avalanche path. However, the strength of this optimization concept is the possibility that it might be adapted to multi-event or multipath analyses, or to other flow models and their related parameters, and it handles reduced or extended sets of optimization variables, i.e. more or less prior information.

Acknowledgements

This research work was carried out with the financial support of the Federal Ministry of Agriculture, Forestry, Environment and Water Management (BMLFUW), Austrian Service for Torrent and Avalanche Control (WLV); the Austrian Research Fund (FWF, project I 1600-N30 (www. avaflow.org), implemented within the D-A-CH framework under the leadership of the German Research Foundation DFG); and in parts by the US National Science Foundation under grant No. NSF PHY11-25915. The related optimization Python source code is planned to be published in the framework of www.avaflow.org. We thank Marc Adams (Austrian Research Centre for Forests, BFW), the reviewers and the editor for valuable comments which helped to improve the clarity and quality of this paper.

References

Ancey, C (2005) Monte Carlo calibration of avalanches described as Coulomb fluid flows. Philos. Trans. R. Soc. London, 363, 15291550 (doi: 10.1098/rsta.2005.1593)Google Scholar
Ancey, C and Meunier, M (2004) Estimating bulk rheological properties of flowing snow avalanches from field data. J. Geophys. Res., 109, F01004 (doi: 10.1029/2003JF000036)Google Scholar
Ancey, C, Meunier, M and Richard, D (2003) Inverse problem in avalanche dynamics models. Water Resour. Res., 39(4), 1099 (doi: 10.1029/2002WR001749)Google Scholar
Barbolini, M, Gruber, U, Keylock, CJ, Naaim, M and Savi, F (2000) Application of statistical and hydraulic-continuum dense-snow avalanche models to five real European sites. Cold Reg. Sci. Technol., 31(2), 133149 Google Scholar
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. J. Geophys. Res., 117(F1), F01015 (doi: 10.1029/2010JF001957)Google Scholar
Beven, K and Binley, A (1992) The future of distributed models: model calibration and uncertainty prediction. Hydrol. Process., 6(3), 279298 CrossRefGoogle Scholar
Beven, KJ, Smith, PJ and Freer, JE (2008) So just why would a modeller choose to be incoherent? J. Hydrol., 354(1), 1532 CrossRefGoogle Scholar
Bozhinskiy, AN (2008) Modelling of snow avalanche dynamics: influence of model parameters. Ann. Glaciol., 49, 3842 Google Scholar
Bühler, Y, Kumar, S, Veitinger, J, Christen, M and Stoffel, A (2013) Automated identification of potential snow avalanche release areas based on digital elevation models. Natur. Hazards Earth Syst. Sci., 13, 13211335 (doi: 10.5194/nhess-13-1321-2013)Google Scholar
Burkard, A and Salm, B (1992) Die Bestimmung der mittleren Anrissmächtigkeit d0 zur Berechnung von Fliesslawinen. (Technical report) SLF, Davos Google Scholar
Buser, O and Frutiger, H (1980) Observed maximum run-out distance of snow avalanches and the determination of the friction coefficients μ and ξ. J. Glaciol., 26(94), 121130 Google Scholar
Christen, M, Bartelt, P and Kowalski, J (2010a) Back calculation of the In den Arelen avalanche with RAMMS: interpretation of model results. Ann. Glaciol., 51(54), 161168 Google Scholar
Christen, M, Kowalski, J and Bartelt, P (2010b) RAMMS: numerical simulation of dense snow avalanches in three-dimensional terrain. Cold Reg. Sci. Technol., 63, 114 Google Scholar
Dent, J and Lang, T (1980) Modelling of snow flow. J. Glaciol., 26(94), 131140 Google Scholar
Dent, J and Lang, T (1983) A biviscous modified Bingham model of snow avalanche motion. Ann. Glaciol., 4, 4246 Google Scholar
Eckert, N, Parent, E and Richard, D (2007) Revisiting statistical–topographical methods for avalanche predetermination: Bayesian modelling for runout distance predictive distribution. Cold Reg. Sci. Technol., 49(1), 88107 Google Scholar
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. Stoch. Environ. Res. Risk Assess., 22(2), 185206 (doi: 10.1007/s00477-007-0107-4)Google Scholar
Eckert, N, Naaim, M and Parent, E (2010) Long-term avalanche hazard assessment with a Bayesian depth-averaged propagation model. J. Glaciol., 56(198), 563586 Google Scholar
Fellin, W and Ostermann, A (2006) Parameter sensitivity in finite element analysis with constitutive models of the rate type. Int. J. Num. Analyt. Meth. Geomech., 30(2), 91112 Google Scholar
Fischer, J-T (2013) A novel approach to evaluate and compare computational snow avalanche simulation. Natur. Hazards Earth Syst. Sci., 13(6), 16551667 (doi: 10.5194/nhess-13-1655-2013)Google Scholar
Fischer, J-T, Kowalski, J and Pudasaini, S (2012) Topographic curvature effects in applied avalanche modeling. Cold Reg. Sci. Technol., 74–75, 2130 (doi: 10.1016/j.coldregions. 2012.01.005)Google Scholar
Fischer, J-T, Fromm, R, Gauer, P and Sovilla, B (2014) Evaluation of probabilistic snow avalanche simulation ensembles with Doppler radar observations. Cold Reg. Sci. Technol., 97, 151158 (doi: 10.1016/j.coldregions.2013.09.011)Google Scholar
Gauer, P (2014) Comparison of avalanche front velocity measurements and implications for avalanche models. Cold Reg. Sci. Technol., 97, 132150 (doi: 10.1016/j.coldregions.2013.09.010)Google Scholar
Gauer, P and Issler, D (2004) Possible erosion mechanisms in snow avalanches. Ann. Glaciol., 38, 384392 Google Scholar
Gauer, P, Kern, M, Kristensen, K, Lied, K, Rammer, L and Schreiber, H (2007) On pulsed Doppler radar measurements of avalanches and their implication to avalanche dynamics. Cold Reg. Sci. Technol., 50(1), 5571 (doi: 10.1016/j.coldregions.2007.03.009)Google Scholar
Gauer, P, Medina-Cetina, Z, Lied, K and Kristensen, K (2009) Optimization and probabilistic calibration of avalanche block models. Cold Reg. Sci. Technol., 59(2–3), 251258 (doi: 10.1016/j.coldregions.2009.02.002)Google Scholar
Gray, JMNT, Wieland, M and Hutter, K (1999) Gravity-driven free surface flow of granular avalanches over complex basal topography. Proc. R. Soc. London, 455(1985), 18411874 CrossRefGoogle Scholar
Gruber, U and Bartelt, P (2007) Snow avalanche hazard modelling of large areas using shallow water numerical methods and GIS. Environ. Model. Softw., 22(10), 14721481 (doi: 10.1016/j.envsoft.2007.01.001)Google Scholar
Gubler, H (1987) Measurements and modelling of snow avalanche speeds. IAHS Publ. 162 (Symposium at Davos 1986 – Avalanche Formation, Movement and Effects), 405420 Google Scholar
Harbitz, C (1998) A survey of computational models for snow avalanche motion. (NGI Report 581220-1, SAME Collaboration 32) Norwegian Geotechnical Institute, Oslo Google Scholar
Iman, R and Conover, W (1982) A distribution-free approach to inducing rank correlation among input variables. Commun. Stat.: Simul. Comput., 11(3), 311334 Google Scholar
Issler, D and Gauer, P (2008) Exploring the significance of the fluidized flow regime for avalanche hazard mapping. Ann. Glaciol., 49, 193198 Google Scholar
Issler, D and 11 others (2005) A comparison of avalanche models with data from dry-snow avalanches at Ryggfonn, Norway. Proceedings of the 11th International Conference and Field Trip on Landslides, Norway, 1–10 September, 2005. A.A. Balkema, Rotterdam; Taylor & Francis Group, London, 173179 Google Scholar
Iverson, R M (2012) Mechanics of debris flows and rock avalanches. (Handbook of Environmental Fluid Dynamics 43) CRC Press, Hobart, 573–587Google Scholar
Johannesson, T, Gauer, P, Issler, P and Lied, K (2009) The design of avalanche protection dams: recent practical and theoretical developments. Directorate General for Research, European Commission,Brussels Google Scholar
Kampstra, P (2008) Beanplot: a boxplot alternative for visual comparison of distributions. J. Stat. Softw., Code Snippets, 28(1), 19 Google Scholar
Maggioni, M and Gruber, U (2003) The influence of topographic parameters on avalanche release dimension and frequency. Cold Reg. Sci. Technol., 37(3), 407419 (doi: 10.1016/S0165-232X(03)00080-6)Google Scholar
McClung, DM and Schaerer, P (2006) The avalanche handbook, 3rd edn. The Mountaineers Books, Seattle, WA Google Scholar
McMillan, H and Clark, M (2009) Rainfall–runoff model calibration using informal likelihood measures within a Markov chain Monte Carlo sampling scheme. Water Resour. Res., 45(4) (doi: 10.1029/2008WR007288)Google Scholar
Mergili, M, Schratz, K, Ostermann, A and Fellin, W (2012) Physically-based modelling of granular flows with Open Source GIS. Natur. Hazards Earth Syst. Sci., 12, 187200 (doi: 10.5194/ nhess-12-187-2012)Google Scholar
Mergili, M, Marchesini, I, Rossi, M, Guzzetti, F and Fellin, W (2013) Spatially distributed three-dimensional slope stability modelling in a raster GIS. Geomorphology.Google Scholar
Meunier, M and Ancey, C (2004) Towards a conceptual approach to predetermining long-return-period avalanche run-out distances. J. Glaciol., 50(169), 268278 Google Scholar
Monaghan, J (1992) Smoothed particle hydrodynamics. Ann. Rev. Astron. Astrophys., 30, 543574 Google Scholar
Montanari, A, Shoemaker, CA and Van de Giesen, N (2009) Introduction to special section on uncertainty assessment in surface and subsurface hydrology: an overview of issues and challenges. Water Resour. Res., 45(12), W00B00 Google Scholar
Naaim, M, Furdada, G and Martinez, H (2002) Calibration and application of the MN2D dynamics model to the avalanches of Las Leñas (Argentina). Natur. Hazards Earth Syst. Sci., 2(3/4), 221226 Google Scholar
Naaim, M, Faug, T, Naaim, F and Eckert, N (2010) Return period calculation and passive structure design at the Taconnaz avalanche path, France. Ann. Glaciol., 51(54), 8997 Google Scholar
Naaim, M, Durand, Y, Eckert, N and Chambon, G (2013) Dense avalanche friction coefficients: influence of physical properties of snow. J. Glaciol., 59(216), 771782 (doi: 10.3189/2013JoG12J205)Google Scholar
Oberguggenberger, M, King, J and Schmelzer, B (2009) Classical and imprecise probability methods for sensitivity analysis in engineering: a case study. Int. J. Approx. Reason., 50(4), 680693 (doi: 10.1016/j.ijar.2008.09.004)Google Scholar
Pitman, E, Nichita, C, Patra, A, Bauer, A, Sheridan, M and Bursik, M (2003) Computing granular avalanches and landslides. Phys. Fluids, 15, 3638 (doi: 1070-6631/2003/15(12)/3638/9)Google Scholar
Platzer, K, Bartelt, P and Kern, M (2007) Measurements of dense snow avalanche basal shear to normal stress ratios (S/N). Geophys. Res. Lett., 34(7), L07501 (doi: 10.1029/2006GL028670)Google Scholar
Pudasaini, SP (2012) A general two-phase debris flow model. J. Geophys. Res., 117(F3), F03010 (doi: 0148-0227/12/2011JF002186)Google Scholar
Pudasaini, SP and Hutter, K (2003) Rapid shear flows of dry granular masses down curved and twisted channels. J. Fluid Mech., 495(1), 193208 (doi: 10.1017/S0022112003006141)Google Scholar
Sailer, R, Rammer, L and Sampl, P (2002) Recalculation of an artificially released avalanche with SAMOS and validation with measurements from a pulsed Doppler radar. Natur. Hazards Earth Syst. Sci., 2, 211216 Google Scholar
Sailer, R and 6 others (2008) Snow avalanche mass-balance calculation and simulation-model verification. Ann. Glaciol., 48, 183192 Google Scholar
Salm, B, Burkhard, A and Gubler, HU (1990) Berechnung von Fliesslawinen: eine Anleitung fuer Praktiker; mit Beispielen. Eidg. Inst. Schnee-Lawinenforsch. Mitt., 47, 137 Google Scholar
Sampl, P and Zwinger, T (2004) Avalanche simulation with SAMOS. Ann. Glaciol., 38, 393398 Google Scholar
Savage, SB and Hutter, K (1989) The motion of a finite mass of granular material down a rough incline. J. Fluid Mech., 199(1), 177215 Google Scholar
Sovilla, B, Burlando, P and Bartelt, P (2006) Field experiments and numerical modeling of mass entrainment in snow avalanches. J. Geophys. Res., 111(F3), F03007 (doi: 10.1029/2005JF000391)Google Scholar
Sovilla, B, Margreth, S and Bartelt, P (2007) On snow entrainment in avalanche dynamics calculations. Cold Reg. Sci. Technol., 47(1–2), 6979 (doi: 10.1016/j.coldregions.2006.08.012)Google Scholar
Sovilla, B, Schaer, 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. J. Geophys. Res.: Earth Surf., 113(F1) (doi: 10.1029/2006JF000688)Google Scholar
Stein, M (1987) Large sample properties of simulations using Latin hypercube sampling. Technometrics, 29(2), 143151 Google Scholar
Teich, M, Fischer, J-T, Feistl, T, Bebi, P, Christen, M and Grêt-Regamey, A (2014) Computational snow avalanche simulation in forested terrain. Natur. Hazards Earth Syst. Sci., 14(8), 22332248 (doi: 10.5194/nhess-14-2233-2014)Google Scholar
Tiefenbacher, F and Kern, M (2004) Experimental devices to determine snow avalanche basal friction and velocity profiles. Cold Reg. Sci. Technol., 38(1), 1730 (doi: 10.1016/S0165-232X(03)00060-0)Google Scholar
Vera Valero, C, Wikstroem Jones, K, Bühler, Y and Bartel, P (2015) Release temperature, snow-cover entrainment and the thermal flow regime of snow avalanches. J. Glaciol., 61(225), 173184 (doi: 10.3189/2015JoG14J117)Google Scholar
Voellmy, A (1955) Über die Zerstörungskraft von Lawinen. Schweiz. Bauztg., 73, 159165, 212–217, 246–249, 280–285Google Scholar
Vrugt, JA, Ter Braak, CJ, Gupta, HV and Robinson, BA (2008) Equifinality of formal (dream) and informal (glue) Bayesian approaches in hydrologic modeling? Stoch. Environ. Res. Risk Assess., 23(7), 10111026 (doi: 10.1007/s00477-008-0274-y)Google Scholar
Zwinger, T, Kluwick, A and Sampl, P (2003) Numerical simulation of dry-snow avalanche flow over natural terrain. In Hutter, K and Kirchner, N eds Dynamic response of granular and porous materials under large and catastrophic deformations. Springer Verlag, Berlin, 161194 Google Scholar
Figure 0

Fig. 1. The proposed framework consists of three parts: observation, simulation and optimization. Monte Carlo simulation runs are performed for a large number of parameter combinations Θ = μ, ξ, eb following the input parameter distributions . The simulation software, SamosAT, is utilized, taking into account boundary and initial conditions, which are obtained from the observations. The set of optimization variables X = {…} (with variables runout r, matched affected area (true) T, exceeded affected area (false) F, maximum velocity umax, average deposition depth d and mass growth G) is defined in terms of the simulation results (flow depths h, velocities u and mass m of the avalanche) as well as observations and their related uncertainty . The simulation–observation correspondence α(Θ) is quantified for each simulation run. Based on a consistently defined correspondence limit αlim, parameter combinations are withdrawn or accepted, yielding the final results: adjusted parameter distributions ΩΘ.

Figure 1

Fig. 2. Wolfsgruben avalanche (images WLV). (a) Wolfsgruben avalanche path (vertical drop Δz = 984 m): central flowline in black, release area (196 225 m2, with mean slope angle 36.5°) in blue, affected area (64 153 m2, at 14.5°) in orange. (b) Destroyed house in the runout area. Back-calculated pressures of damages on infrastructure range from about 7.5 to 17.5 kPa.

Figure 2

Table 1. Observational variables for the Wolfsgruben avalanche

Figure 3

Fig. 3. Sketched simulation results (e.g. simulation outline plim; blue) and affected area (orange), superimposed with avalanche path domain and coordinate system along the central flowline (black). Runout r, , matched area Amatched (green) and exceeded area Aexceeded (red) and their complements (blue and gray area), which lead to the optimization variables for true and false prediction T, F, maximum velocity umax, average deposition depth d and mass growth G.

Figure 4

Fig. 4. Adjusted parameter frequency ΩΘ for the simulation runs with ααlim. The violine plot above summarizes statistical features of the adjusted distributions ΩΘ such as the minimum value Θmin the 25%, 50% and 75% quantiles Θ25%, Θ50% and Θ75% and the maximum value Θmax. To provide a reference to the input distribution the minimum and maximum values are shown.

Figure 5

Table 2. Information on the adjusted output parameter distributions ΩΘ, for different sets of optimization variables X. Listed are minimum and maximum values Θmin, Θmax, 25%, 50% and 75% quantiles for each parameter Θ25%, Θ50%, Θ75% and . The same data are visualized in Figure 6

Figure 6

Table 3. Distribution of simulated optimization variables X = {r, T, F, umax, d, G} with adjusted output parameter distributions ΩΘ and αlim = αdesign. denotes the average, and σX the variance

Figure 7

Table 4. Information on the correlation of the adjusted parameter distributions ΩΘ and the optimization variables X. Shown is the Spearman (rank) correlation coefficient rs of optimization variables X and parameters Θ. Only highly significant (p < 0.01) correlations with ∥rs∥ > 0.25 are shown

Figure 8

Fig. 5. Adjusted parameter frequency ΩΘ for the simulation runs with ααlim: (left) ; (right) WG = Wd = 0. The violine plot above summarizes statistical features of the adjusted distributions ΩΘ such as the minimum value Θmin the 25%, 50% and 75% quantiles Θ25%, Θ50% and Θ75% and the maximum value Θmax. To provide a reference to the input distribution the minimum and maximum values , are shown.

Figure 9

Fig. 6. Comparison of violine plots for the adjusted output parameter frequency ΩΘ for the three cases: (1) all optimization variables (Fig. 6); (2) ; and (3) wG = wd = 0 (Fig. 5). The violine plots summarize statistical features of the adjusted distributions ΩΘ such as the minimum value Θmin the 25%, 50% and 75% quantiles Θ25%, Θ50% and Θ75% and the maximum value Θmax. To provide a reference to the input distribution , the minimum and maximum values are shown. The same data are summarized in Table 2.

Figure 10

Fig. 7. Number of simulation runs , that are assigned to the adjusted parameter distributions ΩΘ, i.e. with sufficiently large simulation–observation correspondence ααlim.

Figure 11

Fig. 8. Comparison of violine plots for the adjusted output parameter frequency ΩΘ for varying αlim. The violine plots summarize statistical features of the adjusted distributions ΩΘ such as the minimum value Θmin, the 25%, 50% and 75% quantilesΘ25%, Θ50% and Θ75% and the maximum value Θmax. To provide a reference to the input distribution the minimum and maximum values are shown.

Figure 12

Fig. 9. Comparison of violine plots for the adjusted output parameter frequency ΩΘ for varying sample size N. The violine plots summarize statistical features of the adjusted distributions ΩΘ such as the minimum value Θmin the 25%, 50% and 75% quantiles Θ25%, Θ50% and Θ75% and the maximum value Θmax. Note the scale of the parameter ranges.