A novel approach to evaluate and compare computational snow avalanche simulation

. An innovative approach for the analysis and interpretation of snow avalanche simulation in three dimensional terrain is presented. Snow avalanche simulation software is used as a supporting tool in hazard mapping. When performing a high number of simulation runs the user is confronted with a considerable amount of simulation re-sults. The objective of this work is to establish an objective, model independent framework to evaluate and compare re-sults of different simulation approaches with respect to indicators of practical relevance, providing an answer to the important questions: how far and how destructive does an avalanche move down slope. For this purpose the Automated Indicator based Model Evaluation and Comparison (AIMEC) method is introduced. It operates on a coordinate system which follows a given avalanche path. A multitude of simulation runs is performed with the snow avalanche simulation software SamosAT (Snow Avalanche MOdelling and Simulation – Advanced Technology). The variability of pressure-based run out and avalanche destructiveness along the path is investigated for multiple simulation runs, varying release volume and model parameters. With this, results of deterministic simulation software are processed and analysed by means of statistical methods. Uncertainties originating from varying input conditions, model parameters or the different model implementations are assessed. The results show that AIMEC contributes to the interpretation of avalanche simulations with a broad applicability in model evaluation, comparison as well as examination of scenario variations.


Introduction
Over the past years simulation software for geophysical mass flows such as snow avalanches has become an important tool in hazard potential estimations (Barbolini et al., 2000;Gruber and Margreth, 2001;Pudasaini and Hutter, 2007).These tools are used to provide an estimate of the avalanche's destructive power and spatial extent.
In European engineering practice block models based on the work of Voellmy (1955) were mainly used (Salm et al., 1990) until the end of the 20th century.These were followed by the first approaches of numerical software (Bartelt et al., 1999;Christen et al., 2002), which operated on a two dimensional simplification of the mountain topography.Similar to statistical models (Lied and Bakkehøi, 1980;McClung and Lied, 1987) the model results include the avalanche run out on a predefined mountain profile.Since the beginning of the 21st century, flow models, mainly based on the two dimensional depth-averaged shallow water equations and the continuum mechanical model for granular flows (Savage andHutter, 1989, 1991), have been implemented in simulation software (Pitman et al., 2003;Patra et al., 2005;Sampl and Granig, 2009;Mancarella and Hungr, 2010;Christen et al., 2010;Mergili et al., 2012).
This modern software usually provides a user interface and functionality for geographic data handling.It is used to handle data input, output and provides an intuitive way of data visualization in natural three dimensional terrain.With the development of new models and software tools, amount and scope of in-and output information has increased from point information up to results in two or three spatial dimensions over time.
A detailed model analysis such as sensitivity analysis, calibration, evaluation with field measurements or the comparison to other tools require a technically comprehensive, standardized way of result processing.Previous work has been done for statistical models or models and software operating in two dimensional terrain.For example, Gauer et al. (2009b) performed a probabilistic calibration of avalanche block models and Gauer et al. (2010) back-calculated over 300 Norwegian and Austrian avalanches with a traditional numerical block model investigating the parameter ranges.Ancey (2005) performed a block model Monte Carlo calibration on seven two dimensional path profiles.Barbolini and Savi (2001) estimated uncertainties in avalanche hazard mapping using a one dimensional avalanche dynamics model.Barbolini et al. (2000) compared two statistical and three hydraulic-continuum models proposing a scheme for avalanche hazard zoning that integrates the statistical and dynamic models.
In the case of modern simulation software the available evaluation methods are limited.Simulation results have mainly been evaluated in predefined cross sections along the avalanche path (Pirulli and Sorbino, 2008;Christen et al., 2010;Buehler et al., 2011;Mergili et al., 2012;Bartelt et al., 2012).Besides missing evaluation approaches, the computational effort restricted the evaluation of modern simulations software to a low number of simulation runs for back calculations of single events (Sailer et al., 2002;Issler et al., 2005).A model independent and objective analysis procedure which is capable of processing large sets of simulation results in three dimensional terrain is in great demand.
Therefore the goal of this work is to evaluate high numbers of simulation runs under the conditions of varying input.Avalanche simulation software results describe the temporal evolution of flow variables, flow depth and slope parallel velocities in two to three space dimensions.To answer the main questions of main importance: how far and how destructive does an avalanche move down slope, the simulation results have to be interpreted by an avalanche engineer, e.g. the position of run out, by determining the point where avalanche velocity or flow depth fall below a certain threshold.This method works for single simulations but is impracticable for scenario evaluations.For this reason the Automated Indica- In the first part of this paper a detailed description of the new AIMEC approach is given.In the second part the AIMEC analysis of two simulation scenarios, variation of release volume and variation of model parameters, are presented.The results are discussed and an outlook on possible further applications is presented.

