Modeling Protective Forests for Gravitational Natural Hazards and How It Relates to Risk-Based Decision Support Tools

Abstract

Simulation tools and their integrated models are widely used to estimate potential starting, transit and runout zones of gravitational natural hazards such as rockfall, snow avalanches and landslides (i.e., gravitational mass flows [GMFs]). Forests growing in areas susceptible to GMFs can influence their release and propagation probabilities (i.e., frequency and magnitude of an event) as well as their intensity. If and how well depends on the GMF type, the topography of the terrain and the forest’s structure. In this chapter, we introduce basic concepts of computer models and state-of-the-art methods for modeling forest interactions with rockfall, snow avalanches and landslides. Furthermore, an example of a protective forest routine embedded in the runout angle-based GMF simulation tool Flow-Py will be presented together with its parameterization for forest-GMF interactions. We applied Flow-Py and two custom extensions to model where forests protect people and assets against GMFs (the protective function) and how forests reduce their frequency, magnitude and/or intensity (the protective effect). The goal of this chapter is to describe protective forest models, so that practitioners and decision makers can better utilize them and their results as decision support tools for risk-based protective forest and ecosystem-based integrated risk management of natural hazards.

1. Introduction

Simulation tools and their integrated models are widely used by the scientific community and practitioners for their predicting power, which is where science and practice overlap. One major advantage of simulation tools is their potential to highlight what is known about a system and where knowledge gaps exist. Every model is a simplification of reality and practitioners can use the output of models as a decision support tool (see chapter [1] of this book). By simplifying a natural system, it is necessary to make assumptions, which goes hand in hand with reducing complexity and loss of detail. To best benefit from using simulation tools, practitioners should therefore be familiar with the basics of the underlying models to understand their limitations.

In this chapter, we discuss how to incorporate forests into gravitational mass flow (GMF; [2]) models to gain an understanding of where forests protect people and assets (their protective function) and to estimate how forests reduce the frequency, magnitude and intensity of gravitational natural hazard (their protective effect; see chapters [3, 4] for definitions and details on forest-GMF interactions). We first introduce basic concepts of computer models regarding protective forests and then summarize some of the state-of-the-art methods used for modeling forest interactions with rockfall, snow avalanches and landslides in their starting zone as well as transit and runout zones. Lastly, an example of a protective forest routine embedded into Flow-Py, a GMF simulation tool based on a runout angle (also referred to as the travel or α-angle) model [2, 5], is presented together with the justifications for the parameterization and implementation. This example highlights the development process of simulation tools and allows the user a deeper look into the assumptions that are necessary for modeling protective forests.

The model development process can be thought of as a 4-part cycle as shown in Figure 1. The cycle starts with obtaining a greater process understanding via laboratory experiments, field measurements, remote sensing data acquisition or existing observations. The next step of the cycle is to incorporate that new knowledge into an existing model or newly developed model. The third step is a model validation where the model results are compared to observations. The last step is a model evaluation where the validation results are used to highlight improvements that have been made during the cycle and, more importantly, to reveal remaining knowledge gaps. These knowledge gaps are opportunities for new research to obtain a deeper process understanding and the cycle repeats.

Figure 1.

The model development cycle.

In general, models can be grouped in several categories. For simplicity, we group GMF models by two main characteristics:

1. with regards to the size of the modeling domain which is linked to a model’s (spatial) resolution, and

2. the general methods and techniques that were used to develop a model.

The domain sizes and model types that will be discussed are the regional (10s to 100s of km²) and the hill slope or path scale (less than 10 km²), and data-driven statistical models (hereafter referred to as statistical models, e.g., the α-β model [6], Flow-R [7], or the topographic approach of [8]) and process-based physical models (hereafter referred to as physical models, e.g., RAMMS [9] or Rockyfor3D [10]).

We use the expression “statistical model” to summarize data-driven, data-based and machine learning approaches whose quantitative outcome originates directly from empirical data without needing a deeper understanding of the underlying processes. The equation that results is based on statistics and the parameters have no physical meaning but are inferred from real-world observations. The major strengths of statistical models are the simplistic parameterization and their flexibility in terms of input data requirements compared to physical models. However, statistical models can only answer questions they were designed for, which strongly depends on the input data used to develop and calibrate them. The major challenge for developing a reliable statistical model is, therefore, the collection of data on which to base the statistical relationships and parameterizations on.

In contrast, physical models break the main modeling question into the governing processes or smaller sub-processes that are expressed in equations. By compiling the sub-processes, information of the system is gained, and the modeling question can be addressed. Physical models usually involve more calculations and require larger computational resources than statistical models. However, the major advantage of physical models is that the variables, parameters, and some intermediate calculations (usually) have physical meanings providing additional information, e.g., the energy of the GMF or the depth of the flow. This strength of physical models is also one of their major drawbacks, because they require very accurate input data and parameterization. That is, physical models can have large parameter sets, whose values can be difficult to measure and are often impractical to obtain for many modeling efforts. Moreover, although physical models split the main modeling question into sub-processes, significant feedback between processes still exists and cannot be ignored, i.e., the results of one sub-process are fed into another process. In practice, the calibration of physical models is often based on empirical parameter adjustments that gives these models a statistical character. In the end, a model can only be as “good” as the data that is used for its development, which is true for any (statistical or physical) model development.

