r.avaﬂow v1, an advanced open-source computational framework for the propagation and interaction of two-phase mass ﬂows

. r.avaﬂow represents an innovative open-source computational tool for routing rapid mass ﬂows, avalanches, or process chains from a deﬁned release area down an arbitrary topography to a deposition area. In contrast to most existing computational tools, r.avaﬂow (i) employs a two-phase, interacting solid and ﬂuid mixture model (Pudasaini, 2012); (ii) is suitable for modelling more or less complex process chains and interactions; (iii) explicitly considers both entrainment and stopping with deposition, i


Introduction
Rapid flows or avalanches of snow, debris, rock, or ice, or processes, process chains, or process interactions involving more than one type of movement or material, frequently lead to loss of life, property, and infrastructures in mountainous areas worldwide.All state-of-the-art methods for anticipating the occurrence, characteristics, and dynamics of such events rely on computer simulations.On the one hand, models attempt to identify those areas where mass flows are likely to release (landslide susceptibility; Guzzetti, 2006;Van Westen et al., 2006).On the other hand, they attempt to anticipate the motion of rapid mass flows once they are released (Hungr et al., 2005a).Whilst conceptual models (Lied and Bakkehøi, 1980;Gamma, 2000;Wichmann and Becht, 2003;Horton et al., 2013;Mergili et al., 2015) are employed to identify possible impact areas at broad scales, physically based dynamic models are used for the detailed back-analysis or prediction of specific events.
Advanced fluid dynamics offer a broad array of physically based dynamic modelling approaches for mass flows, mostly referred to as granular avalanches or debris flows.Such models often centre on two-dimensional "shallow flow" equations, but they vary considerably among themselves in terms of their concept, complexity, and capacity to model specific types of phenomena.Voellmy (1955) pioneered mass flow modelling, followed by the work of Grigoriyan et al. (1967), Savage and Hutter (1989), Takahashi (1991), Iverson (1997), Pitman and Le (2005), and many others (see Pudasaini and Hutter, 2007 for a review).Savage and Hutter (1989) introduced depth-averaged mass and momentum conservation equations which were later utilized, modified, and extended by Mangeney et al. (2003Mangeney et al. ( , 2005)), Denlinger and Iverson (2004), and McDougall and Hungr (2004, 2005).The Savage and Hutter (1989) model was further extended to include the effects of pore fluid by Iverson and Denlinger (2001), Savage and Iverson (2003), Pitman and Le (2005), Pudasaini et al. (2005), Pastor et al. (2009), and Hutter and Schneider (2010a, b).Still, these approaches either represent effectively one-phase models, or do not fully consider the twophase nature of most mass flows.More recently, the software GeoClaw and its extension D-Claw consider shallow water and quasi-two-phase flows (M.J. Berger et al., 2011;Iverson and George, 2016).Pudasaini (2012) introduced a general two-phase mass flow model including several essentially new physical aspects of two-phase solid-fluid mixture flows.In comparison to one-phase models, amongst a few other twophase approaches (e.g.Kowalski and McElwaine, 2013), this appears suitable for the realistic simulation of most types of process chains and interactions such as overtopping of a lake and a subsequent flood or debris flow due to the impact of a landslide into the lake.
Various types of numerical schemes have been used to solve mass flow model equations in order to redistribute mass and momentum (e.g.Davis, 1988;Toro, 1992;Nessyahu and Tadmor, 1990;Tai et al., 2002;Wang et al., 2004).Previously, equations were commonly formulated and solved for predefined types of topographies (Pudasaini et al., 2005(Pudasaini et al., , 2008;;Wang et al., 2004), whereas a mathematically consistent application to arbitrary mountain topographies -and therefore to real-world conditions -still remains a challenge (Mergili et al., 2012).This issue is closely related to the fact that the model equations are commonly expressed in topography-following coordinates hardly compatible with global Cartesian coordinates, which usually appear in geographic information systems (GIS) and are referred to as GIS coordinates in the following.Nevertheless, some of the mass flow models mentioned have been implemented in computational tools used for hazard mapping and zoning, such as DAN (Hungr, 1995), TITAN2D (Pitman et al., 2003b;Pitman and Le, 2005), SamosAT (Sampl and Zwinger, 2004), or RAMMS (Christen et al., 2010a, b).Hergarten and Robl (2015) developed a modelling tool relying on the opensource flow solver GERRIS (Popinet, 2009).
None of these models explicitly consider stopping and deposition, and they offer only basic functionalities for simulating chains or interactions of two-phase mass flows.There is, however, a particular need to appropriately consider process chains and interactions in mass flow simulations: some of the most destructive events in history have evolved from cascading effects, such as the 1970 Huascarán event in Peru (Evans et al., 2009) or the 2002 Kolka-Karmadon event in Russia (Huggel et al., 2005).
The present work addresses some of the needs and issues raised by introducing the multifunctional open-source computational framework r.avaflow, employing an enhanced version of the Pudasaini (2012) two-phase flow model for routing mass flows from a defined release area down arbitrary topography to a deposition area.Next, we introduce the structure and components of r.avaflow (Sect.2).Then, we perform two computational experiments in order to demonstrate the functionalities of the computational framework (Sect.3).We discuss the implementation of r.avaflow and the implications of our findings (Sect.4) and finally conclude with the key messages of the work and a brief outlook on the next steps (Sect.5).
2 The computational framework r.avaflow 2.1 Computational implementation r.avaflow computes the propagation of mass flows from one or more given release areas over a defined basal topography until (i) all the material has stopped and deposited; (ii) all the material has left the area of interest; or (iii) a user-defined maximum simulation time has been reached.r.avaflow is developed along two lines with regard to its software environment and operation, r.avaflow [EXPERT] and r.avaflow [PROFESSIONAL].The latter represents a stand-alone version with still reduced functionalities.It is operated through a graphical user interface (GUI), suitable for practitioners.The present work, however, refers to r.avaflow [EXPERT] which is implemented as a raster module of the open-source software package GRASS GIS 7 (Neteler and Mitasova, 2007;GRASS Development Team, 2016).We use the Python programming language for data management, preprocessing, and post-processing tasks (module r.avaflow).The flow propagation procedure (see Sect. 2.3 and 2.4) is written in the C programming language (sub-module r.avaflow.main).Together with Python, the R software environment for statistical computing and graphics (R Core Team, 2016) is employed for built-in validation and visualization functions.Figure 1 illustrates the logical framework of r.avaflow.
Multiple model runs may be executed in parallel, exploiting all computational cores available (see Sect. 2.5).This speeds up the processing considerably and allows the use of r.avaflow on computational clusters.Parallelization is implemented at the Python level (Mergili et al., 2014(Mergili et al., , 2015)): for each model run, a batch file is produced within the module r.avaflow.This batch file calls the Python-based sub-module r.avaflow.mult,launching r.avaflow.main,which is then executed with the specific parameters for the associated model run.Thereby, the Python library "Threading", a higher-level threading interface is exploited.The Python class "Queue" is employed for handling the queue of items to be processed.
r.avaflow was developed and tested with the operating systems (OS) Ubuntu 12.04 and 16.04 LTS, and Scientific Linux 6.6 (Red Hat).It is expected to work on other UNIX systems, too.A simple user interface is available.However, the tool may be started more efficiently through command line parameters, enabling a straightforward batching on the shell script level.This feature facilitates model testing and the combination with other GRASS GIS modules.
Experiments where parallel processing is not applied are performed on an Intel ® Core i7 975 with 3.33 GHz and 16 GB RAM (DDR3, PC3-1333 MHz), exploring a maximum of eight cores through hyperthreading and using the OS Ubuntu 12.04 LTS.All experiments with parallel processing are performed on the Vienna Scientific Cluster, serving with approximately 2020 nodes (Supermicro X9DRD-iF Board), each equipped with an Intel Xeon processor E5-2650v2 with 2.6 GHz and 8 × 8 GB RAM.The OS for these computations is Scientific Linux 6.6 (Red Hat).

Input and output
The key input parameters of r.avaflow are summarized in Table 1.Essentially, r.avaflow relies on (i) a digital terrain model (DTM) representing the elevation of the basal surface (in the release areas beneath the release mass) before the event under investigation, (ii) raster maps of the spatial distribution of the solid and fluid release heights or hydrographs of solid and fluid release, and (iii) a set of flow parameters (Table 2).Input raster maps of the entrainable solid and fluid heights, and a raster map or value defining the empirical entrainment coefficient (needed for entrainment) are optional.Instead of the solid and fluid release and entrainable heights, the total heights and fixed values of the solid concentration may be defined.
There is no restriction imposed on the arrangement of the release cells.With the term "cell", we refer to a regular, equidistant, square, ground-projected computational/numerical unit, i.e. an element of a GIS raster.Patches of cells where the release height is larger than zero may be defined in various parts of the investigation area.An arbitrary number of release hydrographs -each associated with a given set of coordinates -can be defined alternatively or in addition to the release masses.This allows the simulation of complex interactions between different types of processes (see Sect. 3).Hydrographs are defined through their solid and fluid heights at the centre point of the hydrograph profiles, and by the solid and fluid flow velocities.The flow height distribution along the hydrograph profile -which should be aligned perpendicular to the main flow direction -is derived from the assumptions of a horizontal cross section of the flow table and a maximum profile length (Fig. 2).
Mandatory parameters further include the time interval at which output maps are written t out (s), the maximum time after which the simulation terminates, and the threshold flow height for visualization and validation H t (m; see Table 1).Optional parameters further include raster maps of the observed impact area and deposition height, as well as a set of flow path coordinates (for validation and visualization; see Fig. 1 and Sect.2.6).An exhaustive list of input parameters is provided in the user manual of r.avaflow, available at http://www.avaflow.org/software.html.
If a single model run is executed (see Fig. 1), the output of r.avaflow consists in raster maps of solid, fluid, and total flow heights, flow velocities in x and y direction and in absolute terms, pressures and kinetic energies, and the change of the basal topography (only relevant with entrainment or stopping; see Sect.2.4).All raster maps are produced for each output time step (defined by t out ) and for the maximum over all time steps.Further, a table summarizing the   maximum solid and fluid flow heights and velocities as well as flow volumes and kinetic energies for all output time steps is produced.Optionally, solid and fluid output hydrographs are generated for an arbitrary number of given output hydrograph profiles (see Table 1 and Fig. 2).With multiple model runs, the results of each single run are aggregated to impact or deposition indicator indices (see Sect. 2.5).In the present work, we focus on the output heights, hydrographs, and indices when analysing the results, rather than on velocities or deduced results such as pressures or kinetic energies (see Sect. 3).

Mass and momentum evolution
The core functionality of r.avaflow consists in the redistribution of mass and momentum, employing a dynamic flow model and a numerical scheme.Thereby, the tool offers implementations (i) of a single-phase shallow water model with Voellmy friction relation (Christen et al., 2010a, b;Fischer et al., 2012) and (ii) essentially the Pudasaini (2012) twophase flow model with ambient drag (Kattel et al., 2016) and a set of additional numerical treatments (complementary functions) outlined in Sect.2.4.In the present work, we only consider the implementation (ii).It builds on the conservation of mass and momentum, computed separately but simultaneously for the solid and fluid components of the flow.
A system of six differential equations (expressed in locally topography-following coordinates) represents the basis for a set of six flux and source terms, regarding solid and fluid flow depths (D s , D f ), solid momentum M sx and fluid momentum The Pudasaini (2012) model employs the Mohr-Coulomb plasticity for the solid stress.The fluid stress is modelled as a solid-volume, fraction-gradient-enhanced, non-Newtonian viscous stress.The generalized interfacial momentum transfer includes viscous drag, buoyancy, and virtual mass induced by relative acceleration between the phases.A new generalized drag force is proposed that covers both solidlike and fluid-like contributions.Strong coupling between the solid-momentum and the fluid-momentum transfer leads to simultaneous deformation, mixing, and separation of the phases.Inclusion of the non-Newtonian viscous stresses is important in several aspects.The advection and diffusion of the solid volume fraction play an important role.The model includes a number of innovative, fundamentally new, and dominant physical aspects.Please consult Pudasaini (2012) for the full details of the model, including the corresponding equations.The flow parameters required are summarized in Table 2.
Solving the differential equations and propagating the flow from one cell to the next requires the implementation of a numerical scheme.For this purpose, r.avaflow employs a highresolution total variation diminishing non-oscillatory central differencing (TVD-NOC) scheme, a numerical scheme used to avoid unphysical numerical oscillations (Nessyahu and Tadmor, 1990).Cell averages of all six state variables are computed using a staggered grid: the system is moved half of the cell size with every time step; the values at the corners of the cells and in the middle of the cells are computed alternatively at half and full time steps, respectively.The TVD-NOC scheme with the minmod limiter has successfully been applied to a large number of mass flow problems (Tai et al., 2002;Wang et al., 2004;Mergili et al., 2012;Pudasaini and Krautblatter, 2014;Kafle et al., 2016;Kattel et al., 2016).
The input and output of r.avaflow (see Sect. 2.2) is discretized on the basis of GIS coordinates, i.e. in cells which are rectangular in shape in the ground projection.For the numerical solution, the cell lengths in x and y directions, and the area, are corrected for the local slope in order to maintain consistency with the state variables expressed in the local topography-following coordinates.Gravitational acceleration in the topography-following x, y, and z directions -representing a fundamental input to the Pudasaini (2012) model equations -is computed from the DTM, employing a finite central differencing scheme.All input heights H (m) are expressed in a vertical direction and are converted into depths D (m) expressed in a direction normal to the local topography as in the model equation formulation.The resulting depths are converted into heights for output.The time step length t is dynamically updated according to the CFL condition (Courant et al., 1967;Tai et al., 2002;Wang et al., 2004).
We note that all total (solid plus fluid) heights and depths represent the real-world heights and depths only if all the www.geosci-model-dev.net/10/553/2017/Geosci.Model Dev., 10, 553-569, 2017 pores in the solid material are filled with fluid (pores filled with air are excluded).

Complementary functions
Table 3 summarizes some additional functions of r.avaflow.
The functions with ID 1-3 have been introduced to compensate for deficiencies of the numerical scheme and its implementation experienced with complex real-world flows (see Sect. 4).Entrainment and stopping, in contrast, represent dynamic functions not covered by the Pudasaini (2012) model and are executed at the end of each time step (see Fig. 1).
Even though the separation of the complementary functions from the TVD-NOC scheme, and their treatment in a simple forward Euler manner, can be questioned physically and mathematically, we consider the current implementation a reasonable first approximation (see Sect. 4).We now elaborate the concepts employed for entrainment as well as for stopping and deposition in more detail.Full handling of the evolution of the basal topography within the TVD-NOC scheme is not straightforward and could also produce some diffusion.Therefore, as entrainment is not included in the original Pudasaini (2012) model, entrainment is treated as a complementary function in a first step.We note, however, that the time steps at which entrainment and the change of the basal topography are updated are identical to the time steps of the numerical scheme.The potential solid and fluid entrainment rates q E,s and q E,f (expressed perpendicular to the basal topography) build on the user-defined empirical entrainment coefficient C E (see Table 2) and the solid and fluid momenta.We assume a vertically homogeneous solid fraction α s,Emax within the entrainable material, which is reflected in the ratio between q E,s and q E,f : (1) The fact that the basal velocities, which are relevant for entrainment, are lower than the depth-averaged velocities is not explicitly considered, but has to be reflected in the value of C E .q E,s and q E,f are always positive.The solid and fluid changes of the basal topography, H E,s and H E,f , due to entrainment are where H E,s (t − t) and H E,f (t − t) are the change of the basal topography at the start of the time step, H Emax,s and H Emax,f are the maximum entrainable depths at the given cell, t is the time passed at the end of the time step, t is the time step length, and β is the local slope of the basal surface.The division by cos β approximates the conversion from depths to heights.The solid and fluid entrained depths cos β are added to the solid and fluid flow depths.We further assume that entrainment increases the solid and fluid momentum of the flow in each direction by the product of the entrained solid and fluid depth and the velocity in the given direction (M E ; Fig. 3a).The basal topography and, consequently, the x and y cell sizes, cell areas, and gravitational acceleration components in x, y, and z direction are updated after each time step.
The changes in gravitational acceleration also influence the magnitude of the frictional terms (Pudasaini and Hutter, 2003), which are important for stopping processes.In the literature, few approaches explicitly consider stopping processes directly in their numerical scheme by operator splitting methods coupled with the determination of admissible stresses (e.g.Mangeney et al., 2003;Zhai et al., 2015).Here, in order to consider stopping which occurs at a spatial scale that is not numerically resolved, we choose a different approach by proposing the dimensionless factor of mobility (FoM), relating the distance required for stopping s stop to the numerical spatial resolution s in the direction of move- ment.The flow stops if s stop ≤ s, i.e.FoM ≤ 1 (see Fig. 3b): To estimate s stop , we formulate the energy balance considering that the initial kinetic energy at an initial velocity v 0 and the change of potential energy while travelling the distance s stop have transformed in dissipative energy due to Coulomb friction, which dominates close to stopping.With this, the energy balance estimate yields Consequently, where δ is the basal friction angle, β v is the slope angle in the direction of movement, and g is gravitational acceleration (see Table 2).According to Eq. ( 6) the stopping distance s stop is positive for δ > β v , meaning that stopping is possible when the friction angle is higher than the slope angle, i.e. in particular at flat or even counter slopes.We note that, by a simple transformation of Eq. ( 6), FoM can alternatively be derived by relating the stopping time to the time step length.The stopping criterion is only relevant for v 0 > 0.
FoM can relate to various spatial units, such as (i) a single cell; i.e.FoM is computed separately for each cell (it may happen that stopping of the flow occurs at a certain cell, but not at its neighbour cells); (ii) v 0 and β v are averaged over a certain cell neighbourhood to compute FoM, so that stopping occurs at patches of adjacent cells; and (iii) β v and the associated component of v are averaged over the entire flow area.This means that the entire flow stops at once.
The third possibility is currently implemented with r.avaflow as an optional function.If activated, the simulation terminates as soon as stopping occurs and the entire flow material is deposited.Note that, in the current implementation, stopping and deposition always consider the total mass, without differentiating between the solid and the fluid components.This simplification is reasonable for flows characterized by a relatively small fluid volume fraction.The change of basal topography due to entrainment H E after the last time step is subtracted from the height of the deposited material H D in order to derive the change of basal topography (or net deposition) H C at the end of the simulation (positive for an increase, negative for a decrease of terrain elevation).

Multiple model runs
r.avaflow includes a built-in function to perform multiple model runs at a time with controlled or random variation of uncertain input parameters between given lower and upper thresholds.Essentially, this concerns the flow parameters (see Table 2) but also the solid concentration of the release mass α s0 .Multiple parameters can be varied at a time.This procedure serves two purposes: -It facilitates multi-parameter sensitivity analysis and optimization efforts.
-The results of all model runs are aggregated to an impact indicator index (III) and a deposition indicator index (DII), each in the range 0-1.III represents the fraction of model runs where H Max ≥ H t at a given cell, whilst DII represents the fraction of model runs where H C ≥ H t at a given cell.III and DII are used to directly account for uncertain input parameters in the simulation result.
The model runs can be assigned to multiple computational cores (parallel processing), enabling the exploitation of highperformance computational environments (see Sect. 2.1).

Validation and visualization
r.avaflow can be used to produce map layouts and animations of the key results (see Fig. 1).It further includes builtin functions to validate the model results against observations.Validation relies (i) on the availability of a raster map of the observed impact or deposition area of the event under investigation, (ii) on a user-defined profile along the main flow path (see Table 1), or (iii) on measurements of H or v Table 4. Validation criteria used in r.avaflow (see also Fig. 4).S: single model run, binary simulation result; M: multiple model runs, simulation result in the range 0-1.The concepts of CSI and D2PC are taken from Formetta et al. (2016).All validation parameters are computed for H Max (OIA as reference) and/or H C (ODA as reference), depending on which of the reference data are available.at selected coordinates and time steps.Those cells with observed impact or deposition are referred to as observed positives (OPs), those without observed impact or deposition as observed negatives (ONs).When using the observed impact area (OIA) as reference, all cells with H Max ≥ H t are considered as predicted positives (PPs), all cells with H Max < H t are considered as predicted negatives (PNs).When using the observed deposition area (ODA) as reference, all cells with H C ≥ H t are considered as PPs, all cells with H C < H t are considered as PNs.Intersecting ONs and OPs with PPs and PNs results in four validation scores: true positive (TP), true negative (TN), false positive (FP), and false negative (FN) predictions (Fig. 4).TN strongly depends on the size of the area of interest.It is normalized to 5 • (TP + FN) − FP in order to allow a meaningful comparison of model performance among different case studies.These scores build the basis for most of the validation parameters described in Table 4.Only the excess travel distance L relies on the observed and simulated terminal points of the flow, based on a user-defined longitudinal profile.We note that this profile is only needed for validation but is not used for the mass flow simulation itself.

Scope
Values of L > 0 and FoC > 1 indicate conservative results (simulated impact or deposition area is larger than observed impact or deposition area) whilst values of L < 0 and FoC < 1 indicate non-conservative results.CSI, D2PC, and AUROC do not allow to conclude on the conservativeness of the results.L, FoC, CSI, and D2PC as defined in Table 4 target at the validation of H Max or H C derived with one single model run.With multiple model runs (see Sect. 2.5), those validation parameters are computed separately for each run, allowing to conclude on the sensitivity of the model performance to given input parameters, or to optimize input parameter values.In this sense, optimum parameters always refer to one particular criterion, and different criteria may suggest different optimum parameter values.
In contrast, ROC (receiver operating characteristic) curves are used to test the performance of the overall output of multiple model runs.Such curves are produced for III (OIA as reference) and/or DII (ODA as reference): the true positive rate is plotted against the false positive rate for various levels of III or DII.The area under the curve connecting the resulting points, AUROC, is used as an indicator for model performance (AUROC ≈ 1 indicates an excellent performance; see Fig. 4 and Table 4).
Further, the difference between observed and simulated values of H and v at selected sets of coordinates and points of time can be analysed.This function is mainly useful for very well-documented case studies, such as laboratory experiments, and is not further used in the present work.

Topographic setup
In a first step, the potential of r.avaflow for simulating process chains is demonstrated, considering the interaction between one or more landslides, a reservoir, and the dam impounding the reservoir.This experiment represents a followup to the work of Pudasaini (2014), Kafle et al. (2016), and Kattel et al. (2016).We construct a generic landscape of size 3200 m × 2000 m, illustrated in Fig. 5a.This landscape consists of the following elements: (i) west-east stretching trough-shaped valley with an amphitheatre-shaped head, inclined towards the east in its lower part; (ii) dam with a trapezoidal cross section running across the valley, consisting of 100 % solid material; (iii) reservoir impounded by the dam; (iv) landslide release mass near the northwest corner of the area of interest (Landslide 1); (v) landslide release mass directly north of the dam (Landslide 2); (vi) hydrograph release of landslide near the southwest corner of the area of interest; (vii) measurement profile for output hydrograph downstream from the dam.Both landslide release masses assume the shape of a distorted hemi-ellipsoid imposed on the basal topography (see Fig. 5a).The algorithm for exactly reproducing the generic landscape in GRASS GIS is available at http://www.avaflow.org/casestudies.html.

Modelling strategy and parameterization
Landslides 1 and 2 consist of 75 % solid and 25 % fluid by volume (uniformly mixed); the input hydrograph I1 (see Fig. 5b) consists of 50 % solid and 50 % fluid per volume.The parameters and settings applied are summarized in Tables 2 and 3.
Three computational experiments are performed, with increasing complexity from A to C: -Experiment 1A: Landslide 1 is released and interacts with the reservoir.The dam is assumed stable and may therefore not be entrained.
-Experiment 1B: Again, Landslide 1 is released and interacts with the reservoir.However, dam material is allowed to be entrained in this experiment.
-Experiment 1C: Landslide 2 is released and interacts with the dam and the reservoir.The release from the input hydrograph I1 starts after 10 s and continues for a period of 130 s (see Fig. 5).Dam material is allowed to be entrained at all stages of the computational experiment.
All experiments are performed at a cell size of 10 m and for a duration of t term = 300 s; t out = 5 s.The solid and fluid discharges are continuously recorded at the output hydrograph profile O1 downstream.The stopping function is deactivated (see Table 3).

Results
Animations illustrating the time evolution of the flow heights in all three experiments are enclosed in Animations 1A, B, and C in the Supplement.
Figure 6a-f illustrates the flow heights at selected points of time during Experiment 1A.Landslide 1 (see Fig. 5a) impacts the backward portion of the reservoir after few seconds and generates a water wave -oblique and perpendicular to the impact -that overtops the dam from t = 50-55 s onwards.The output hydrograph O1 starts recording discharge at t = 65 s, with the peak of the first major flood wave passing at t = 75 s (Q f = 8 × 10 4 m 3 s −1 ; Fig. 6g).We note that the discharge and the flow height recorded by the hydrograph do not strictly follow the same pattern, as the discharge relates to a profile and the flow height relates to a point (see Fig. 2).Meanwhile, the impact wave is deflected at the dam and alleviates slowly.Further overtopping events caused by multiple deflections of the alleviating wave occur mainly at the marginal parts of the dam at t = 110, 150, 160, 200, and 270 s, leading to smaller peaks in the output hydrograph (Q f = 1.5×10 4 m 3 s −1 at t = 175 s; Q f = 2.2×10 3 m 3 s −1 at t = 285 s).The solid content passing the hydrograph profile is almost negligible as all solid landslide material remains in the reservoir basin.At t = 300 s, the impact wave in the lake has almost alleviated (see Animation 1A in the Supplement).
Experiment 1B (Fig. 7) is identical to the Experiment 1A until the point when the impact wave reaches the dam at t = 50 s.Entrainment of the dam starts with overtopping, which sets on at the lateral portions.Part of the dam is entrained during overtopping by the initial impact wave.Whilst massive outflow from the reservoir occurs due to the decreased level of the dam crest, part of the wave is deflected at the dam and pushed back towards the rear part of the reservoir, inducing a system of secondary waves.The remaining dam material is entrained when hit by those secondary waves.At t = 200 s, the entire dam has disappeared and the reservoir starts emptying completely.In contrast to Experiment 1A, due to the emptying process, the system does not approach a static equilibrium after t = 300 s (see Animation 1B in the Supplement).
The temporal patterns of the simulated entrainment and wave propagation are clearly reflected in the discharge recorded at the output hydrograph O1 (see Fig. 7g).As a consequence of dam overtopping, fluid discharge at O1 starts increasing at t = 65 s and reaches a first peak at t = 80 s (Q f = 5.1 × 10 4 m 3 ).Solid discharge -a consequence of entrainment of the dam -starts slightly delayed, reaching a first peak roughly 10 s later (Q s = 2.1 × 10 4 m 3 s −1 ).A depression in both of the discharge curves at t = 155-160 s indicates that the initial impact wave has passed through.A second, larger peak of fluid discharge is simulated at t = 195 s (Q f = 1.0 × 10 5 m 3 s −1 ).It occurs synchronously with a second, comparatively smaller peak of solid discharge (Q s = 2.1 × 10 4 m 3 s −1 ), indicating a high degree of mixing of the solid and fluid components of the flow.The pronounced second peak of Q f is a consequence of the secondary waves in combination with the lowered level of the dam.After the peak, Q s slowly and unsteadily decreases (the entire dam has been entrained and the material has passed through), whilst Q f remains high.Due to the entrainment of the dam, the simulated discharges are much higher than those computed in the Experiment 1A (see Fig. 6g).
In Experiment 1C (Fig. 8), Landslide 2 impacts the dam and the frontal part of the reservoir less than 10 s after release.The proximal portion of the dam is entrained rapidly.The eastern part of the landslide moves outside of the reser- voir in downstream direction.Consequently, the solid discharge at the output hydrograph O1 starts at t = 40 s, reaching a peak of Q s = 2.4 × 10 4 m 3 s −1 5 s later (see Fig. 8g).Due to the high (75 %) solid fraction of the landslide, the fluid discharge is lower at that time (Q f = 1.5 × 10 4 m 3 s −1 ).The western part of the landslide interacts with the reservoir, causing overtopping at the south (distal) portion of the dam.This results in the increase of fluid discharge recorded at O1, culminating at t = 60 s when the solid discharge has already passed its peak (Q f = 3.7 × 10 4 m 3 s −1 ).The resulting impact at O1 has reduced after t = 105 s in terms of discharge, even though the total flow height remains at H > 15 m.This means that the flow material moves slowly at O1.
From t = 35 s onwards, the flow released through the input hydrograph I1 (see Fig. 5b) pushes the reservoir water towards the northeast.The southern part of the remaining dam is overtopped by the resulting inhomogeneous solid-fluid mixture (including material originating from Landslide 2), leading to substantial further entrainment.In contrast to Experiment 1B, however, the dam is not completely entrained.The wave increasingly influences the discharge recorded at O1, leading to a peak at t = 180 s (Q s = 6.9 × 10 3 m 3 s −1 ; Q f = 1.7 × 10 4 m 3 s −1 ).At that time the hydrograph indi- cates a well-mixed flow with α s ≈ 0.25, composed of fluid from the reservoir, solid-fluid mixtures from the landslide and the hydrograph release, and solid material from the dam (see Fig. 5a).The solid and fluid discharges remain at an almost constant level thereafter, reflecting a steady emptying of the reservoir.

Event description
The Acheron rock avalanche in Canterbury, New Zealand (Fig. 9), occurred approximately 1100 years BP (Smith et al., 2006).It is characterized by sharp bending of the flow path, a limited degree of spreading into the lateral valleys, and a high mobility (travel distance: 3550 m; measured angle of reach: 11.62 • ).It was used as a test event for the computational tool r.randomwalk (Mergili et al., 2015).We employ a 10 m resolution digital elevation model (DEM) derived by stereo-matching of aerial photographs.ODA and OIA are derived from field and imagery interpretation as well as from data published by Smith et al. (2006).The OIA possibly underrepresents the real impact area, as it might exclude some lateral and run-up areas of the rock avalanche which are not any more recognizable as such in the field.The distribution of release and deposition heights and an estimated release volume of 6.4×10 6 m 3 are deduced from the reconstruction of the pre-event topography.According to this reconstruction, the maximum release height is 78.5 m whilst the maximum deposition height is 25.9 m.

Modelling strategy and parameterization
Preliminary tests have shown that the simulation results of r.avaflow are potentially sensitive to variations in the initial solid fraction α s0 and the basal friction angle δ, parameters which are uncertain in many real-world applications.We perform two computational experiments for the Acheron rock avalanche: 1. Experiment 2A: III and DII are computed from a set of 121 model runs.Thereby, α s0 is varied from 0.5-0.9, and δ is varied from 15 to 25 • (see Table 2).The variation is done in a controlled way assuming a uniform probability density function; i.e. a regular grid with 11 grid points in each dimension is laid over the two-dimensional parameter space.III is then evaluated against the OIA, and DII is evaluated against the ODA.α s0 and δ are optimized in terms of L, FoC, CSI, and D2PC derived from H C and the ODA.
Both experiments are conducted at a cell size of 20 m.Entrainment is not considered whilst stopping and deposition are included (see Table 3).All flow parameters except for δ are kept constant (see Table 2).

Results
Figure 10 illustrates III and DII derived with the parameter settings shown in Tables 2 and 3 (Experiment 2A).AUROC is 0.830 with regard to III and 0.838 with regard to DII.In general, those areas with high values of III coincide with the OIA, whilst those areas with lower values of III lie close to the margins or outside of the OIA.The performance of III suffers from the motion of small portions of the simulated avalanche in the wrong (northern) direction and from excessive lateral spreading and run-up in the upper part, observed for all tested combinations of α s0 and δ (high values of III; see Fig. 10a).However, one has to consider that the event occurred hundreds of years ago and run-up may have occurred even though it is not any more recognizable in the field and therefore excluded from the OIA.High values of DII are fairly constrained to those cells within the ODA (see Fig. 10b) which is most probably better defined than the OIA.Those areas with lower, but non-zero values of III or DII both extend well beyond the reference areas.Particularly the travel distance appears highly sensitive to the choice of α s0 and δ.
We now focus on the DII map and evaluate the performance of the deposition maps simulated with the various combinations of α s0 and δ against the ODA. Figure 11 illustrates the dependency of the model performance (defined by the parameters summarized in Table 4) on the combination of α s0 and δ employed for a given model run.All four parameters clearly indicate that, within the ranges tested, the model results are sensitive to both δ and α s0 .L, CSI, and D2PC display their optima near δ = 17 • as long as α s0 ≥ 0.7.With higher fluid content, the optimum value of δ increases, arriving at 20 • with α s0 = 0.5 (see Fig. 11a, b, and d).This pattern appears plausible as far as a higher fluid content is supposed to increase the mobility of the flow, compensating for higher values of δ.However, values of α s0 < 0.7 are not plausible at all for rock avalanches of this type.For α s0 ≥ 0.7, FoC displays its optimum of 1.0 at δ ≥ 21 • , depending on α s0 .FoC ≈ 1.25 for the value of δ where the other parameters reach their optimum (see Fig. 11c).This would be fine for many applications in practice where slightly conservative results are desirable.
Consequently, we consider δ = 17 • and α s0 = 0.8 -in addition to the parameter values given in Table 2 -one possible combination for back-calculating the Acheron rock avalanche.The simulation is repeated with exactly this combination (Experiment 2B).We note, however, that these parameter values do not necessarily represent the real-world conditions, as the fluid content of rock avalanches may be much lower.Figure 12 shows the maps of H Max and H C , both corresponding reasonably well to the OIA and the ODA, respectively.The slightly larger simulated than observed deposit (see Fig. 12b) corresponds to FoC ≈ 1.25; the almost perfect correspondence of the observed and simulated termini corresponds to L ≈ 0. This means that the fact that the result is rather conservative than non-conservative (FoC > 1) relates to lateral spreading rather than to the travel distance of the rock avalanche.Animation 2 in the Supplement illustrates the time evolution of the flow height in Experiment 2B.

Discussion
The key purpose of the present article is to provide a general introduction to the key functionalities of the computational tool r.avaflow.Thereby, the simulated patterns of flow height in Experiment 1 (see Sect. 3.1) are plausible, and the agreement of the observed and simulated deposition areas in Experiment 2B (see Sect. 3.2) appears reasonable.Yet, these experiments can neither replace model validation with observed process chains or interactions, nor can they replace thorough multi-parameter sensitivity analysis and optimization efforts, which will both be the subjects of future research.Fully documented two-phase process chains with readily available pre-and post-event DTMs are scarce.Preliminary r.avaflow results for the 2012 Santa Cruz multilake outburst flood in the Cordillera Blanca, Peru (Emmer et al., 2016), are however promising.
Experiment 2 serves for the demonstration of the parameter sensitivity analysis and optimization functions of r.avaflow.The outcomes may be different when changing the cell size or any of the flow parameter values (see Table 2).Making r.avaflow fit for forward predictions will require a thorough multi-parameter sensitivity analysis and optimization campaign involving a large number and variety of well-documented events.Thereby, we aim at obtaining guiding parameter values -or, more appropriately, guiding parameter ranges -for mass flow processes of different types and magnitudes.Approaches to perform such analyses are readily available, and some of them can be directly coupled to r.avaflow (Fischer, 2013;Fischer et al., 2015;Aaron et al., 2016;Krenn et al., 2016).However, due to the complex nature of two-phase mixture flows, r.avaflow depends on a relatively large number of flow parameters, a fact that rep- resents a particular challenge in terms of the computational resources as well as in terms of visualization and interpretation of the results of multi-parameter studies.
r.avaflow represents a modular framework, allowing for the future enhancement of its particular components.One issue concerns the numerical implementation of the twophase model equations, combining topography-following coordinates with the quadratic cells of the raster data given in GIS coordinates (see Sect. 2.3).As in comparable simulation tools (e.g.Christen et al., 2010a, b;Hergarten and Robl, 2015), approximations are currently used for coordinate transformation in r.avaflow.This issue is closely related to the fact that the model equations that are commonly exwww.geosci-model-dev.net/10/553/2017/Geosci.Model Dev., 10, 553-569, 2017 pressed in topography-following coordinates are hardly compatible with the data given in GIS coordinates.
A detailed and fully discrete description of the TVD-NOC scheme exists in the literature (Wang et al., 2004), and the scheme served well for various theoretical test cases (e.g.Pudasaini et al., 2014;Kafle et al., 2016;Kattel et al., 2016).However, we also identify two major drawbacks: -Although the numerical scheme itself should be shock capturing and volume preserving (Tai et al., 2002;Wang et al., 2004), these properties may not fully hold in practical applications (i.e.bounded gravitational mass flows with well-defined margins over complex topography).
The complementary functions with ID 1-3 introduced in Sect.2.4 partly compensate for the issues raised.
-For real flow applications, full handling of the evolution of the basal topography is not straightforward: the TVD-NOC scheme may introduce diffusion even though the evolution of the basal topography is not a standard transport equation.Entrainment is therefore, as a first step, included as a complementary function.
The numerical scheme employed will have to be enhanced to directly and effectively incorporate the complementary functions outlined in Sect.2.4 in a fully consistent way.Extensions of similar schemes have been tested for generic examples (e.g.Zhai et al., 2015) and could serve as a valuable basis also to implement a mechanical model for erosion, entrainment, and deposition (Pudasaini and Fischer, 2016).On the one hand, such an erosion model may build on existing concepts (e.g.Fraccarollo and Capart, 2002;Sovilla et al., 2006;Medina et al., 2008;Armanini et al., 2009;Crosta et al., 2009;Hungr and McDougall, 2009;Le and Pitman, 2009;Iverson, 2012;Pirulli and Pastor, 2012).On the other hand, it may further require some fundamentally new ideas with regard to deposition.

Conclusions and outlook
We have introduced r.avaflow, a multifunctional open-source GIS application for simulating two-phase mass flows, process chains, and interactions.The outcomes of two computational experiments have revealed that r.avaflow (i) has the capacity to simulate complex solid-fluid process interactions in a plausible way, and (ii) after the optimization of the basal friction angle and the solid content of the release mass, reasonably reproduces the observed deposition area of a documented rock avalanche.However, it was out of the scope of the present work to validate the results obtained for complex process interactions against observed real-world data or even to conduct a comprehensive multi-parameter optimization campaign.Such efforts will be the next step towards making r.avaflow ready for the forward prediction of possible future mass flow events.Thereby, we will attempt to estab-lish guiding parameter values for different types of processes and process magnitudes.
At the same time, we have identified a certain potential for the future enhancement of some the components of r.avaflow.
The key challenges will consist in (i) integrating the model equations in an up-to-date numerical scheme, allowing to directly include the complementary functions, and (ii) replacing the empirical entrainment model with a mechanical model for entrainment and deposition.

Code availability
The model codes along with a user manual are available at http://www.avaflow.org/software.html(Mergili et al., 2017).

Data availability
The scripts, the text file, and the GRASS locations with the spatial data required for reproducing the computational experiments described in Sect. 3 are available at http://www.avaflow.org/casestudies.html(Mergili and Krenn, 2017).
The Supplement related to this article is available online at doi:10.5194/gmd-10-553-2017-supplement.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Logical framework of r.avaflow.The transformations and retransformations refer to the conversion of heights and GIS coordinates to depths and topography-following coordinates, and vice versa (see Sect. 2.3).

Figure 2 .
Figure 2. Sketch of a hydrograph profile.The flow surface of input hydrographs is defined by H P and is extended in cross-profile direction either to the edge of the profile or until it intersects with the basal topography.
and M sy and M fy in y direction (M sy = D s •v sy , M fy = D f •v fy ), where v is flow velocity.

Figure 3 .
Figure 3. Interactions of the flow with the basal topography: (a) entrainment, assuming that H Emax,s and H Emax,f are not limiting; D i : total initial flow depth (s + f); M i : total initial momentum (s + f); D E : entrained depth; M E : total increase in momentum due to entrainment (s + f).Panel (b) indicates stopping and deposition.Both panels represent sections along the steepest slope of the basal topography.Note that stopping and deposition usually occur on less inclined slopes than drawn in (b) which represents upslope movement.

Figure 4 .
Figure 4. Validation of r.avaflow results.(a) Validation scores for single model run; (b) multiple model runs: threshold levels of III or DII, employed to produce (c) ROC curves.

Figure 5 .
Figure 5. Generic landscape used for Experiment 1A-C.(a) Oblique view illustrating the topography and elements of the landscape.(b) Input hydrograph I1 employed for Experiment 1C.

Figure 6 .
Figure 6.Key results of Experiment 1A.(a-f) Sequence of simulated flow heights and solid ratios at selected points of time; see Animation 1A in the Supplement for animations of flow height sequences; (g) output hydrograph O1 (see Fig. 5a).

Figure 7 .
Figure 7. Key results of Experiment 1B.(a-f) Sequence of simulated flow heights and solid ratios at selected points of time; see Animation 1B in the Supplement for animations of flow height sequences; (g) output hydrograph O1 (see Fig. 5a).

Figure 8 .
Figure 8. Key results of Experiment 1C.(a-f) Sequence of simulated flow heights and solid ratios at selected points of time; see Animation 1C in the Supplement for animations of flow height sequences; (g) output hydrograph O1 (see Fig. 5a).

Figure 9 .
Figure 9.The Acheron rock avalanche.(a) Oblique view; the view point is indicated in (b) illustrating the location and the main elements of the rock avalanche; ORA: observed release area.

Figure 10 .
Figure 10.Results of Experiment 2A: (a) impact indicator index III and (b) deposition indicator index DII derived for the Acheron rock avalanche.

Figure 12 .
Figure 12. Results of Experiment 2B.(a) Maximum flow height H Max ; (b) height of final deposit H D (as entrainment is not considered, H D = H C ).Note that, due to the predominance of solids, the bluish and greenish colours indicated in the legend do not appear in the map (see Figs. 6-8).

Table 1 .
Key input and output parameters of r.avaflow -s: solid; f: fluid; t: total.Remarks: 1: mandatory; 2: one of the input data sets A, B, or C + D is mandatory, C + D may also be provided in addition to A or B; n D ≥ n C , if n D > n C the remaining sets of D are output hydrographs; 3: either A or B may be provided if entrainment is activated, otherwise all values of H Emax = ∞; C is mandatory with entrainment; 4: at least one of the data sets A, B, and C is mandatory for validation.
Output (excluding validation and visualization output; see Sect.2.6) Maximum flow height, kinetic energy, and pressure (each for s, f, t) H Max , T Max , p Max m, J, Pa Raster maps Always Flow height, flow kinetic energy, and flow pressure at each output time step t out (each for s, f, t) H tout , T tout , p tout m, J, Pa Raster maps Always Flow velocities in x and y direction, and in absolute values (each for s, f) v x , v y , v m s −1 D -n C output hydrograph tables: flow heights, velocities, and discharges at defined points of time (s, f)

Table 2 .
Flow Pudasaini (2012)ntrainment coefficient required with the enhanced version of thePudasaini (2012)two-phase flow model.Exp. 1 and 2 refer to the values used for the computational experiments introduced in Sect.3.