Methods
In the present study the proprietary snow avalanche simulation software SamosAT (Snow Avalanche MOdelling and Simulation -Advanced Technology, Zwinger et al., 2003;Sampl and Zwinger, 2004;Granig and Oberndorfer, 2008;Sampl and Granig, 2009) is used.An avalanche scenario for the avalanche path Ryggfonn, Norway provides the input data.SamosAT has been extended, such that a multitude of simulation runs are automatically performed.The simulation results are analysed in a standardized way using the AIMEC approach.A coordinate transformation in an avalanche path dependent system provides the basis to define the indicators, a new metric describing the main avalanche features.AIMEC has been developed in the MATLAB environment.The following sections provide a complete description of the in-and output as well as the theoretical basis and its implementation.

Simulation software -SamosAT
The software consists of two basic models, a dense flow avalanche (DFA) model and a powder snow avalanche (PSA) model in order to describe the descent of dry snow avalanches (Sampl, 1999;Zwinger et al., 2003;Sampl and Zwinger, 2004;Granig and Oberndorfer, 2008;Sampl and Granig, 2009).In SamosAT avalanches are assumed to consist of a three layer structure with different flow densities, see Fig. 1.The DFA is reflected by the dense flow layer and modelled as a shallow flow in two dimensions above the mountain surface.The PSA model includes all three avalanche layers where the powder snow layer is computed as three dimensional flow.The transition and powder snow layer are assumed to develop above the dense flow layer over time.This includes air layers a few hundred metres above ground.The transition layer is not spatially resolved in the model and mainly determines the flux of snow mass from dense to powder snow layer or vice versa (Sampl and Granig, 2009).
In the DFA model the released volume is discretized in a large number of mass elements and the depth-averaged model equations are solved with a smoothed particle hydrodynamics (SPH) scheme (Monaghan, 1992).To compute a PSA the three dimensional air flow is computed with a SIMPLE (semi-implicid pressure linked equations) scheme (Patankar, 1980;AVL, 2009), coupled to the ice particle equations of motion which are solved in a Lagrangian framework.In this work the SamosAT DFA model is used.However, the new approach can also be applied to the SamosAT PSA model or any other models of the same type.

Simulation input
To perform an avalanche simulation the mountain topography represented by a digital elevation model (equally spaced, discrete raster of 5 m ×5 m), an assigned area of potential avalanche release, the release depth and the flow model parameters have to be provided.The dense flow snow avalanche model of SamosAT includes seven model parameters specifying the bottom friction τ b relation: with with the specifications of snow density ρ, Coulomb friction µ, coefficients of the turbulent friction term κ, R, B, minimum shear stress τ 0 and the phenomenological fluidization factor R 0 s .K accounts for the terrain curvature in flow direction, h for the flow depth and g z is the slope normal component of the gravitational acceleration (Gray et al., 1999;Pudasaini and Hutter, 2003;Pudasaini et al., 2005a,b;Fischer et al., 2012).
SamosAT has been calibrated for extreme avalanches with a 150 yr return period (Granig and Oberndorfer, 2008).A set of reference parameters has been determined in order to best reproduce avalanche run outs based on Austrian avalanche data.

Simulation output
The output of SamosAT DFA and other avalanche simulation models are flow depth and slope parallel velocities (h(x, y, t), v(x, y, t)) at a constant density (ρ).x, y denote the two dimensional Cartesian coordinates.Other field quantities such as flow pressure, momentum, front velocities, etc. have to be deduced from the simulation output.For engineering questions the most important simulation results are the maximum values over the simulation duration, i.e. the peak velocities.Here peak corresponds to the observed maximum value over the computation time at a given location, v peak (x, y) = max t {v(x, y, t)}.With the peak velocities and the snow density the peak pressure field is defined, corresponding to the maximum impact pressures over time.Peak pressures are an important instrument for the simulation result interpretation.Figure 2a shows a three dimensional view of the test site with the peak pressure results.This formulation of impact pressure assumes a classical drag coefficient of 2, corresponding to a flat plate perpendicular to the flow, which is often used in avalanche simulation.Recent experiments have shown that measured impact pressures can significantly exceed the classically expected ones depending on flow type and shape of the structure (Baroudi et al., 2011).