GMF and protective forest models have often elements of both statistical and physical models (e.g. [11, 12]). The sub-processes that are lacking a strong understanding or that demand an excessive amount of data, such as forest interactions with a GMF, are frequently represented with statistical models. GMF runout models identify the spatial extent or how far a hazard reaches down slope. To account for the vegetation effect in such a model, it requires the adjustment to the mechanisms for modeling the two main effects on gravitational natural hazards (see chapter [4] of this book):

1. the starting zone identification (release susceptibility or probability) in forested areas, and

2. the runout model (the movement of the GMF) in forested locations.

The mechanisms to adapt a GMF model to a protective forest model are different depending on the type of model (statistical or physical), the size of the modeling domain and the research question being investigated. In statistical models, the mechanisms to account for interactions between a GMF and forest have to be rooted in empirical data, while these mechanisms must adjust the equations of physical models. However, the size of the modeling domain often determines the type of model that is applied. For large areas or regional scale modeling, statistical models are regularly used due to the lack of existing and detailed input data. On single paths or the hill slope scale precise input data and parameterizations can often be collected more easily. Compared to regional-scale studies, investigations at the hill slope scale provide finer details relevant to individual properties in terms of the resolution and accuracy of simulated outputs (see chapter [13] of this book). In contrast, detailed information is often less useful for studies at the regional scale, and input data and simulation results are often simplified and presented as summaries of higher resolution data (e.g., average values, standard deviations or trends). As a practitioner it is important to choose the appropriate model depending on the question at hand.

2. Protective forest models

Protective forest models can be applied to model (1) forests’ protective functions, and/or (2) forests’ protective effects. That is, they can be used to identify locations of forests with an (object) protective function (see chapter [3] of this book for definitions), and/or to quantify the degree to which these forests protect a location.

Protective forest models are used by different user groups to support decision making. For example, the road administration could use a protective forest model to investigate the degree of protection a forest provides in a single avalanche starting zone or avalanche flow path that endangers a section of a road. Or a local or state government could use protective forest models to investigate the extent of protective forest located in a region (e.g., municipality, state, or province) to determine a budget for managing these forests. These examples require very different amounts of information and different mechanisms that must be implemented into the GMF model, so that it can be applied as a protective forest model. In this subsection the main mechanisms currently used to include protective forest into different types of GMF models are introduced.

2.1 Identification of protective forests’ function

To address the question where object protective forest is located, a union between the spatial distribution of the forest and the spatial distribution of the GMF areas (consisting of starting, transit and runout zones) that endanger assets must be established. Therefore, the applied GMF runout model must be able to discriminate between parts of the GMF that endanger infrastructure and parts that do not, because not all parts of the transit zone and/or runout zone will necessarily reach infrastructure for a given starting zone. This adaption to GMF runout models is relatively straight forward. The union between the spatial distribution of forest and the sub-set of GMF areas that endanger infrastructure is done in a post-processing step. The resulting maps highlight forests that have a direct object protective function, which has been mapped, for example, in Switzerland (project: SilvaProtect-CH [14, 15]) and in Austria (projects: PROFUNmap for data integration [16]; GRAVIMOD II for shallow-seated landslides [17, 18]; DAKUMO for snow avalanches [19]; GRAVIMOD I for rockfall and snow avalanches [20, 21]).

2.2 Quantification of protective forests’ effect

Two interactions between forests and GMFs must be considered to answer how much a forest reduces the frequency, magnitude and/or intensity of a GMF (see section 1). That is, the initial release of mass can be influenced by forest and the forest can also interact with the movement in the GMF’s transit and runout zones dissipating energy by mass reduction or increasing friction, which differs depending on the GMF type.

The process for releasing a GMF is highly dependent on the material type and composition. Therefore, each GMF type requires different models for identifying its starting zones and/or quantifying associated release and propagation probabilities. Deeper process understanding of how the forest interacts with the movement and release mechanics is available for some types of GMFs, while for others the process understanding is weak (see chapter [4] of this book). A lack of process understanding will result in less precise model parameterizations for identifying starting zones and simulating forest-GMF interactions. Therefore, hazard-specific parameterizations and mechanisms are needed for the different GMF types such as rockfall, snow avalanches or landslides.

3. State-of-the-art forest-GMF interaction modeling approaches

3.1 Rockfall

Two main categories of rockfall models exist: 1) models that identify and characterize the potential starting (i.e., release or disposition) zones of boulders, and 2) models that simulate rockfall trajectories. That is, the effects of forest must be considered separately for the starting zone, and transit and runout zones, since forest might even increase rockfall activity in release areas due to wedge effects of the roots, and snow or wind effects on stem movements [22]. In the transit and runout zones, forest can shorten the runout lengths of boulders if the forest on the slope has a minimum length and no large gaps. For example, a minimum forested slope length of 250 m with openings below 40 m in length is one target [23] to potentially provide protection against rockfall. Currently, no rockfall model exists that combines disposition and trajectory modeling [24].

3.1.1 Rockfall release models

