Flow-Py v1.0: A customizable, open-source simulation tool to estimate runout and intensity of gravitational mass ﬂows

. Models and simulation tools for gravitational mass ﬂows (GMF) such as snow avalanches, rockfall, landslides and debris ﬂows are important for research, education and practice. In addition to basic simulations and classic applications (e.g., hazard zone mapping), the importance and adaptability of GMF simulation tools for new and advanced applications (e.g., automatic classiﬁcation of terrain susceptible for GMF initiation or identiﬁcation of forests with a protective function) are currently 5 driving model developments. In principle, two types of modeling approaches exist: process-based physically motivated and data-based empirically motivated models. The choice for one or the other modeling approach depends on the addressed question, the availability of input data, the required accuracy of the simulation output, and the applied spatial scale. Here we present the computationally inexpensive open-source GMF simulation tool Flow-Py. Flow-Py’s model equations are implemented via the Python computer language and based on geometrical relations motivated by the classical data-based runout angle concepts 10 and path routing in three-dimensional terrain. That is, Flow-Py employs a data-based modeling approach to identify process areas and corresponding intensities of GMFs by combining models for routing and stopping, which depend on local terrain and prior movement. The only required input data are a digital elevation model, the positions of

We present the innovative and educational Flow-Py simulation tool, which employs a data-based motivated approach to 70 predict the magnitude, i.e., runout (spatial extent including starting, transit and runout zones) and intensity (effects of a GMF at a specific location) of GMF processes. Flow-Py builds on the ideas and algorithms from existing data-based GMF models.The Flow-Py algorithm is based on a flow path identification in three-dimensional terrain (routing) and concepts for runout and intensity estimates along this path (stopping). To determine the GMF's runout and intensity we utilized well known runout (travel) angle concepts (Heim, 1932), and derived corresponding geometrical quantities to motivate the Flow-Py model equa-75 tions. These geometric relations serve further as reference to validate the Flow-Py implementation and results. In addition to runout and intensity predictions, Flow-Py simulations results are also a measure of how exposed a location in the flow path is regarding the number of starting zones and associated transit zones, which route flux through that location.
This contribution is structured as follows: In Sect. 2 we describe the motivation and implementation of our GMF model, which is further explained in the code repository . A validation experiment is presented in Sect. 3 80 which shows simulation results from three simple generic slopes. The performance of Flow-Py is tested via a regional scale simulation of snow avalanches in Sect. 4. The customization of Flow-Py is described in Sect. 5 and shows how flexible the simulation tool is and that it can be easily adapted with extensions to specific modeling questions.
With this contribution we enable the reader to reproduce and understand the basic GMF model concepts and their implementation in the Flow-Py code.

85
The main objectives of the Flow-Py simulation tool are to compute the spatial extent (hereafter referred to as runout) of GMFs, which consists of starting, transit and runout zones, and the intensity of the GMF. Flow-Py is based on data-driven empirical modeling ideas (Heim, 1932) with automated path identification (Holmgren, 1994;Horton et al., 2013;Huber et al., 2016;Wichmann, 2017) to solve the routing and stopping of GMFs in three-dimensional terrain. Data-based models often require less 90 input data, a less complex parameterization and solution (e.g. no time-dependent equations are usually solved) than processbased models. The Flow-Py simulation tool has been designed as a computationally inexpensive data-based model, which facilitates its application on regional scales, including a large number of GMF paths. Simulations of single starting cells take 1 to 10 seconds where process based, depth average simulations usually operate under the order of minutes. This can be attributed to the fact that no time depend equations, which process based models are built on, are solved in the underlying model equations 95 of Flow-Py. The Flow-Py code is written in the Python computer language taking advantage of Pythons object-oriented class method. The well-structured model implementation allows users to address GMF specific modeling questions by keeping the parameterization flexible and enabling to include customized model extensions and add-ons. Flow-Py has already been applied to dry snow avalanches, rockfall and shallow-seated landslides by adapting the parameterization. Experience from similar studies also suggests that the model may also be suitable for other GMFs such as debris flows and wet snow avalanches 100 Holmgren (1994); Gamma (1999).
The development philosophy to maximize the applicability of Flow-Py builds on: 1. flexible yet minimal input data requirements, 2. simple parameterizations which can describe a range of GMFs, and 3. a highly adaptable and customizable source code.