Avalanche scenario
The presented approach is conceptually designed such that it is applicable to any avalanche path.To demonstrate its potential an avalanche scenario at the avalanche test site Ryggfonn, Norway, is used.This full-scale avalanche test site has been in operation since 1980 (Gauer et al., 2007(Gauer et al., , 2008(Gauer et al., , 2009a)).
The test site has a vertical drop of about 900 m and a horizontal length of approximately 2000 m.The inclination of the main track is about 30 • on average and 17 • in the run out area.In the run out zone, a 75 m wide (crown length) and 16 m high catching dam is situated.Its slope angle is 40 • .Beyond the dam, the path crosses a creek and ends in a counter slope.The path is slightly channelized.At Ryggfonn, observed avalanche sizes range between class 2 (mass of 10 5 kg) and class 4 (mass of 10 7 kg), and may even reach class 5, according to the Canadian snow avalanche size classification (Gauer et al., 2008).
The main goal of this work is to investigate the dynamic response of the simulation software to variable simulation scenarios.For this reason input conditions and model parameters are varied and the corresponding simulation results are compared.As basis for any optimization or comparison an initial guess or reference has to be provided.Here a reference scenario is defined by choosing representative model parameters and boundary conditions.In this study the reference parameter set is based on the calibration guidelines of SamosAT (Granig and Oberndorfer, 2008, ρ = 200 ).The corresponding release area and depth are chosen according to field observations at the Ryggfonn test site and classical mapping strategies in hazard mapping (Maggioni and Gruber, 2003).A release area of ≈ 55 × 10 3 m 2 with a reference release depth of 1 m are assigned, compare Fig. 2.This corresponds to a total release mass of ≈ 11 × 10 6 kg (with 200 kg m −3 snow density), which is in the upper range of observed Ryggfonn avalanches.

Simulation scenarios
For the purpose of this analysis two simulation scenarios are performed, each with J simulation runs: The simulation software has been adapted to perform a multitude of simulation runs automatically.For each simulation run predefined results are exported as a uniform, discrete, two-dimensional raster.With this approach it is guaranteed that all simulation runs are performed in the same way.Furthermore, time and data space is saved compared to single simulations performed manually.
The size of the simulation scenario is determined to provide a sufficient representation of the possible result spectrum.Preanalysis showed that for a Monte Carlo sample size of 1000, the result variations are small comparing different samples.Increasing the sample size did not lead to significantly different values.A relatively high number of simulation runs (1000) were computed in an acceptable time range, considering 2-5 min per simulation run.

Coordinate transformation
Avalanche characteristics such as run out are directly related to information about the main flow direction, which defines the avalanche path.Simulation results do not directly include this information.For this reason quantitative statements about important characteristics derived from the two dimensional results cannot be made without additional information.A coordinate system following the avalanche path automatically provides this information.Therefore the coordinate transformation is a crucial step of AIMEC.It allows to define indicators, which represent metrics in an avalanche path dependent framework.Coordinate systems aligned with the avalanche path, so called curvilinear coordinate systems have already been used in the derivation of model equations (Pudasaini and Hutter, 2003;Pudasaini et al., 2005b).Here similar coordinate systems are used for the analysis of simulation results.
The simulation results are given in a global x, y Cartesian coordinate system.To investigate the results (here the peak pressure P (x, y)) with respect to a specific avalanche path, the coordinate transformation is performed.For simplicity the ˜are dropped in the following.In Fig. 2b the Cartesian coordinate systems superimposed with path domain dependent coordinate system are shown for the Ryggfonn avalanche path.The horizontally projected coordinate in the main flow direction s -the longitudinal coordinate, and in the cross flow direction l -the lateral coordinate, is introduced.
Technically the coordinate systems are represented by discrete two dimensional rasters.The original raster is aligned with the Cartesian x, y coordinates.The new irregular raster is in alignment with the main flow direction coordinate s and the cross flow coordinate l, given by the central flow line.The implementation of the coordinate transformation is described in the following sections and includes four steps: -choice of a central flow line z(x, y), -creation of the avalanche path domain of width w, -creation of the new discrete two dimensional coordinates (s, l), -result transformation (x, y) → (s, l).
If the avalanche path flow direction s is directly aligned with the Cartesian coordinate x, x = s, y = l no coordinate transformation is necessary and indicators can be directly defined (Sect.2.7).

Central flow line and path domain
The new coordinate system is aligned with the central flow line and embedded in the path domain.In the present case the choice of the central flow line (z(x, y)) is based on field observation and historical data and fulfills the following requirements: -alignment with the main flow direction of the observed avalanche, -longitudinal extent large enough to cover the expected flow extent, z(x, y) defines the l = 0 line of the new coordinate system s, l.It consists of N straight segments, defined by a finite number of N + 1 points S i in the x − y plane: Depending on the research question, different ways to choose the central flow line are possible and encounter different challenges.For example, a central flow line oriented along the talweg (following the maximum slope gradient, starting at the release area) may not always be aligned with main flow direction in natural terrain.For avalanche paths splitting into two channels one main channel or an averaged flow direction has to be defined.For avalanche paths where no prior knowledge about the main flow direction is available an automated central flow line detection from simulation results would be desirable.However, in the presented case the central flow line is defined according to corresponding field observations.
The path domain of width w is created along the central flow line.To do so, domain boundaries O l,r i are defined for each domain segment along the points S i of the central flow line.For the N + 1 points S i the angle of deviation is computed.It is and With this, the left and right

Coordinate system and result transformation
In order to transform the simulation results from the global x, y Cartesian coordinate system into the path dependent s, l system, bounded by the path domain, a discrete allocation is performed.The original x, y Cartesian coordinate system is discretized on a regular n x ×m y raster with a fixed cell size of 5 m ×5 m.The avalanche path dependent coordinate system s, l is represented by an irregular raster of size n s × m l .The aim of this transformation is to prevent any data loss, therefore the spatial resolution of the new raster is at least equal to or higher than the resolution of the original raster.To ensure this the following steps are performed.
The number of elements m l i on each lateral path segment line O l i O r i is determined using the Bresenham line algorithm (Bresenham, 1965).This algorithm originates from the field of computer graphics and is used to determine the number of discrete elements by producing a close approximation of straight lines on a regularly spaced two dimensional discrete raster.The new raster size in lateral direction is obtained as a global maximum of all i lateral segment lines, In the longitudinal direction the number of elements n s,l i along the left O l i O l i+1 , and the number of elements n s,r i of the right O r i O r i+1 domain segment lines parallel to the central flow line are determined.The longitudinal number of elements for each segment is defined as local maximum of the left and right n s i = max{n s,l i , n s,r i } number of raster elements.Finally each domain segment is divided into n s i − 1 sub segments and the total longitudinal raster size, is computed, compare conceptual sketch in Fig. 3b.The content of the corresponding cell of the original raster is assigned to each cell of the newly obtained raster.To transfer the cell-centered values, a nearest neighbour method is sufficient.In Fig. 3b the discrete representation of the new coordinate systems is shown.Although an overlapping of cells might occur in heavily curved avalanche paths, no simulation results are lost, in particular there are no holes in the new raster.The new variable cell size has to be taken into account when performing the coordinate transformation.The avalanche simulation results are now available on a new raster aligned with the avalanche path.

Indicators, path dependent metrics
The main goal of defining indicators is to introduce metrics which describe the main characteristics of the simulation results.Generally indicators have to be consistent and comparable.They can be based on any simulation result, e.g.flow depths, velocities or any deduced quantities such as pressures.Two pressure-based indicators representing information about the spatial extent and destructiveness along the path of the avalanche are presented below.The choice of these indicators is motivated by considerations regarding the application and behaviour of the underlying flow model.
Impact or peak pressure results are of major interest in snow avalanche simulation.They are estimates for the avalanche's destructive potential and allow multiple applications in simulation result interpretation, engineering practice or hazard mapping.In most countries criteria for avalanche hazard mapping and run out are mainly based on pressure estimates (Johannesson et al., 2009).In modern risk assessment or cost-benefit analysis, detailed information on the vulnerability of various elements at risk is required (Eckert et al., 2012).Vulnerability in turn is directly related to impact pressures, which are an indicator for the magnitude of a process.
When compared to field measurements, a pressure-based run out length can be defined as a demarcation line of destruction and thus be estimated by back calculations of damages.Past avalanche run outs are mainly remembered on account of their destructive impact, which can be interpreted as a certain peak pressure level.
Due to the importance and applicability of peak pressure results, they provide the basis to introduce scalar indicators, that represent metrics of the avalanche dynamics along the path, as well as the run out behaviour of one specific simulation run.To define these indicators, the two dimensional peak pressure distribution P (s, l) is reduced to a distribution along the avalanche path coordinate s.Average Pcross (s) and maximum P max cross (s) peak pressures are defined for every cross section along the central flow line, In Fig. 4 the cross-sectional average and maximum peak pressures are displayed along the avalanche path.The central flow line profile z(s, 0) provides a topographic reference for the result interpretation.The peak pressure originates from the peak velocities and considering Eq. (3) all statements regarding peak pressures also apply to peak velocities.

Spatial indicator -run out
In the following, run out is referred to as a pressure-based run out limit.It is defined via a pressure limit of the crosssectional average peak pressure Pcross (s).In order to detect the location of final avalanche deceleration, the user-defined pressure threshold P limit should be in the range of impact pressures expected in the run out zone.To examine the relative response or sensitivity of the simulation results to input variations, consistent definitions of indicators are essential.The exact definition of the pressure-based run out may be adapted to the research question.Different pressure limits or other run out definitions, e.g. based on a pressure limit of the maximum cross-sectional pressure P max cross (s) or other simulation results would lead to different absolute run out values.Intuitively pressure limits may be chosen according to guidelines of hazard mapping, which vary for different zones and countries, see Johannesson et al. (2009) for additional information.Here P limit = 1 kPa is chosen, which corresponds to the lowest pressure limit according to the Austrian hazard  14)) and average Pcross (s) (red, compare Eq. ( 13)) peak pressure along the avalanche path coordinate s.Additionally the averaged maximum peak pressure AMPP (Eq.( 15)) and the pressure-based run out is shown.The profile of the central flow line z(s, 0) is shown for orientation on the right ordinate.The pressure values t belong the reference simulation run with a release volume of 55 × 10 3 m 3 .mapping guidelines.Following the longitudinal path coordinate s the first point where the peak pressure limit is exceeded ( Pcross (s) > P limit ) is referred to as starting point s = s start .The last point, where the cross-sectional average peak pressure Pcross (s) falls below the limit ( Pcross (s) < P limit ) is respectively referred to as run out s = s runout , compare Fig. 4.
With this definition of run out the final position of the avalanche head is detected in the cross section of the central flow line, as long as the result extent is in rough alignment to and covered by the path domain.A pressure threshold dependent run out definition (P limit > 0 kPa) is reasonable.For many flow models in avalanche simulation practice no realistic stopping behaviour may be modelled in flat natural terrain.Therefore run out criteria based on flow depth or zero pressure might not yield satisfying results within a finite simulation time.Moreover different model approaches and their implementation might result in diffusive run out behaviour.Therefore additional stopping criteria are often defined in the practical application.By using a run out criterion based on a sensible pressure threshold, these problems are avoided and no additional criteria are necessary.