The impacts of trees in starting zones are mainly destabilizing due to the blasting and leverage effects of root systems and/or windthrow, which have not been considered in release models. The choice of modeling technique to identify rockfall starting zones depends on data availability and the desired resolution of the output. Physical rockfall release models consider the internal friction of bedrock, bedrock types, slope inclination, foliation and fractures [25]. The input data requirements to physically model the rockfall release mechanisms limit the application of this type of models to one or few starting zones on a hill slope scale. In contrast, statistical rockfall starting zone models have been developed using digital elevation models (DEMs), which are often the only reliable data source for regional-scale modeling. This type of models applies a threshold based on the slope calculated from a DEM and may consider some local geology [26]. The slope threshold can be applied to DEMs of different resolution by adjusting it accordingly [27].

3.1.2 Statistical rockfall propagation models

Released blocks move mainly by sliding, rolling, jumping or bouncing downslope; however, knowing the exact type of movement is not necessary for statistical rockfall models, which use a runout-angle approach to predict rockfall runout lengths [28]. Two angles are used for rockfall modeling: first, the classical runout angle (also travel or α-angle), which is the angle from the top of the starting zone to the furthest runout of a block. This runout angle assumes that the starting zone (release or source area) is known. The second angle is the shadow angle, which is the angle of the line that is drawn from the bottom of a cliff face (or the lowest possible point of a rockfall release) to the furthest reach of the block runout [22]. The shadow angle describes the maximum travel distance of blocks by intersecting the topography with an energy line starting at the base of the rock face. Empirical studies on α-angles of rockfall trajectories where boulders can still bounce, roll or slide suggest a range between 51.2° and 28.5° [29, 30, 31]. A regression approach was used to specify α-angles for rockfall modeling by means of observed α-angles with an optimal solution of 32.9° [29], indicating that values between 30° and 35° provide useful α-angles to model large rockfall distances.

Individual trees can dissipate the energy of a falling rock by the impact of the boulder, deformation of the stem, rotation or translation of the root, or rebound of the boulder [32, 33]. Therefore, in addition to the length of the forested slope, the protective effect of a forest stand depends strongly on its structural properties, mainly average stem diameter, stem density and/or basal area (see chapter [4] of this book), and has been incorporated into both statistical and physical rockfall trajectory models. That is, forests with a high basal area and high stem density have been identified as a particularly effective measure against rockfall by reducing the kinetic energy and velocity of falling and bouncing blocks in the transit zone [34], which translates to an average increase of the runout angle by 6° [35, 36]. For coherent modeling results, calibrating statistical rockfall models with observations of past events is key; however, physical rockfall models that account for the energy dissipating effect of forests are also being used when such data is not available.

3.1.3 Rockfall trajectory models at the slope scale and integration of forest effects

Physical rockfall trajectory models calculate the runout length and trajectories of blocks and may consider the block’s shape and type of movement as well as its interactions with different underlying surfaces (depending on roughness and soil cover) and the vegetation. Many physical rockfall models provide statistical distributions of the block jump height, velocity and kinetic energy at each point along a slope; however, only few also distinguish the block movement types. Current research is dedicated to also integrating block fragmentation during the fall process, which has not yet been implemented in rockfall models [37, 38].

Physical rockfall trajectory models have initially been developed in a two-dimensional framework. Similar to statistical models, forest-block interactions are implemented by a kinetic energy dissipation of a block in forested terrain, which is here translated by an equivalent friction coefficient [39]. Two-dimensional physical rockfall trajectory models are often employed in simulation experiments for parameterizing statistical rockfall models to quantify the effects a forest has on the dispersion of possible block trajectories. However, these models are not able to represent the variability in the spatial distribution of trees, tree species and stem diameters, and how they affect rockfall trajectories.

Spatially explicit simulations of forest effects on rockfall require three-dimensional models that account for the effects of individual trees on kinetic energy dissipation, lateral rebound and impact [40]. Therefore, three-dimensional trajectory models have been developed to simulate block dispersion depending on the topography of the terrain rather than just on a set of several two-dimensional trajectories [25, 41, 42, 43]; however, only few consider the protective effect of forest, e.g., Rockyfor3D by a reduction in kinetic energy dependent on stem diameter and the percentage of coniferous trees [44]. In RAMMS::ROCKFALL, non-smooth mechanics are coupled with hard contact laws in the framework of a discrete element model; the plugin to simulate the effect of forest on rockfall trajectories will be released in fall of 2021 [45].

3.2 Snow avalanches

The primary effect of forest on slab avalanche is to reduce the probability of release by stabilizing the snowpack on the ground [46] (see also chapter [4] of this book). A secondary effect of forest is its capability to stop or significantly decelerate small-to-medium avalanches starting within forests or close to the upper timberline [47]. However, breakage, uprooting and overturning of trees also reduce runout lengths of medium-to-large avalanches starting above the timberline [48, 49, 50, 51]. Both protective effects are influenced by forest structure in terms of canopy cover, stem density, species composition and size and distribution of forest gaps [52, 53]. To model these two effects usually requires different approaches and even different models.

The models used to identify an avalanche starting zone or to quantify the probability of an avalanche release in forests are primarily statistically based and can be applied to regional or hill slope scales [54, 55]. Models that quantify the forest effect on avalanche dynamics are mainly physical models for regional or hill slope scales [56]; however, the processes that describe forest-avalanche interactions can be expressed physically or statistically depending on the resolution of the avalanche dynamics model. Simpler statistical models are used more often on regional scales [57], while physical models are often limited by the required detail of the input data and are more applicable on the hill slope scale [58].

