com4FlowPy theory
Note
This section is currently under development/review!
… as such it might still contain some broken links, missing references, typos, etc.
- For additional reference on the concepts and motivation behind
com4FlowPy, we kindly refer to: D’Amboise, C.J.L.; Neuhauser, M.; Teich, M.; Huber, A.; Kofler, A.; Perzl, F.; Fromm, R.; Kleemayr, K.; Fischer, J.-T.(2022). Flow-Py v1.0: a customizable, open-source simulation tool to estimate runout and intensity of gravitational mass flows. Geoscientific Model Development. 15, 2423–2439. (https://doi.org/10.5194/gmd-15-2423-2022)
- Initial information pertaining to the application of
com4FlowPywith consideration of forest-effects can be found in: D’Amboise, C.J.L; Teich, M.; Hormes, A.; Steger, S. and Berger, F. (2021). Modeling Protective Forests for Gravitational Natural Hazards and How It Relates to Risk-Based Decision Support Tools. in: Protective forests as Ecosystem-based solution for Disaster Risk Reduction (ECO-DRR), Teich et al., 2021. (https://www.intechopen.com/chapters/78979)
Background and motivation
Fig. 48 Definition of angles and geometric measures for the calculation of \(Z^{\delta}\), where s is the projected distance along the path and \(z(s)\) the corresponding altitude.
The model equations that determine the run out in three dimensional terrain are mainly motivated by simple, geometric, two dimensional concepts [Hei32, HFKK16, JL11, Korner80, Meissl98] in conjunction with existing algorithms for flow routing in three dimensional terrain [Hol94, HJRZ13, Wic17], controlling the main routing and final stopping of the flow. Fig. 48 summarizes the basic concept of a constant run out angle (\(\alpha\)) with the corresponding geometric relations in two dimensions along a possible process path.
The local travel angle gamma is defined by the altitude difference and projected distance along the path, from the release point to the current location.
The angle delta is the difference between the local travel angle gamma and the runout angle alpha and is related to \(Z_{\delta}\), so when \(Z_{\delta}\) equals zero or \(Z_{\gamma}\) equals \(Z_{\alpha}\), the maximum runout distance is reached.
\(Z_{\alpha}\) can be interpreted as dissipation energy.
\(Z_{\gamma}\) is the altitude difference between the starting point and the current calculation step at the projected distance \(s\).
\(Z_{\delta}\) is the difference between \(Z_{\gamma}\) and \(Z_{\alpha}\), so when \(Z_{\delta}\) is lower or equal zero the stopping criterion is met and the flow stops. \(Z_{\delta}\) is associated to the process magnitude and can be interpreted as the kinetic energy or velocity of the process.
The major drawback of implementing the geometric runout angle concepts is that they require a predefined flow path in two dimensional terrain. To allow for an enhanced routing in three dimensional terrain without prior knowledge of the flow path we combine these concepts [HFKK16] with extensions of existing algorithms [Hol94, HJRZ13, Wic17] that are described in the following sections.
Spatial input and iterative calculation steps on the path:
In nature a GMF has one or more release areas that span over single or multiple release cells. Flow-Py computes the so called path, which is defined as the spatial extent of the routing from each release cell. Each release area (single raster cell in release area layer) has it’s own unique path (collection of raster cells), and a location on the terrain (a single raster cell) can belong to many paths. Flow-Py identifies the path with spatial iterations starting with a release area raster cell and only iterating over cells which receive routing flux. The corresponding functions are implemented in the code in the flow_class.calc_distribution() function.
To route on the surface of the three dimensional terrain, operating on a quadrilateral grid, we implement the geometric concepts that have been sketched in the model motivation utilizing the following cell definitions:
Fig. 49 Definition of parent, base, child and neighbors, as well as the indexing around the base.
Each path calculation starts with a release cell and operates on the raster, requiring the definition of parent, base, child and neighbor cells (see Fig. 49). The base cell is the cell being calculated on the current spatial iteration step. The 8 raster cells surrounding the base cell are called neighbor cells (n, i) which have the potential to be parents (supplying flux to base cell), or a child (receive flux from the base cell). In 2d the base cell corresponds to the cell/location at the distance s along the path in Fig. 48.
Every base has at least one parent cell, except in the first calculation step from the release cell, where we start our calculation, this would be at \(s = s_0\) in Fig. 48.
During an iteration step a raster cell from the iteration list is identified as the current base cell. The routing flux is calculated across the base cell from the parent cell to possible child cells. The goal is to keep the spatial iteration steps to a minimum, which is achieved by only adding neighbor cells to the iteration list that have flux routed to them from the base cell and do not meet either of the stopping conditions. These cells are called child cells. Child cells that are not already on the iteration list are added to the list and flow_class python object is created for the raster cell. The child cells flow_class has the parent added to it as a source for routing flux. By being added to the iteration list the cell has been recognized as being part of the GMF path and will be the base cell for a future iteration step.
When the iteration list is empty and all potential children fulfill one of the stopping criteria:
\(Z^{\delta}\) has to be smaller than zero: \(Z^{\delta}<0\)
Routing Flux has to be smaller than the flux cut off: \(R_i < R_{Stop}\)
the path calculation is finished. The required information is saved from the cell class to the summarizing output raster files. Then the calculation starts again for the next release cell and respective flow path. The spatial extent and magnitude for all release cells are summarized in the output raster files, which represent the overlay of all paths.
Every path is independent from the other, but depending on the information we want to extract, we save the highest values (e.g. \(Z_{\delta}\)) or sums (e.g.Cell Counts) of different paths to the output raster file.
\(Z^{\delta}\)
For each base cell in a path we solve the equations (Eq.38, Eq.39, Eq.40) for every neighbor n, if \(Z_{bn}^{\delta}\) is higher than zero, this neighbor is defined as a potential child of this base, and routing in this direction is possible.
Here \(S_{bn}\) is the projected distance between the base and the neighbor.
As \(Z_{bn}^{\delta}\) can be interpreted as process magnitude (and kinetic energy or velocity respectively) it is possible to limit this value to a maximum. In comparison to process based modeling approaches this would correspond to maximum velocity induced by a velocity dependent turbulent friction term.
The projected path lengths, or total travel distance to one of the neighbors (\(S_n\)) equals the path length to the base (\(S_b\)) plus the path from base to the neighbor (\(S_{bn}\)), which reads:
As there are many possibilities for the path from the starting point to the actual cell or base, the shortest path is taken into account, corresponding to the highest \(Z^{\delta}\) in the base. If \(Z^{\delta}_{max}\) is set to infinity, or as in the code to 8848 m (= Mount Everest), we can calculate the shortest path from the starting point to the base and yields the total projected travel distance:
This equations determine the routing and corresponding run out distance for the process, the next steps demonstrate how spreading is handled on the surface of the three dimensional terrain.
Persistence based routing:
The persistence contribution \(P_i\) aims to capture inertia related effects, and takes the change in flow direction into account [3]. The direction contribution is scaled with the process magnitude \(Z^{\delta}_{parent}\), such that the direction from a parent cell with higher process magnitude has more effect on the path routing and direction.
The direction contributions \(D_n\) are defined by the cosine of the angle between parent, base and child/neighbor minus \(\pi\).
Therefore the direction contribution limits the maximum number of potential children to three, getting input via the persistence function from one parent.
In the first calculation step, at the release or start cell no parent cells are defined and the persistence is set to one. So the first calculation step is solely determined by the terrain contribution.
Terrain based routing
The terrain based routing is solely dependent on the slope angle phi. The exponent exp allows to control the divergence of the spreading. The Holmgren (1994) algorithm [Hol94] is used in different kind of models and works well for avalanches but also rockfall or soil slides. For avalanches an exponent of 8 shows good results. To reach a single flow in step terrain (rockfall, soil slides, steepest descend), an exponent of 75 is considered.
To overcome the challenge of routing in flat or uphill terrain, we adapted the slope angle \(\Phi\) for the normalized terrain contribution to:
Routing flux:
The routing flux summarizes the persistence and terrain contributions according to Eq.48:
where i is the direction and n are the neighbors from 1 to 8. \(R_i\) is then the routing flux in direction \(i\). \(R_b\) is the flux in the base, for a release cell or starting cell the flux of the base equals one. The result of Eq.48 is a 3 x 3 array with assigned flux values. A normalization stage is then required to bring the sum of the \(R_i\)’s to the value of \(R_b\). This aims at avoiding loss of flux.
Flow chart / overview:
In Fig. 50 the algorithm of the computational implementation is sketched, including function and files names with respect to the code in the repository.
The file main.py handles the input for the computation and splits the release layer in tiles and saves them in a release list. Then the main.py starts one process per tile, which calls the flow_core.py and starts the calculation for one release cell and the corresponding path. The number of processes is depending on the hardware setting (CPU and RAM). Whenever a new cell is created flow_core.py calls flow_class.py and makes a new instance of this class, which is saved in the path. When the calculation in flow_core.py is finished it returns the path to main.py which saves the result to the output rasters.
Fig. 50 Flow chart of the Flow-Py computational process and an overview of the files and what they manage.
Acknowledgments
This study was carried out in the framework of the GreenRisk4Alps project ASP635, funded by the European Regional Development Fund through the Interreg Alpine Space programme. Additional financial support from the AvaRange (www.AvaRange.org, international cooperation project “AvaRange - Particle Tracking in Snow Avalanches” supported by the German Research Foundation (DFG) and the Austrian Science Fund (FWF, project number I 4274-N29) and the AvaFrame (www.AvaFrame.org, AvaFrame - The open Avalanche Framework is a cooperation between the Austrian Research Centre for Forests (Bundesforschungszentrum für Wald; BFW) and Austrian Avalanche and Torrent Service (Wildbach- und Lawinenverbauung; WLV) in conjunction with the Federal Ministry Republic of Austria: Agriculture, Regions and Tourism (BMLRT)) projects are greatly acknowledged.
Citation
If you are using com4FlowPy in your work, we are happy if you include the following reference in any resulting
publications:
D’Amboise, C.J.L.; Neuhauser, M.; Teich, M.; Huber, A.; Kofler, A.; Perzl, F.; Fromm, R.; Kleemayr, K.; Fischer, J.-T. (2022). Flow-Py v1.0: a customizable, open-source simulation tool to estimate runout and intensity of gravitational mass flows. Geoscientific Model Development. 15, 2423–2439. (https://doi.org/10.5194/gmd-15-2423-2022)
com4FlowPy applications
com4FlowPy has been applied in different studies in recent years. Here are the ones we know about:
autoATES
Toft, H.B.; Sykes, J.; Schauer, A.; Hendrikx, J. and Hetland, A. (2024). AutoATES v2.0: Automated Avalanche Terrain Exposure Scale mapping. Nat. Hazards Earth Syst. Sci., 24, 1779–1793 (https://doi.org/10.5194/nhess-24-1779-2024)
Sykes, J.; Toft, H.B.; Haegeli, P. and Statham, G. (2024). Automated Avalanche Terrain Exposure Scale (ATES) mapping – local validation and optimization in western Canada. Nat. Hazards Earth Syst. Sci., 24, 947–971 (https://doi.org/10.5194/nhess-24-947-2024)
Spannring, P.; Hesselbach, C.; Mitterer, C. and Fischer, J.-T. (2024). Classification of avalanche terrain: a open-source model chain for the Avalanche Terrain Exposure Scale. Proceedings Interpraevent 2024, Vienna, Austria (https://www.interpraevent.at/en/proceeding/proceedings-ip-2024)
von Avis, C.D.; Sykes, J. and Tutt, B. (2023). Development of large scale automated Avalanche Terrain Exposure Scale (ATES) ratings in collaboration with local avalanche experts. Proceedings, International Snow Science Workshop, Bend, OR, USA (http://arc.lib.montana.edu/snow-science/item/2998)
Huber, A.; Hesselbach, C.; Oesterle, F.; Neuhauser, M.; Adams, M.; Plörer, M.; Stephan, L.; Toft, H.B.; Sykes, J.; Mitterer, C. and Fischer, J.-T. (2023). AutoATES Austria - Testing and application of an automated model-chain for avalanche terrain classification in the Austrian Alps. Proceedings, International Snow Science Workshop, Bend, OR, USA (http://arc.lib.montana.edu/snow-science/item/2989)
Hesselbach, C. (2023). Adaptaion and application of an automated Avalanche Terrain Classification in Austria. Masters’ Thesis. University of Life Sciences (BOKU), Vienna (https://forschung.boku.ac.at/de/publications/175549)
Schumacher, J.; Toft, H.B.; McLean, J.P.; Hauglin, M.; Astrup, R. and Breidenbach, J. (2022). The utility of forest attribute maps for automated Avalanche Terrain Exposure Scale (ATES) modelling. Scandinavian Journal of Forest Research, 37:4, 264-275 (https://doi.org/10.1080/02827581.2022.2096921)
other
Perzl, F.; Huber, A.; Fromm, R. and Teich, M. (2024). Estimation of potential snow avalanche hazard probability in areas below protective forests in Austria. Proceedings Interpraevent 2024, Vienna, Austria (https://www.interpraevent.at/en/proceeding/proceedings-ip-2024)
D’Amboise, C.J.L; Teich, M.; Hormes, A.; Steger, S. and Berger, F. (2021). Modeling Protective Forests for Gravitational Natural Hazards and How It Relates to Risk-Based Decision Support Tools. in: Protective forests as Ecosystem-based solution for Disaster Risk Reduction (ECO-DRR), Teich et al., 2021. (https://www.intechopen.com/chapters/78979)