Destructiveness indicator -rAMPP
The relative averaged maximum peak pressure (rAMPP) is introduced as a relative measure to evaluate the avalanche's destructiveness along the avalanche path comparing different simulation runs.The absolute value AMPP is defined as the average of the cross-sectional maximum peak pressure along the avalanche path, from initiation s start to run out For the comparison of j = 1, . . ., J different simulation runs with the relative value: with respect to AMPP 0 , the result corresponding to the reference simulation (compare Sect. 2.4), are used.The rAMPP value of the reference simulation itself is rAMPP 0 = 1, compare Figs. 5 and 6b.With the rAMPP indicator the range of result variability in the dynamic behaviour (destructiveness) of the avalanche along the path is investigated.It is a global metric independent of a certain position along the path and contains valuable information for the interpretation of avalanche simulations.Especially in combination with the pressure-based run out indicator, the rAMMP highlights the influence of input variations on simulation results.

AIMEC -analysis
The aim of the AIMEC analysis is a clear and well-arranged representation of result variations for multiple simulation runs.For the present analysis the input of the analysis consists of: -peak pressure results of various simulation runs with SamosAT, -the central flow line, see Fig. 2, -the avalanche path domain width, w = 500 m, -the pressure limit P limit = 1 kPa.

Simulation scenario 1: variation of release volume
For this analysis, J = 100, simulation runs with continuously increasing release volume (2.75 − 275 × 10 3 m 3 , constant re-lease area with varying release depth) are performed, the corresponding run out and rAMMP indicator results are shown in Fig. 5. Generally an increase in release volume leads to an increase in both indicators, run out and rAMPP.A statistical measure for this monotonic dependency is Spearman's rank correlation coefficient r X Y .It is defined as the covariance divided by the product of the standard deviations of the ranks of two variables X, Y .The correlation coefficients are computed for the variables' release volume, run out and rAMPP and tested non-directional for significance against the hypothesis of no correlation (p values < 0.05).The analysis shows a significant correlation for all investigated variables.The correlation coefficient of run out and rAMPP, r runout rAMPP = 0.995, indicates that high run outs are accompanied by a high avalanche destructiveness.The correlation coefficients of release volume, run out and rAMPP indicate a slightly higher correlation for volume and rAMPP, r V rel rAMPP = 0.999, than volume and run out, r V rel rAMPP = 0.995.One characteristic of the Ryggfonn avalanche path is the catching dam in the run out area (at path position s = 1670 m).For small release volumes the influence of the dam is clearly observable.The rAMMP is constantly increasing, while the run out stays rather unchanged.At a volume of around 27.5 × 10 3 m 3 the avalanche flows over or around the dam.For higher volumes a continuous increase in run out and rAMMP is observed, although the rate of change decreases with increasing volume.The range of the rAMMP indicator is ≈ 0.14−2.57and of the run out indicator around ≈ −94 to 137 m projected run length s.For approximately one-fourth of the reference release volume (13.75 × 10 3 m 3 ) the run out indicator is 1658 m and the rAMPP value is 0.38, indicating average maximum velocity values of about 40 % less compared to the reference simulation run.Table 1 summarizes the results.The general and important observation is that the impact on the indicators decreases from low to high volumes in the investigated range.This indicates, that the relative volume variation is determining the rate of indicator change.
One important observation is that while the rAMMP continuously increases, the run out variations partly contradict the intuition that an increase in volume correlates with an increase in run out, compare run out variations in the range of s = 1670 − 1800 m with release volumes 24.75 − 137.5 × 10 3 m 3 .This variation also contributes to a slightly lower correlation coefficient of volume and run out.Reasons for the observed run out variations could be related to the validity of a depth-averaged model in significantly curved terrain.In the dam area the shallowness assumption of the underlying model could be violated (Domnik and Pudasaini, 2012).Further reasons for this counter intuitive result variations are believed to be related to the numerical model implementation and need to be further investigated.
However, the example considered here shows that the ability to detect such unexpected variations originating from flow