3.2.1 Snow avalanche starting zone and release models

The most basic release area models locate potential avalanche starting zones with statistical relations to terrain features [59, 60, 61]. In these models, forest is often oversimplified with no consideration of its type or structure. More sophisticated models, which usually require additional information on the snow climate, terrain roughness and/or vegetation type, can also quantify the avalanche release probability [62, 63]. The integration of forest in these models is necessary to quantify the forest effect in dependence on its structure, so that it can be integrated into risk assessments [64].

3.2.2 Statistical snow avalanche-forest interaction models

Statistical avalanche models are widely used to predict maximum runout lengths of extreme avalanches. These models are developed from topographic parameters of paths and observed runout lengths of representative avalanche events. The two mainly applied statistical models are the α–β model [6], and the runout ratio model [65], which uses a runout ratio or density probability function to fit a distribution to observed runouts. However, modeling forest-avalanche interactions with statistical models is challenging since (1) the existing approaches are designed to predict extreme avalanches where forest has a very limited effect on runout [47], and (2) reliable avalanche observations in forested terrain, which are needed for model parameterization and calibration, are rare (see also chapter [66] of this book). For example a terrain and runout ratio analyses for 45 forest-penetrating avalanches, which showed that all but one avalanches stopped up-slope of or at the β-point (the point at which the slope first becomes 10°) [49]. A plot of runout ratio probabilities was derived that can be used in the field to determine the likelihood of an avalanche travelling through forest to a given point on a slope.

In section 4, we present another example of a protective forest routine that was embedded in the GMF simulation tool Flow-Py [5], but is also applied in the online “Protective Forest Assessment Tool – FAT” [67, 68] – a risk-based decision support tool to estimate the value forest has for protecting buildings and infrastructure against GMFs and to compare it to technical and avoidance measures (see also chapter [1] of this book).

3.2.3 Physical snow avalanche dynamics-forest interaction models

Early approaches to model forest-avalanche interactions apply an increase in friction in forested areas to a Voellmy-type relation accounting for the decelerating effect of forest on the avalanche flow and thus a reduction in runout length [69, 70, 71]. Voellmy-type models split the total friction into a velocity-independent Coulomb friction term and a velocity-dependent “viscous” or “turbulent” friction [72]. The Coulomb friction is thought to summarize snow properties, whereas the velocity-dependent turbulent friction term expresses the topography and roughness of an avalanche flow path [73, 74]. To model the braking effect of forests on avalanches the turbulent drag of the basal friction is increased in forested areas compared to open unforested terrain [50, 55]. This increase in friction is supposedly caused by a combination of different forest-avalanche interactions such as breaking, overturning, uprooting or entrainment, which are thought to act mostly on the velocity-dependent turbulent friction [48]. However, this interpretation of the Voellmy friction terms is not ideal because it is based on expert judgment rather than on measurements [12]. Recently, a small-scale experiments of granular flows traveling through regularly spaced ‘trees’ was used to show that the overall deceleration rate can be predicted by applying a stem density-dependent effective basal friction coefficient [75]. The friction approach has been tested for large-scale fast-moving avalanches where the braking effects are small and occur over longer runout lengths [51, 76]; however, this method may not be valid for small-to-medium avalanches [77]. That is, snow detrainment which is the main process of forest-avalanche interactions in small-to-medium avalanches is not well represented with a frictional relationship and the local braking effect of forests on avalanche flow is difficult to model at the grid scale [56].

Based on field observations, Feistl et al. [56] developed an additional one-parameter function (detrainment function) to include forest-avalanche interactions in physical avalanche dynamics models. This function accounts for the snow that is deposited behind trees, groups of trees or remnant stumps by a combination of impact, rubbing dissipation, deflection, cohesion and jamming. The stopped mass is extracted from the avalanche volume and the corresponding momentum is removed from the total momentum of the moving snow leading to a reduction in runout length. The braking effect of forests on avalanche flow and, therefore, the mass to be removed from the avalanche volume is summarized in one parameter (the forest detrainment coefficient K) representing different forest structures [77].

3.3 Landslides

A multitude of possibilities exist to build spatial models for landslides, which are considered crucial for implementing efficient risk reduction strategies and for managing landslide risk [78, 79, 80]. Spatial landslide models are regularly applied to provide estimates on landslide-prone terrain, slope (in)stability or landslide hazard and are employed at different scales to obtain insights into the predisposing, preparatory or triggering factors of slope instability [81, 82, 83].

Spatial landslide models are generally based on qualitative (expert judgment/heuristic) or quantitative (statistically or physically based) approaches [78, 84]. The results of heuristic procedures can be considered subjective as they are primarily based on expert-driven weighting schemes. Instead, statistical and physical models provide numerical outputs and are extensively applied within current scientific studies. The wealth of currently available landslide modeling tools has certainly facilitated the spatial analysis of landslide processes (e.g., [7, 85, 86, 87, 88, 89]). However, despite the recent technical advancements, the explanatory power and applicability of a quantitative spatial landslide model still relies on a non-trivial balance between the quality of the available input data, the model complexity and the envisaged spatial coverage [78, 90, 91]. The following two subsections focus particularly on quantitative spatial landslide models for starting zones of shallow slope (in)stabilities and the implementation of forest’s protective effects.