105
In the following sections the model motivation, implementation, input data and Flow-Py results, and underlying model equations are explained in detail.

Model motivation
The Flow-Py's routing and flow path identification in three-dimensional terrain was inspired by the gravitational process path model CPP (Wichmann, 2017), which introduced a weighting factor for the flow direction, and the programming architecture 110 and persistence equations of Flow-R (Horton et al., 2013), combined with an adapted version of the flow direction algorithm (Holmgren, 1994) to appropriately model movement in flat and uphill terrain. The routing is based on local terrain and prior movement (flow direction and process intensity), which determines the flow path from starting to transit and runout zones and simultaneously describes the flow concentration, including lateral spreading. To estimate the process intensity along the identified path and the runout by introducing a stopping criterion we utilize the well-known runout angle (α) concept (Heim, 115 1932;Lied and Bakkehøi, 1980;Bakkehøi et al., 1983;Körner, 1980) and derived corresponding geometrical quantities to Figure 1. GMF path with altitude z(s), projected travel distance s and local slope angle ψ with starting point s0, z(s0) and runout point sα, z(sα). The corresponding geometric quantities are directly related to the runout angle α concept and include the local travel angle γ, with corresponding total altitude change zγ and the process intensity measure z δ with angle δ. motivate the Flow-Py model equations. Figure 1 depicts the runout angle along with the corresponding geometric relations in a two-dimensional representation along a GMF path, building the foundation for the underlying model equations.
The geometric relations are directly deduced from the runout angle α and allow to motivate the stopping and intensity estimates. Additionally, the geometric solution, represented by the α-line from starting (s 0 , z(s 0 )) to the runout (s α , z(s α )) 120 points, serves as reference for a model validation. Important quantities include the local travel angle γ: which, at the end of the GMF path, corresponds to the total travel angle (i.e., the so-called runout angle α (Heim, 1932)), that can be expressed as: The local travel angle height z γ corresponds geometrically to the total elevation drop z γ from the starting point s 0 to the currently projected runout length s along the path: The total elevation drop z γ splits into z α 130 which is associated with the dissipation kinetic energy height, and z δ which is a measure of the process intensity, corresponding to the kinetic energy height (based on the principles of energy conservation, assuming a block movement with frictional dissipation associated to a Coulomb friction, Heim, 1932).