Simulation scenario 2: variation of model parameters
A simulation scenario varying the seven model parameters is performed with a Monte Carlo simulation including J = 1000 simulation runs.For each parameter (snow density ρ, Coulomb friction µ, turbulent friction coefficients κ, R, B, minimum shear stress τ 0 and fluidization factor R 0 s ) 1000 random values equally distributed over an interval based on the calibration standard values of SamosAT (Granig and Oberndorfer, 2008) ±50 % were assumed.This approach allows rather extreme combinations of parameter set-ups.Parameter specifications such as negative friction are not allowed.
One goal of this analysis is to estimate the scope of possible simulation results.Figure 6 shows the density distributions of J = 1000 indicator results (rAMMP vs run out) of two Monte Carlo simulations with different release volumes (release volume 13.75 × 10 3 m 3 and 55 × 10 3 m 3 ).Shown is the result density, which increases from blue to red, the quantiles are marked at the colour bar.To allow a direct comparison of the destructiveness indicator rAMMP the reference simulation run is kept the same in each scenario.Although this approach might seem counter-intuitive in the case of the small release volume it simplifies the comparison of the performed simulation scenarios, compare Figs. 5 and 6, pink (+) and gray (X) mark the simulation run with large and small release volume and reference parameter set.With this the rAMMP results of the parameter variation are consistent and comparable to the release depth variation.
Considering the magnitude of the input parameter variation a high result variability is observed.In the case of the smaller release volume V rel ≈ 13.75 × 10 3 m 3 the run out variability reaches from 500 to 1764 m including run outs in the upper part of the avalanche path, compare the large one, respectively.It is important to note that in each case the rAMMP is relative to the reference simulation result with release volume V rel ≈ 55 × 10 3 m 3 .Table 2 summarizes minimum, maximum, median and average value with according standard variation of the marginal run out and rAMMP distribution.Comparing the decrease of the run out standard deviations from the small (σ = 246 m) to the big release volume (σ = 56 m) the importance of the initial conditions is obvious.Again Spearman's rank correlation coefficient r X Y is calculated as statistical measure for all investigated variables X and Y .Table 2 lists the correlation coefficients with high significance (p values < 0.05) for all investigated indicators and parameters.Besides the parameter R 0 S all correlations are at a significant level.Correlation of run out and rAMMP is slightly smaller for the large release volume r rAMPP runout = 0.796 than for the small one r rAMPP runout = 0.828.The run out indicator is mainly influenced by bottom friction (coulomb friction µ and minimum shear stress τ 0 ) and turbulent friction (turbulent friction parameter κ), while destructiveness shows highest correlation to the turbulent friction The influence of the minimum shear stress τ 0 decreases for an increase in volume.
The results reveal that different parameter sets may lead to the same run out indicator but with large variations in the rAMMP indicator.This means that the simulated avalanche reaches the same run out result but the avalanches destructive power along the path differs significantly.Especially for models or simulation tools optimized or calibrated only with regard to run out, this information is essential in order to perform a meaningful result interpretation.
A high run out density is observed in the dam area.This is a characteristic of the Ryggfonn avalanche path and points out the topographic dependence of the run out in SamosAT.So-called starving avalanches (Bartelt et al., 2007) that stop in the upper part of the avalanche path, which are also observed in nature are only possible to simulate for the lower release volume.Generally, the position of the run out is connected to terrain parts with flat slopes, which is in correspondence to field observations (Sovilla et al., 2010).