3.3.1 Statistical models for landslide release: regional scale

Statistical landslide models are usually based on a classification algorithm or are regression-based models that link landslide inventory data (e.g., landslide presence/absence information, counts on landslides) to a variety of environmental variables that are supposed to represent the static or dynamic causes of slope instability. The resulting relationships can then be transferred to each spatial unit of a study site (e.g., raster cell or catchment) to derive a static landslide susceptibility map or spatio-temporal landslide predictions [92, 93]. Statistical models are mainly applied for larger areas, as they are not reliant on a specific set of subsurface parameters that are hard to measure (e.g., soil properties), and because their predictive power increases with the number of observations [94, 95]. Several studies have shown that statistical models can provide insights into the contribution of different forest types, land cover changes and timber harvesting on landslide occurrence [81, 96, 97]; however, recent research conducted within the framework of the project GreenRisk4ALPs [98] also highlighted that an interpretation of statistical models must be done with great care [92], since widespread errors in the underlying input data may lead to distortions in the modeling results and a wrong inference on the causes of slope stability. For instance, a model based on land cover variables and landslide data that is incomplete in forested terrain is likely to underestimate the true stabilizing effect of the forest cover [99].

3.3.2 Physical models for landslide release: catchment scale

Spatial physical slope stability models are based on physical laws and formally allow to analyze causes and effects. Such models are frequently entitled the most advanced approaches to spatially assess slope instability [90]. Physical slope stability models calculate the ratio between stabilizing and destabilizing forces (i.e., factor of safety or probability of failure) for a planar sliding plane [100, 101] or non-planar shear surface [89], while the integration of dynamic components (e.g. hydrological modules) allows to trace the stability state over time. In theory, the (de)stabilizing effects of trees and/or forest cover on slope stability can be considered within such physical models by accounting for hydrological processes (e.g., evapotranspiration, canopy interception) and mechanical forces (e.g., cohesion, surcharge) [82, 102, 103, 104, 105].

The high potential of physical models to elaborate causes of slope instability is seldomly questioned while recent technical advancements even allow to run these models efficiently for larger study areas [85]. In practice and especially when evaluating the effects of vegetation on slope stability, the actual explanatory power of such analyses is restricted by a limited availability of geotechnical data and an appropriate spatio-temporal description of surface and subsurface biomass-related parameters (e.g., root system description, surcharge, evaporation) [82, 95, 106]. For instance, a detailed elaboration of the mechanical effects of roots (i.e., root reinforcement) requires additional area-wide information on potential slip surface depths, because anchoring becomes particularly relevant as soon as the roots penetrate the sliding plane. At the same time, sliding surfaces might not be equally distributed in space and be co-determined by prevalent soil depths and topography (e.g., slope angles). Further spatial information on species distribution, tree age and soil properties might be necessary to approximate root distribution (penetration depth) while geotechnical soil parameters are also known to influence the cohesion forces of a plant. Ultimately, this example also highlights that sophisticated physical models are a simplification of reality and challenging to parameterize, even for smaller and data-rich study sites [78, 84].

4. Protective forest modeling with the Flow-Py simulation tool

The Flow-Py simulation tool contains a GMF runout model for regional scales that allows users to customize GMF simulations by adjusting the parameterization or developing extensions to adapt the calculations, input data and outputs of the model [2, 5]. Flow-Py integrates a statistical model where the runout is solved by splitting the modeling question into two sub-questions, which are addressed in the two modeling routines:

1. The stopping routine: is a stopping criterion met?

2. The routing routine: if not, to where does the flux of the GMF (a portion of the mass flow) gets distributed?

The solver operates on a DEM by calculating the flux from one raster cell to the next until a stopping criterion is met. The cell-by-cell calculation results in the runout (magnitude) and intensity of the GMF without the need to predefine a flow path.

GMFs can flow very differently depending on the material of the mass [107]. Flow-Py overcomes this challenge by adapting the parameterization to fit the type of GMF being modeled using only four parameters: (1) runout angle, (2) divergence exponent, (3) flux cutoff, and (4) maximum kinetic energy height. The runout angel needs to be derived from observations and describes the average reach of the GMF measured from the top of the starting zone to the end of the furthest runout (Figure 2). The divergence exponent adjusts the concentration and spreading of flux across the terrain. The flux cutoff is a limit to describe the amount of flux needed to further propagate the GMF. Lastly, the maximum kinetic energy height limits the kinetic energy of the GMF, which is essentially a limit of the maximum velocity. Two criteria stop the flux, if (1) the simulated GMF runout reaches further than the predefined runout angle, and/or (2) the amount of mass, which is spreading across a slope, reaches a critical value that is required for further propagation. The routing routine calculates to which cell(s) the flux is distributed and uses the DEM to favor directions with steeper downward slopes, as well as persistence, which is similar to momentum but derived by statistical means and does not consider the GMF’s mass (for further information and details see [2]).

Figure 2.

The concept behind the forest extension. The solid brown lines indicate the runout length with α being the runout angle for a specific GMF. The dashed lines show the runout length with the forest being considered and α_forest, the runout angle in forested areas. Artwork: Karl Kleemayr.