135
The Flow-Py simulation tool is implemented based on object-orientated programming ideas, which allows for easy model customization . Flow-Py is written in the freely available modern programming language Python3 (Van Rossum and Drake, 2009), which is widely used and supported by an active online community. The simulation tool is highly adaptable and different routing and stopping routines can be easily implemented, which enables the user to adjust the parameterization, also for multi model runs, and the equations that govern the movement of the mass down slope, and to im-140 plement Flow-Py in model chains. Flow-Py can be run either by command line allowing it to be called by external programs or in a BASH file, or with a simple GUI, which guides the user through choosing input files and the parameterization.
A GMF usually has one or more starting zones that span over a single or multiple starting cells. Flow-Py computes the socalled path, which we define as the spatial extent of the routing from each starting cell to the stopping cells. Each starting zone 145 is associated with its own unique path; however, a certain location in the terrain can belong to many paths. Flow-Py identifies the path with spatial iterations on the cell level, starting with a single cell of a starting zone and then transferring the final results of the cell and path levels to the output raster level (see Fig. 2). To route on the three-dimensional terrain operating on a quadrilateral grid, we implemented the geometric concepts that have been introduced in Sect. 2.1. That is, each path calculation starts with a starting cell, operating on the cell level, requiring the definition of parent, base, child and other neighbor cells (see  queue will be the base cell (center cell) for subsequent calculations. Information on the iteration step is temporally stored to the respective cell's Flow-Class. When the path calculations are finished, values from each cell's Flow-Class are updated to their respective location in the result array such that either a maximum value (e.g., Z δ max ,the maximum Z δ for a cell over all path calculations) or a running sum (e.g.,Z δ sum , the sum of all Z δ values for a cell over all path calculations) of all calculated paths 160 that route flux through that cell are stored. The Flow-Class can be extended to store additional information that can be used to adjust stopping and routing calculations, e.g., the runout angle α is saved in the Flow-Class and could be adapted and scaled with Z δ b to account for an energy dependent friction. Using Python's object-oriented class method is a major advantage for advanced users since they can easily develop custom extensions or add-ons. We present an example for a back-tracking extension, which saves information of infrastructure located 165 in GMF flow paths in the Flow-Class adapting the Flow-Py output in Sect. 5.

Input data
Flow-Py's core function loads and handles all input data, which are a digital elevation model (DEM) and a release raster in .asc or .tiff format. The release raster shows observed or potential GMF starting zones containing one or several starting cells. The release raster can be created by vector-to-raster conversion of polygon mappings by expert or by onset-susceptibility 170 modelling :::::::: modeling. Flow-Py employs parallel processing for short model run times by splitting the release raster and DEM into tiles. Each tile is solved independently and sequentially in its own dedicated computer core and processing threads. Multiprocessing is set as the Flow-Py default, i.e. the number of free cores and the amount of RAM are first checked before splitting the starting cells and spreading runout calculations among the free computer cores, making sure that the amount of RAM is not exceeded. The individual calculations are merged by updating the result arrays, which are transformed into output raster.

175
The release raster shows potential GMF starting zones containing one or several starting cells. The DEM and the release raster must be in the same extent and resolution with no resolution limit; however, 5 m and 10 m raster resolutions have been tested. The major differences between different types of GMFs in regional runout modeling is the behavior of the movement and its runout, which can be summarized by the runout length and the convergence or divergence of the spreading movement.
These behaviors are controlled by the parameterization of the stopping and routing routines in Flow-Py.

Model equations and path identification
The Flow-Py model equations are formulated with respect to an equidistant quadratic grid with the same resolution and extent as the input raster. During each spatial iteration, calculations are made on a 3x3 cells subset of the raster, where the flux across the base cell (subscript b) is solved (see Fig. 2). The eight neighbor cells to the base cell (subscript i and n) can be parent cells (subscript p) during an iteration step acting as flux source cells, or child cells (subscript c) acting as flux sink cells.

185
The governing runout modeling question is broken down into two subquestions: 1. Where does the GMF move to? 2. Where does the GMF stop?
These questions are addressed in two dedicated modeling routines called the routing routine and the stopping routine.

190
The routing routine considers a terrain contribution T i and a contribution accounting for prior motion called persistence P i (Horton et al., 2013); the flux is solved from parent cells through the base cell to child cells. Eq. (6) is the basis of the routing algorithm and shows how the terrain contribution T i and the persistence contribution P i are combined to distribute the routing 195 R i is the routing flux from the base cell to neighbor cell i and R b is the total routing flux into the base cell (for starting cells To conserve R i , the amount of R b must be equal to R i unless a stopping criteria is met (see Sect. 2.4.4). To conserve flux, T i and P i to cell i are normalized across all neighboring cells n. The normalized direction is then scaled with R b .

Terrain-based routing 200
The terrain-based routing accounts for the guiding effect of the slope on the movement. To distribute the flux we utilize the terrain routing function: where T i is the normalized terrain based routing from the base cell i and To control the concentration of routing flux an approach based on the multiple flow direction algorithm for runoff has been 210 employed (Holmgren, 1994). The exponent exp together with the flux cutoff (see Sect. 2.4.6) controls the lateral spreading of the flow (Horton et al., 2013). When exp increases, the terrain based routing flux is concentrated to the steepest decent.
Together with the flux cutoff > 0, this results in the path's lateral spreading to be reduced. As exp → ∞ the divergence results in a single flow direction (block movement) and as exp → 1 wide spreading is encouraged (fluvial movement). However, other terrain-based routing approaches can be easily implemented in Flow-Py (see Horton et al. (2013) for summary).

Persistence-based routing
The persistence-based routing contribution aims to account for the influence or prior GMF movement on the subsequent routing. It must be noted that persistence is empirically derived and may be conceptually comparable to momentum; however, Flow-Py's underlying model equations do not account for mass (and hence momentum).
Equation (8) shows the persistence routing function P i for neighbor cell i: which is consists of two components, the direction D n Eq. (9) and the intensity Z δ p , which has classically been called energy line height (Körner, 1980). Because a base cell can receive flux from many parent cells p the persistence routing function is calculated over all neighbor cells n considering the incoming flux from each parent cell.
The direction D n maintains the flow direction from a parent cell (p, flux source) to the base cell b. Weights are used to define 225 the flow direction and are expressed as: where θ = ∠ pbn − 180 • is the resulting deviation angle between the projected incoming and outgoing direction, with ∠ pbn as angle formed from the directions of parent ::: cell p to base ::: cell b and from base ::: cell : b to neighbor ::: cell ncells, compare  Horton et al. (2013). The reason that the persistence function passes flux through three cells and not only one is to compensate for the restriction that there are only eight directions to move on a raster grid. A weight of 0 is given to all other cells including the parent cell.
The intensity Z δ p is stored in the Flow-Class of the parent cell p from a previous iteration step, and the value of Z δ is saved in the Flow-Class of each child cell. If one child cell has more than one parent cell, then Z δ max,path (maximum value of Z δ for 235 the many combinations of routes to a cell on a single path) is stored in its Flow-Class.
The intensity Z δ n at the neighbor cell n is cell-wise calculated, i.e. cell to cell throughout the spatial iterations. The intensity Z δ bn refers to the iterative part of Z δ n that is associated with the spatial step from the base cell b to the neighbor cell n. Equation (10) shows the calculation of Z δ n , where Z δ b is the intensity of the base cell b, which is stored in its Flow-Class. Z δ b was calculated on a previous spatial iteration when the current base cell was a child cell: where Z δ bn is calculated with respect to:

245
where the subscript bn refers to base cell b to neighbor cell n, with the distance S bn and the iterative energy quantities Z α bn , Z γ bn , Z δ bn (see Fig. 1). The total projected distance along the GMF path S n is expressed as The parent cell further away from the stopping condition (larger Z δ ) will have more influence on the routing flux. After all 250 n parent cells are calculated for each neighbor i the persistence-based routing P i is combined with the terrain-based routing T i , as seen in Eq. (7). When the parent cell has a large Z δ p the persistence-based routing P i will be the dominant term in Eq. (7); however, if Z δ p is small, then the terrain based-routing T i will dominate the routing direction. There are two limits that are imposed in the persistence routing routine: First, any cell that has previously been a base cell cannot be a child cell (a parent cell can not be a child cell). The disadvantage of this limit exerts on half-pipe shaped terrain 255 in which the mass moves up a slope and back down on the same path but in the opposite direction. This limit is necessary to keep small amounts of flux from routing back and forth in terrain shaped like a bowl. The major advantage of this limit is the reduction of iteration steps by not calculating further flux for child cells resulting from flux oscillating in a bowl feature.
The second limit is imposed on the maximum value of Z δ i , which is a limit of the process intensity (Z δ lim ), corresponding to a kinetic energy height or GMF velocity limit, respectively: which is important for some GMF types, because it is analogous to introducing a turbulent friction coefficient in a processbased model (Horton et al., 2013).
In the examples used in Sect. 3, 4 and 5, no such limits are imposed.

Stopping
Two stopping criteria are employed: First, a runout angle criterion that limits how far the GMF runout goes. The second is a flux 265 cutoff stopping routine, which, together with the divergence control (exp) in the routing routine, limits the lateral spreading of the path. The GMF will not propagate further, if either stopping criteria is met; however, the runout angle mainly determines the total travel distance in the main flow direction, while the flux cutoff influences the lateral spreading.

Runout angle-induced stopping
The runout angle-induced stopping routine is based on the geometric quantities derived with the α angle concepts (c.f. Eq. (1 270 to 5); see Fig. 1). The local travel angle γ n is the inclination of the line formed from the top of the starting zone to the current neighbor cell n. The stopping condition is reached when γ n < α, i.e. when Z δ n < 0.
When the stopping condition is met, no child cells are assigned in the next iteration step.

Routing flux-induced stopping 275
The second stopping criterion is based on the assumption that a GMF must have a critical amount of routing flux R i to continue its propagation. If the GMF has an excessively divergent flow concentration that dilutes down and across a slope, then the flow concentration (that can be associated to GMF mass) disappears at a critical amount of spreading, corresponding to the critical routing flux threshold R stop .
The routing flux stopping criteria is met when and the runout angle stopping condition is also met and neighbor cell i is not a child cell. If R i ≥ R stop and the runout angle stopping condition is not met, then neighbor cell i is a potential child cell, and is added to the calculation queue and a Flow-Class is accepted.
The routing flux-induced stopping mainly limits the width or spreading of the path. The magnitude of the routing flux of the 285 potential child cell R i relates to the percentage of initial routing flux from the start cell, where the starting flux R start = 1. As default R stop = 3 · 10 −4 has been adopted in Flow-Py and is shown in the examples in Sect. 3, 4 and 5.
2.5 Flow-Py Outputs ::::::: outputs The outputs of Flow-Py are a set of raster in the same resolution and extent as the input DEM providing information about the runout of the GMF and different measures of the intensity:

290
-Z δ max is the local maximum Z δ for a cell over all path calculations. This is a geometric measure of highest intensity in terms of Z δ for all starting cells which can be associated to maximum kinetic energy that is expected at each location (raster cell).
-R max is the local maximum routing flux for a cell over all path calculations. This is a measure of intensity in terms of the maximum of flow concentration from a single start cell that is expected at each location (raster cell).

295
-Z δ sum is the sum of all Z δ for a cell over all path calculations. This is a measure of intensity in terms of Z δ combined with the number of starting cells that route flux through a location (raster cell).
-CC is path cell counts, which is the number of paths that route flux through a location (raster cell). Together with Z δ sum an average of Z δ can be formed.
γ max is the local maximum flow path travel angle for a cell over all path calculations. This is a measure of how exposed 300 a location is with regards to how close the highest GMF intensity in terms of Z δ is to the runout angle stopping criteria.

Model testing and validation on generic slopes
This first computational experiment demonstrates the Flow-Py routing and stopping algorithms for GMF modeling on simple but increasingly complex generic topographies. We highlight how GMFs interact with different terrain features and show the influence of different parameterizations on the flux; however, we do not perform a detailed parameter study, which is 305 beyond the scope of this contribution. First, we describe the scenarios (terrain and model parameterizations) and present the simulations results. For each scenario, we altered the model parameterization or terrain complexity. Then the behavior of the simulations and a comparison to the geometrically expected results, which allows for validation of the models implementation, are summarized and discussed.
The generic topographies used for Flow-Py testing were generated using the generate topography functions provided within 310 AvaFrame (Wirbel et al., 2021). Terrain data was saved in ASCII raster format (.asc) with 10m resolution. The release raster consisted of three 100 m 2 -neighboring starting cells (starting zone = 300 m 2 , 3 raster cells) located close to the top of the generic terrain model at an elevation of 982 m. The three starting cells are centered on the y-plane.

325
Comparing the top and bottom panels of Fig. 3, it can be seen that keeping the terrain, the runout angle and R stop the same but reducing the exp value, increases the spreading of the GMF, yet the runout length does not change. In the low spreading example in Fig. 3a (exp = 100), the behavior of the downhill flow is restricted to a single flow direction in steeper terrain.
Once the slope flattens out the path diverges with very limited spreading. The small amount of spreading in flatter terrain can be explained by the low Z δ max , which results from the persistence-based routing being dominated by the terrain-based routing.

330
The front of the GMF runout is defined by the runout angle-induced stopping routine with Z δ max = 0 (black). The sides of the GMF process path are defined by the routing flux-induced stopping routine and because Z δ max > 0 the runout angle-induced stopping condition is not met.

Parabolic, channelized slope
This topography has the same extent, center line profile and configuration as the parabolic slope in Sect. 3.1; however, an 335 hour glass shaped channel is added, which begins wide and becomes narrow returning to a wide channel in the runout zone again. The parameterization used for this scenario is α = 25 • , exp = 8, R stop = 3 · 10 −4 and Z δ lim = 8, 849m, such that one can compare it to the simulation results shown in Fig. 3b.
This example highlights the routing flux-induced stopping and the terrain-based routing (Fig. 4). The GMF travels down the channel and does not spread like in the previous example (Fig. 3). That is, the routing algorithms acts on the channelized terrain and concentrates the flux in the center of the channel. The GMF does not spread outside of the channel because the flux that is routed up the channels walls does not exceed the flux cutoff R stop , hence the routing flux-induced stopping criteria is met. Cooler colors indicate areas where the process has a relatively low intensity with regards to Z δ max and warmer colors show areas where the process has a relatively high intensity with regards to Z δ max , which is associated to maximum kinetic energy.. . GMF runout modeled with Flow-Py on a parabolic slope with a channel with α = 25 • , exp = 8, Rstop = 3 · 10 −4 and Z δ lim = 8, 849m. The colors show the value of Z δ max , which is associated to maximum kinetic energy. The topography is a simple parabolic slope connected to a flat plane.

Parabolic, channelized slope with superimposed dam
The topography used in this scenario is the same as in the last Sect. (3.2) including a superimposed obstacle that crosses 345 the terrain such that the GMF must travel uphill to overcome it. We refer to this obstacle as a dam as it could resemble a dam built in the GMF path. This example highlights how the Flow-Py simulation responds to flat or uphill terrain, which is where persistence-based routing will dominate over the terrain-based routing. The parameterization used is α = 25 • , exp = 8, R stop = 3 · 10 −4 and Z δ lim = 8, 849m, so that the result can be directly compared with the spreading example shown in Fig.  3b and the channelized example (Fig. 4). The dam has a shape of a Gaussian function with a width of 75 m and a height of 75 350 m which is added on top of the topography of the parabolic slope with a channel. The center of the dam (maximum height) is located at 1350 m (Fig. 5).

Discussion on model testing and validation
The Flow-Py simulation tool is based on a simple model that allows for regional application and was not specifically designed to model a singular GMF. However, simulations on generic topographies and of single paths provide a visual description of 360 how the implemented routing and stopping routines react to different terrain features and parameterizations. The parameters R stop and exp are primarily responsible for limiting the spreading of the path, where α and Z δ lim are primarily responsible for limiting the runout distance. R stop and exp are dependent on the resolution of the DEM, where α and Z δ lim are not. Figure 6 shows Z δ max values for the center line of all scenarios presented in this Sect., which allows to quantify and validate the model implementation. All simulations yield the same values for Z + Z δ max (where Z is the terrain height) along the center 365 line, although topographies and associated three-dimensional runout extents differ significantly. This is particularly interesting for the third scenario (Fig. 5) where not only Z + Z δ max values are matched but the routing and propagation of the GMF continued beyond the obstacle, where it usually would prohibit any propagation, e.g., with an often employed steepest descent routing approach.
In addition, the model motivation allows to predict the geometrically and theoretically expected solution in terms of runout 370 and z δ . By comparing the geometrically correct solution z δ with the simulation results of Z δ max for each scenario we obtain a match with root mean squared error of 4 · 10 −5 for each simulation result compared with the geometric solution (see Fig. 6).
This in turn validates the discretized model equations and their correct implementation. That is, the cell by cell approach to the solution. This validation however is only relevant for the intensity Z δ max and runout length along the center line. It was not the aim to fully validate the implementation of the spreading algorithm; however, the scenarios show satisfying results where single flow and divergent flow behavior, i.e. ranges from block to fluvial GMF behavior, can be reproduced by changing the Flow-Py parameterization (Fig. 3).
This section is dedicated to highlight the performance of the Flow-Py simulation tool in real terrain and on a regional scale by applying it to the snow avalanche GMF.

Study area description and experimental setup
The study area is located in the mountains surrounding the Austrian villages Vals and Gries am Brenner in Tyrol close to the Italian border (Plörer and Stöhr, 2021). The area of the study area is 104.5 km 2 . The input DEM is freely available from Land 385 Tirol (data.tirol.gv.at) issued under a Creative Commons Attribution 4.0 International (CC BY 4.0) license.
The computation time is dependent on the number of starting cells and the extent of the paths (how divergent/concentrated). We developed an overly simple starting zone model to test the performance of Flow-Py on a regional scale. There are many models for identifying potential avalanche starting zones that use a range of slope inclinations such as 28 • to 60 • Veitinger et al. (2016); Maggioni and Gruber (2003); Pistocchi and Notarnicola (2013) :::::::::::::::::::::::::::::::::::::::::: (Veitinger et al., 2016;Pistocchi and Notarnicola, 2 390 . More information such as terrain curvature, forest cover and average maximum snow depth are used to further restrict the number and size of potential starting zones. The starting zone model employed is based solely on the slope inclinations derived from the 10 m DEM with the goal to provide a sufficient number of starting cells with potentially long runout lengths for performance testing. To achieve this we used two criteria for identifying starting cells: first, starting cells must be located above 1800 m, second, the starting cell must have a slope inclination between 31 • to 34 • . The range of slope inclinations used 395 is much smaller than used in more sophisticated models, this method was used to reduce the number of starting cells with out introducing more information such as forest area, or average snow depth, but rather relying solely on the 10 m DEM. The parameterization for this simulation was α = 25 • and exp = 8 and R stop = 3 · 10 −4 and Z δ lim = 8, 849m, which have successfully been used to model large to very large avalanches . For snow avalanches an exp of 8 on a 10 m resolution DEM has produced good results in the past studies (Huber et al., 2016).

Results
The study area contains 1045311 raster cells 104531100 m 2 ) and starting cells comprise 5.4% of the total study area 56969 raster cells or 5696900 m 2 ), which can be seen in Fig. 7. The simulation took 3 h and 45 min with multi-processing on 16 cores.
Flow-Py identified 642630 cells or 61.5 % of the total study area as part of the avalanche starting, transit and runout zones 405 (see Fig. 8 ). Many of the these cells belong to multiple paths and are, therefore, base cells for many calculations which is

Model customization and adaptability
The third computational experiment highlights the adaptability of the Flow-Py simulation tool with an example for a custom extension that was designed to answer a specific question, but additional information and calculations can be easily added into 445 the Flow-Classes .

Experimental setup and methods
This specific model customization experiment addresses the research question: What areas on the terrain are associated with endangering a location containing infrastructure by a GMF? For this example the parabolic slope from Sect. 3.1 and release raster from Sect. 3 were used as input data with the parameterization α = 25 • , exp = 8, R stop = 3 · 10 −4 and Z δ lim = 8, 849m 450 as well as an additional input raster that contains the location of infrastructure. With the experimental setup defined the initial question can be refined: What raster cells on the synthetic parabolic slope are associated with routing flux of an GMF through a specific set of raster cells that have been identified as locations with infrastructure?
A custom extension called the back-tracking extension was implemented to change the runout model to a model that highlights terrain associated with endangering infrastructure. The back-tracking results will be a spatially explicit subset of the 455 GMF path. Three steps are required to adjust the Flow-Py simulation tool: 1. Load the infrastructure raster as an additional input raster.
2. Adjust the calculation and store new information in the Flow-Class.
3. Save the back-tracking information as a raster and discard the default outputs.
Steps 1 and 3 are simple tasks when using the existing input and outputs of Flow-Py as an example. For the version of Flow-initiate the back-tracking extension and suppress the normal Flow-Py outputs (ee :: see : the Flow-Py repository  for more information on implementing these steps).
Step 2 is the more challenging adjustment that highlights the adaptable nature of the Flow-Class organization. Since the goal of back-tracking is to find the avalanche starting, transit and runout zones that are associated with endangering infrastructure, 465 a new back-tracking variable must be added to the Flow-Class storing information about a cell's parents. If the back-tracking variable of a cell is 0, then this cell is not associated with endangering infrastructure; if it is 1, then the cell is associated with endangering infrastructure. After a path is calculated and before updating the result raster, the back-tracking routine can start.
Starting with the cells identified as a location with infrastructure a family tree can be constructed by looking at which cells acted as parent cells to these infrastructure cells. For each parent cell the back-tracking variable is changed from 0 to 1. After 470 looping over all cells identified as parent cells that are related to cells containing infrastructure, the result raster can be updated with the back-tracking results and the next GMF starting cell can be calculated.
To optimize the back-tracking extension with regards to model run time, the starting cells have been ordered by elevation.
If a starting cell is located in the path of a previously calculated starting cell at a higher elevation, then the lower elevation starting cell is removed from the queue of starting cells that must be calculated. This will greatly reduce the run time of the 475 model as less starting cells and process paths need to be computed; however, no information about the back-tracking is lost because the process path of the lower starting cell will be a subset of the upper starting cell's path, but other output raster such as cell counts (CC) and the Z δ sum are no longer be valid since some starting cells are ignored for optimization. To test the back-tracking extension two types of infrastructure were considered: a linear infrastructure such as a road, railroad or walking path that crosses the terrain, and a single pixel which could represent a building or utility pole.

Figure 9a
shows linear infrastructure (vertical red line) crossing the parabolic terrain, with the areas identified by the backtracking extension. Most of the path located uphill of the infrastructure has been identified by the back-tracking extension except for few cells that lay on its edges. This is because these cells were not parent cells due to the routing flux-induced stopping criterion.

485
The bottom panel of Fig. 9 shows how the back-tracking extension behaves for a single infrastructure cell (red), e.g. a building, in the center of the GMF path. A wedge shape subset of the process path starting at the infrastructure cell and extending up slope is identified by the back-tracking extension, and it is clearly shown that not all uphill cells of the infrastructure cell route flux through the infrastructure cell. In both linear and building infrastructure cases, all cells that lay below the infrastructure cells are not identified by the back-tracking extension.
We used the back-tracking extension to exemplify and highlight the adaptability of the Flow-Py simulation tool. By making small adaptions Flow-Py was changed from a runout model to a model that identifies endangered infrastructure, which 505 demonstrates how Flow-Py can be used to investigate questions related to specific GMFs.

Conclusions and outlook
Flow-Py is an open-source simulation tool for data-based gravitational mass flow (GMF) runout and intensity modeling, which is suitable for spatially explicit applications on a regional scale. GMF is a term that generalizes the flow of various materials in different ratios of solids, water (ice) and air down a slope. The GMF behavior, the runout length and the amount of lateral 510 spreading are all partially dependent on the composition of the material (Pudasaini and Mergili, 2019). Flow-Py handles diverse flow behaviors by providing an adjustable parameterization that acts to control the spreading and runout lengths of the simulated GMF path.
Flow-Py's basic model equations and well-organized solver split the GMFs runout modeling into two routines: 1) routing of the GMF, and 2) stopping of the GMF. The routing routine is further broken down into terrain-based and persistence-based 515 routing, and the stopping routine is further broken down into two stopping criteria based on runout length and the amount of flux. With this, Flow-Py provides and educational GMF model development environment, which combines computational efficiency with low entry barriers for adaptations/extension, such as the presented back-tracking extension.
Besides the local topography two factors influence the spreading of the simulated GMF, namely the exp parameter and the routing flux threshold R stop . However, the four main parameters (runout angle α, divergence control exp, flux cutoff R stop 520 and the limit of the process intensity Z δ lim ) have to be defined based on ones experience or corresponding guidelines to obtain the desired range of motion behaviors corresponding to different materials and their compositions. However, further studies are needed for in depth parameter investigations, including the development of parameter sets which can be used for specific GMF types such as rockfall, different types of snow avalanches or landslides. Part of these parameter studies should include a sensitivity study on the DEM resolution used.

525
The implementation of the model equations to route the flux on the cell level has the advantage that the flow path does not need to be predetermined in contrast to some similar statistical runout methods (Lied and Bakkehøi, 1980). Therefore, Flow-Py combines the simplicity of a runout angle-motivated model with the advantages associated with process-based modeling, providing a corresponding intensity measure and allowing for routing in flat or uphill terrain as we demonstrated in a computational experiment. The results of a second experiment show that the run time of the model is suitable for regional modeling 530 (several 100 km 2 ) (see Sect. 4). The main factor that control model run time is the parameterization, especially the amount of spreading which is controlled by exp and R stop ; however, the number of starting cells and the number of available computer cores are also important factors influencing model run time.
One of the major benefits of Flow-Py compared to existing GMF simulation tools is its well-organized code that allows easy adaptations and extension development. A custom extension was developed for Flow-Py for taking into account terrain 535 complexity with regards to snow avalanches, where automated avalanche terrain exposure scale (ATES) maps were created (Larsen et al., 2020). Future work is being carried out to develop a custom extension which will adapt the stopping criteria to other statistical models (Lied and Bakkehøi, 1980;Barbolini et al., 2011). We presented the back-tracking extension in Sect. 5 to demonstrate adaptability of the simulation tool, which required adjustments to the input data, calculations and output raster.
The additional calculations took advantage of Flow-Py's programming in Python's object-oriented class method, the Flow-

540
Class. The output can be used to identify forests with a direct object protective function by combining the back-tracking results with a map of the current forest cover in a post-processing procedure. More simple extensions have already been developed and used, e.g., the forest extension which has been applied to quantify the forest's protective effects in transit zones of rockfall and starting, transit and runout zones of snow avalanches, adapting the runout angle stopping criteria dependent on forest structure and the intensity (Z δ max ) of the GMF (D' Amboise et al., 2021). 545 We have shown that Flow-Py is an innovative GMF simulation tool that can be applied for basic simulations (e.g., for hazard zone mapping) but also for more sophisticated custom applications such as identifying areas that potentially endanger specific infrastructure. Furthermore, not only presenting Flow-Py in this contribution but also the modeling concepts that motivated its model equations and their implementation in the Flow-Py code enables one to reproduce and understand the basic concepts of GMF modeling and to also use Flow-Py as an educational tool.

550
Code and data availability. The Flow-Py code and user manual can be found on the repository at Neuhauser et al. (2021) Input data and simulation results can be found at