Conclusions
In this work, the Automated Indicator based Model Evaluation and Comparison (AIMEC), an innovative approach to analyse two dimensional results of snow avalanche simulation software was introduced.SamosAT was used to compute avalanche scenarios at the Norwegian avalanche test site Ryggfonn with variations of release volume and model parameters.The simulation results were analysed with respect to a spatial run out indicator and the rAMMP destructiveness indicator.These peak pressure-based indicators were defined in an avalanche path dependent coordinate system with respect to the questions of main interest: how far and how destructive does an avalanche move down slope?
With two scenarios, the simulation software result variability and correlations to the input are quantified.A wide range of results are obtained, which is also observed in nature.Both release volume and model parameter variations showed a significant influence on the investigated indicators.The simulation scenario with varying release volume revealed unexpected and counter intuitive result variations.They are considered to be connected to the model's numerical implementation or a violation of the model's basic assumption which could be solved with a full dimensional description in the dam area.
Over all results the rAMMP destructiveness indicator varies between 5 and 280 % compared to the reference simulation run.The possible run out indicator ranges from 500 to 1918 m (30-116 %) measured along the flow direction coordinate.The variability of the rAMMP destructiveness indicator is similar in all investigated scenarios.The variability of the run out indicator is strongly dependent on the release volume and shows the highest variability in the case of parameter variation for the small release volume.Highest correlations are observed for run out and destructiveness together with release volume.The model parameter variation revealed high correlations for destructiveness and the turbulent friction coefficient, independently of release volume.The run out is mainly correlated to coulomb friction and the turbulent friction coefficient for the large release volume and in the case of small release volume to the minimum shear stress.
As shown, the presented approach can be successfully applied to investigate model behaviour on a single path by running multiple scenarios.Moreover it provides a broad applicability and new opportunities in avalanche simulation.The presented modelling and analysis approach serve as basis to evaluate simulation results with field data and measurements such as velocity data (Fischer et al., 2013).Additionally, AIMEC allows for a direct and quantitative comparison between different models and avalanche paths.This is of great interest for a wide spectrum of users and avalanche simulation applications.In engineering practice it offers the possibility to easily and objectively compare simulation results.It is a new approach to evaluate large simulation scenarios assessing the result variability arising from input uncertainties.This can be helpful especially in complex terrain, since additional result variations may occur for different software due to their individual numerical implementation.Researchers can apply the AIMEC approach to investigate the influence of changes in model parameters, and an insight into the uncertainties of the simulation results is given independently of the software.It was proven to be highly valuable to produce input data needed for parameter optimization of simulation software.For example, Teich et al. (2013) apply the run out indicator to evaluate and implement forests' retarding influence on run out distances of 40 small-to medium-scale avalanches observed in forested terrain.
The presented analysis approach can be further adapted to different problems introducing new indicators (e.g.growth index for entrainment and deposition or deposition morphology).This includes an integrated approach comparing different types of gravitational flows (landslides, debris, mud flows, rock fall, etc.) using different type of models (e.g.Pudasaini, 2012).It provides a comprehensive basis for a conceptual approach to running and evaluating deterministic simulation software in a probabilistic framework.