To address specific questions regarding protective forest, we developed two custom extensions for Flow-Py to:

1. Identify locations of forests with an object protective function (back-tracking extension), and

2. Quantify their protective effect (forest extension).

From a modeling perspective these are two very different questions that generally require specific and significant adaptations of GMF runout models.

4.1 Identifying the location of forests with a protective function

The custom extension called “back-tracking” extension was implemented to adapt Flow-Py from a pure runout model to one that highlights the terrain associated with endangering infrastructures. Therefore, an additional input raster with the spatial distribution of infrastructure in the modeling domain is needed. Furthermore, minor changes must be implemented in the Flow-Py computer code, such that the path that the GMF traveled to reach the infrastructure is stored in computer memory.

4.2 Quantifying the protective effects of forests

To model the interactions between the movement of a GMF and forest, a custom “forest” extension was implemented in Flow-Py. The forest extension accounts for the increase in energy dissipation (or friction) in the parts of the GMF path that are located in forest by adjusting the runout angle to a steeper angle (Figure 2). The amount of change in the runout angle is dependent on the forested slope’s length, the forest’s structure and the kinetic energy height (relatable to velocity) of the GMF. That is, the forest extension adapts Flow-Py’s runout-angle stopping criterion to account for the forest’s ability to dissipate energy from GMFs.

To characterize the effect of different forest structures on the GMF movement, the forest extension uses the so-called forest structure index (FSI), which summarizes how developed a forest is with regards to its optimal protective effect and ranges between 0 (no protection) and 1 (optimal protection). The FSI needs to be provided as input for the Flow-Py simulation in form of a FSI forest raster. Because Flow-Py uses a statistical model, the parameterizations of forest-GMF interaction should ideally be based on observations quantifying the forest’s protective effect. In absence of data on how forests interact with different GMFs, the parametrization was developed in this first model development cycle based on a literature review and can be further refined with future observations.

We developed and applied Flow-Py in the Interreg Alpine Space project GreenRisk4ALPs [98] to model forest’s protective functions and effects on different GMFs. In the next subsections, the parameterizations established to model forest effects on rockfall, snow avalanches and shallow landslides using 10-m resolution raster (DEM and forest raster) will be presented. Since Flow-Py requires an already prepared starting zone raster as input, forest effects on avalanche and rockfall release were not explicitly considered. In terms of rockfall, we assumed that forest in the starting zone does not affect rockfall release. A method to quantify the relevance of the potential forest effect on the release of landslides has been developed as a pre-processing routine, which is an adaption to the input data showing the areas that are associated with potential landslide starting zones.

4.2.1 Modeling forest effects on rockfall propagation with Flow-Py

The recommended slope angle of ≥ 45° was used to identify rockfall release zones [108, 109]. The Flow-Py parameterization used to describe large single block rockfall can be found in Table 1.

Parameter	Value
runout angle	32°
divergence exponent	75
flux cutoff	0.03% of the initial flux (0.0003 m)
velocity limit	50 m s^−1 (maximum kinetic energy height ~ 130 m)

Table 1.

Flow-Py parameterization used to model large single block rockfall. The resulting behavior of the simulation is a single flow in steeper terrain and minimal spreading when the terrain flattens. The velocity limit keeps the speed of the rockfall in a range that has been observed in nature [32, 110].

The average energy loss along a rockfall trajectory can be modeled by friction parameters and is linear along a slope due to bouncing, rolling and sliding of falling rocks [111]. This linear relationship and normally distributed block stopping points along sections of the curve are displayed in Figure 3. The maximum increase to the runout angle due to forest is 13° (when FSI = 1 and kinetic energy height ~ 0), which leads to a maximum runout angle of 45° (mean runout angle of 32° + 13°). However, if the forest structure is not optimal the increase to runout angle is scaled to the FSI value (13° × FSI). That is, as the kinetic energy height increases, the less the protective effect of the forest. To incorporate this into the model the runout angle is reduced linearly by the kinetic energy height until the rock has a velocity of 30 m s⁻¹, which is well above 20 m s⁻¹ that have been found for forest to have an influence on rockfall propagation [112]. The optimal protective forest for rockfall is a mixed or broad-leaved forest (FSI = 1) with a high stem density while dense evergreen forest has a slightly reduced effect (FSI = 0.8). We assume that at a rockfall velocity ≥ 30 m s⁻¹ forest has no energy dissipating effect on a falling block (see Figure 3 for justification).

Figure 3.

Relationship between rockfall velocity (lower x-axis, red line) and the amount of forest-rock interaction, which is expressed in an increase to the runout angle (left y-axis). The upper x-axis (blue line) shows the relationship between kinetic energy height and the amount of forest-rock interaction. The green and orange bars are velocity ranges where forest has been found to affect rockfall propagation (green [112]) or rockfall velocities measured in forests (orange, *[39, 108, 113, 114, 115, 116]). The red dashed lines show measured rockfall velocities outside of forest [110, 117]. The blue dashed line shows that at 20 m kinetic energy height a forest can have a significant influence on the runout of a rockfall [112].

4.2.2 Modeling forest effects on snow avalanches with Flow-Py

Avalanche starting zones were delineated based on a slope inclination ranging between 28° and 55°, which is most often applied to identify slab avalanche starting zones [20]. The parameterization used to describe large dry-snow slab avalanches with Flow-Py is given in Table 2. We applied a mean runout angle of 25°, which was determined from measurements of 89 documented large avalanche events in Austria [19]. The maximum kinetic energy height limit imposed keeps the speed of the avalanche in a range that has been observed in nature [118], and which can be thought of as describing the turbulent friction of the avalanche debris.

Parameter	Value
runout angle	25°
divergence exponent	8
flux cutoff	0.03% of the initial flux (0.0003 m)
velocity limit	70 m s^−1 (maximum kinetic energy height ~ 250 m)

Table 2.

Flow-Py parameterization used to model large dry-snow slab avalanches. The resulting behavior of the simulation is a moderate spreading in steeper terrain and a more divergent flow in flat terrain.

Similar to rockfall, we could not find publications addressing the relationship between avalanche kinetic energy height (or velocity) and the increase to the runout angle (energy dissipation) in forested terrain and, therefore, assumed a linear relationship (see Figure 4 for justifications). The kinetic energy height is related to the square of the velocity which results in the relationship between the runout angle and velocity shown in red in Figure 4. The maximum increase to the runout angle in forest is 10° (when FSI = 1 and kinetic energy height ~ 0), which leads to a maximum runout angle of 35° (mean runout angle of 25° + 10°). However, if the forest structure is not optimal regarding avalanche protection the increase to the runout angle is scaled with the respective FSI value (10° × FSI). Furthermore, as the kinetic energy height increases the forest has less of a protective effect and, therefore, the runout angle is reduced linearly by the kinetic energy height until the avalanche has a velocity of 30 m s⁻¹ (kinetic energy height ~ 105 m). This threshold is reasonable when compared to the Swiss classification according to avalanche impact pressures and potential damages [121], i.e., at 30 m s⁻¹ an avalanche with a snow density of 200 kg m⁻³ would have an impact pressure of 170 kPa, which is a higher impact pressure than is needed to destroy large forest areas and uproot evergreen conifer trees. At a kinetic energy height of ~ 105 m (or velocity of 30 m s⁻¹) the forest is no longer capable of reducing the avalanches energy [51, 122].

Figure 4.

Relationship between avalanche velocity (lower x-axis, red line), and the amount of forest-avalanche interaction, which is expressed in an increase to the runout angle (left y-axis). The upper x-axis (blue line) shows the relationship between kinetic energy height and the amount of forest-avalanche interaction. The dashed red lines show values of documented avalanche velocities [118, 119, 120]. Colored points are calculated velocities (and kinetic energy heights) based on the Swiss classification according to impact pressures and related damages [121].

An evergreen conifer forest with a high stem density and dense canopy cover (FSI = 1) can significantly reduce runout lengths of small-to-medium avalanches [47], where broad-leaved forest (FSI = 0.8) has a reduced effect as well as a forest with lower stem densities and less dense canopy cover, which needs to be reflected in the choice of FSI-values. The applied method of adjusting the runout angle not only accounts for increasing the energy dissipation in avalanche transit zones but also somewhat for forest effects in starting zones since the maximum runout angle increase in well-developed evergreen conifer forest is set to 10°, which leads to an effective runout angle of 35°. That is, applying an increase to the runout angle dependent on the FSI adjusts the starting zones in forests since an avalanche might not be released on forested slopes up to 35°, which is a conservative estimate for regional modeling and based on the assumption that small forest openings might still be present in forests growing in steep avalanche terrain [46, 52].

4.2.3 Modeling the relevance of potential forest effects on shallow landslides with Flow-Py

There are major differences between the protective forest model used for shallow landslides, and the one used for rockfall and snow avalanches. This is due to the lack of data and process understanding on how forest interacts with landslides (i.e., open hillslope debris flows) in the transit zone at the resolution and/or detail required for regional modeling (see chapter [3] of this book). We therefore assumed no forest effect in the transit and runout zones of landslides, but forest in the starting zones reduced the likelihood (probability) of landslide release. Flow-Py was then used to highlight locations on a slope and in the valley bottom that are being protected by reducing the landslide’s release probability (the applied Flow-Py parameterization can be found in Table 3).

Parameter	Value
runout angle	22° [17]
divergence exponent	75 [18]
flux cutoff	0.03% of the initial flux (0.0003 m) [7]
maximum kinetic energy height	~ 12 m [7]

Table 3.

Flow-Py parameterization used to model large shallow-seated landslides. The resulting behavior is low spreading in steeper terrain and a more divergent flow when the terrain flattens.

Landslide starting zone modeling was performed by building upon the principle ‘the past is the key to the future’. A binary supervised classification algorithm was used to link past landslide locations and landslide-absence locations (topographically corrected random sample; > 50 m from past landslides) with a variety of topographical and thematic features (e.g., slope angle, topographic wetness index, relative topographic position, aspect, landform variables, and geology). The result is a raster that highlights potential landslide starting zones and breaks the areas down into likelihood categories based on how similar the areas are to past landslide release areas in the study region. This was done by a Generalized Additive Model that accounts for the previously observed non-linear relationships among landslide occurrence and the explanatory variables [81, 123]. The derived statistical model was applied to spatially predict the likelihood of class-membership of each raster cell that contains information on the environmental explanatory variables. This method of starting zone modeling is more sophisticated than that used for rockfall or snow avalanches; however, the applicability of this type of models is limited to the region of the input data and therefore is not a general model and not easily transferable to other locations.