Fig. 1 .
Fig. 1.Avalanche layer structure assumed in the SamosAT model.The dense core of the snow avalanche flows on the sliding surface.Snow particles are exchanged in the resuspension transition layer to enter the powder snow layer.

Fig. 2 .
Fig. 2. Avalanche simulation at Ryggfonn, Norway.(a) Three dimensional view of the topography (100 m main contour lines) with the peak pressures field P (x, y).The central flow line is marked in black, the release area in pink.(b) Avalanche path dependent domain of width w = 500 m and coordinate system along the central flow line.The simulation results given in the global Cartesian coordinate system (x, y) are transformed in an avalanche path dependent (s, l) coordinate system, with the longitudinal coordinate s in direction of the main avalanche flow and the lateral coordinate l.Shown is the peak pressure field of an avalanche simulation run at Ryggfonn.Blue and red shadings mark the different path domain segments, the black line marks the central flow line of the avalanche dependent coordinate system.
domain boundaries at distance w/2 are determined for each point S i .This corresponds to segment boundary lines perpendicular to the central flow line at the first O l 0 O r 0 and last segment O l N +1 O r N +1 .Any other segment boundary line O l i O r i , i = 1, . . ., N corresponds to the bisecting line of the respective angle between the two central flow line segments.This leads to the fact that l is not necessarily locally perpendicular to s all over the path.A conceptual sketch of the central flow line and the corresponding path domain is shown in Fig. 3a.The central flow line z(x, y) and the path domain width w should be determined carefully.Simulation results outside of the domain, overlapping raster segments in the outer region of the domain due to sharp bends, too short path segments in the central flow line, or a too large path width w should be avoided.But, simulation results close to the path domain borders are generally low and their influence on the indicators results are negligible.These issues could also be investigated by implementing an adaptive domain width w.