The potential forest effect was considered at the level of the classified landslide release susceptibility (likelihood). Areas identified as favorable for landslide release (different classes) were grouped into forested terrain and non-forested terrain using a matrix-based approach. Areas susceptible to landsliding due to their topography and geology were considered more stable (without assigning a specific degree of stabilization) in case they were covered by forest compared to their non-forested counterparts. This builds upon the literature-supported assumption that stabilizing effects of trees usually outweigh destabilizing effects for shallow-seated landslide processes (e.g., [82, 96, 104, 124, 125]). Furthermore, we assumed that a forest located on a landslide-prone slope is more “relevant” compared to a forest on unsusceptible terrain (e.g., flat terrain or too steep terrain). For more details on how the forest effect on shallow landslide release was considered can be found in the GreenRisk4ALPs project report [126].

4.3 Post-processing of Flow-Py simulation results for risk-based decision support

To locate protective forests or quantify their protective effects by reducing impacts of gravitational natural hazards on people and assets several post-processing steps of the Flow-Py simulation results are needed.

4.3.1 Forests with an object protective function

The workflow to identify forests with an object protective function is shown in Figure 5 (see also subsection 4.1). The output of Flow-Py simulations with the back-tracking extension are spatially explicit subsets of the GMF starting, transit and runout zones that are associated with endangering infrastructure. A union between the back-tracking results and the spatial distribution of forested areas is then performed in a GIS resulting in a gridded dataset of forests with an object protective function, i.e., “Direct Object Protective Forest” for each hazard type. In chapter [1] of this book, we further explain these maps and discuss their use as decision support tool in risk-based protective forest and ecosystem-based risk management of natural hazards.

Figure 5.

Post-processing workflow to identify forests with an object protective function producing “Direct Object Protective Forest” maps.

4.3.2 Forests’ protective effects and regional impacts on gravitational natural hazards

The term protective effect describes the forest-GMF interaction, which we model with Flow-Py and its forest extension by increasing the runout angle in forested terrain. In contrast, we introduce the term forest impact(s) to describe how forest affects the spatial distribution or intensity of GMFs on a regional scale and their impacts on, e.g., infrastructure.

There are many types of forest impacts such as reductions in runout length, energy or mass of the GMF, or of its release probability, which can be quantified with different methods (see chapters [66, 127, 128] of this book). Within the GreenRisk4ALPs project, we developed the Impact Reduction Index (IRI), which shows the difference in average kinetic energy height between Flow-Py simulations with (employing the forest extension) and without forest effects. The post-processing workflow to calculate the IRI is outlined in Figure 6, and further examples and explanations of how the IRI can be applied for supporting decisions in protective forest and natural hazard risk management are provided in chapter [1] of this book.

Figure 6.

Post-processing workflow to calculate the Impact Reduction Index (IRI). Adapted from [126].

It is important to note that several starting zones may route flux through the same forest area and, since the forest effect depends on the kinetic energy height (GMF velocity), the same forest area can have a high protective effect associated with one of these GMF starting zones but provide almost no protection for the others. To reflect this in the IRI, we normalized the difference in average kinetic energy height between Flow-Py simulations with and without forest effects with the maximum value of kinetic energy height of “no-forest” Flow-Py simulation results at a given location.

5. Conclusions

Models and simulation tools for gravitational natural hazards combine different bits of process understanding in a way that is useful for scientists and practitioners. The predictive power of protective forest models results from combining the state-of-the-art process understanding of GMFs with the state-of-the-art process understanding of forest-GMF interactions in one modeling approach. Knowledge about GMF processes and their interactions with forest can be derived by physical or statistical means; however, it is often a combination of both that is included in one model (e.g. [129]). By combining these model types scientist can quantify how well the current process understanding describes reality and practitioners can use GMF and protective forest models and their results as decision support tools. Examples of the results of modeling forests’ protective functions and effects with Flow-Py that we described in this chapter, and which were obtained within the GreenRisk4ALPs project for rockfall, snow avalanches and landslides, as well as their application for decision support in natural hazard risk management are summarized in chapter [1] of this book. However, many of the presented parameterizations that were established in this first round of the model development cycle would benefit from further empirical parameter studies. The general lack of model parameterizations describing forest-GMF interactions is a major knowledge gap and should be addressed in future studies with the goal of obtaining a deeper process understanding that can be translated into enhanced parameterizations of, e.g., the four parameters that describe GMFs’ behavior as well as the forest structure index (FSI) applied in Flow-Py.

Acknowledgments

This work was conducted in the context of the GreenRisk4ALPs project (ASP635), which has been financed by Interreg Alpine Space, one of the 15 transnational cooperation programs covering the whole of the European Union (EU) in the framework of European Regional policy. We are grateful to the late Karl Kleemayr, who significantly contributed to the ideas behind the modeling of forest’s protective effects with Flow-Py as well as to the development of the Impact Reduction Index (IRI).

Conflict of interest

The authors declare no conflict of interest.