Fig. 3 .
Fig. 3. Sketched simulation results (e.g. 1 kPa; green) superimposed with avalanche path domain and new coordinate system along the central flow line z(x, y).(a) Points of the central flow line S i and domain of width w with domain boundaries O l i and O r i and corresponding deviation angle θ i are shown.(b) Original raster (gray) with coordinates x, y, superimposed by irregular, path dependent raster (yellow) with coordinates s, l.Blue and red shadings mark the different path domain segments.

Fig. 4 .
Fig.4.Cross-sectional maximum P max cross (s) (blue, compare Eq. (14)) and average Pcross (s) (red, compare Eq. (13)) peak pressure along the avalanche path coordinate s.Additionally the averaged maximum peak pressure AMPP (Eq.(15)) and the pressure-based run out is shown.The profile of the central flow line z(s, 0) is shown for orientation on the right ordinate.The pressure values t belong the reference simulation run with a release volume of 55 × 10 3 m 3 .

Fig. 5 .
Fig. 5. Results of the AIMEC analysis, J = 100 simulation runs with continuously increasing release volume V rel .Shown is rAMPP vs run out.Reference point is the simulation run with V rel ≈ 55 × 10 3 m 3 (pink +) release volume.Additionally the result of the simulation run with V rel ≈ 13.75 × 10 3 m 3 (gray X) is marked.

Table 1 .
Runout and rAMPP for reference V rel = 55.5×10 3 m 3 and V rel = 13.75 × 10 3 m 3 simulation runs.Additionally the minimum, maximum, median, average and standard variation of the result indicators are shown.V rel = 13.75 × 10 3 m 3V rel = 55.5 × 10 3 implementation is one of the strengths of the presented approach.

Fig. 6 .
Fig. 6.AIMEC results of the J = 1000 Monte Carlo simulation runs with release volumes V rel ≈ 13.75 × 10 3 m 3 (a) and 55 × 10 3 m 3 (b).Shown is the result density, which increases from blue to red, the quantiles are marked at the colour bar.The reference simulation with standard parameters is marked in each case (gray X and pink + ).
parameter κ.A negative correlation indicates that an increase in friction parameters is accompanied with an increase and friction and thus a decrease in run out and destructiveness.