1. Introduction
Mountainous areas are prone to many mass movement processes, rockfalls being one of the most common [1]. Rockfalls are defined as the detachment of one or several fragments of rocks (blocks) from rock cliffs, followed by rapid down-slope movements by falling, bouncing, rolling, and sliding [2]. Rockfalls can be an important threat to infrastructure, human life, and their property, since they are rapid processes with long runouts [3], and their instantaneous occurrence makes their temporal prediction practically impossible [4]. Rockfall models, especially process-based ones, can be an efficient tool for predicting potential rockfall hazard areas, making it possible to identify rockfall source, propagation, and runout areas as well as model rock trajectories, the kinetic energy of rocks, rock rebound heights and propagation and reach-out probability [1,5,6,7]. By quantifying the potential rockfall hazard, simulation models can be used for planning different protection measures (e.g., technical measures, natural-based solutions) that can significantly reduce the potential risk of rockfall occurrence in high-threat areas [7,8,9,10].
Forests provide a natural solution for protection against rockfalls in alpine regions [11,12], since they can significantly reduce the intensity (kinetic energy) and propagation probability of falling rocks [13,14,15]. Although several rockfall modelling approaches have been proposed in the last two decades [16,17,18,19,20,21], only a few consider the protective effect of forest [19,21,22]. Models that do consider the protection effect of forest can additionally be used for mapping protection forest and quantifying its protection function against rockfalls [23]. In simulation models [19,21], the protective effect of forest is considered via the spatial distribution of the forest stand, stem diameter at breast height (DBH) distribution, and tree species. Simulation models enable rockfall hazard mapping with and without the effect of forest, and with different forest scenarios [23,24], and thus, they can be used to quantify the protective effect of forest and different forest stands.
Some rockfall models are intended for use on a local scale (individual rockfall events), while others are more appropriate for regional scale (multiple rockfall events on different slopes) [5,18,25,26]. The scale of modelling is an important factor when determining the purpose of the modelled results, as it will affect the choice of the grid cell size of the input data. The most common input data for rockfall simulation models is the digital terrain model (DTM), which carries information about the main morphological properties of the surface. The grid cell size of the DTM can have a significant impact on modelling both potential rockfall source and runout areas, primarily on the lateral dispersion of rockfall trajectories [7,27,28] and on the kinetic energy of rocks [29,30]. In order to achieve the most realistic results, the highest grid cell size is desired in rockfall modelling; however, since individual rockfall models are designed for use at a particular spatial level, a higher DTM grid cell size does not necessarily provide improved modelling results. Different studies have shown that the preferred DTM grid cell size lies between 2 and 10 m [31], since with DTMs that have a grid cell size larger than 10 m, the modelled rockfall trajectories will not be modelled correctly [32]. The impact of the grid cell size will be the most significant in areas with meso-scale roughness, while in the smooth areas, the difference is not as evident [30]. With larger DTM grid cell size, the delineation of rockfall source areas becomes rougher, overestimated, and enhanced, while the majority of the cliffs will still be detected [33]. Modelling with DTMs grid cell size smaller than 5 m can result in too much detail and artefacts in DTM [6,34,35], while with larger grid cell sizes, the artificial features such as roads and railway tracks will not be captured [29]. When considering the cost–benefit of rockfall modelling [36], it is also important to abide by the availability of DTM data and the computational time needed for modelling. The computational time will be longer with DTMs with a small grid cell size, with the difference being the largest when modelling rockfall propagation and runout areas [35]. Although the above-listed studies have investigated the influence of different grid cell sizes of DTM on rockfall modelling, none of them focused on analyzing the changing protection effect of forest.
Accordingly, this paper analyses the impact of grid cell size of DTM on the protective effect of forest in rockfall propagation and runout areas. The aim was to quantify the effect of forest on rockfall runout by considering several DTM grid cell sizes, and to observe how it influences modelling of the rockfall propagation and runout area. For modelling the rockfall propagation area, the Rockyfor3D model was used [21], and modelling was performed at different DTM grid cell sizes (1, 2, 5, 10 m). The model was applied with and without the effect of forest, while the forest stand structure remained the same. Based on this analysis, we investigated: (i) the effect of DTM grid cell size on the modelling of rockfall propagation and runout areas, (ii) the extent to which forest reduces rockfall propagation at different DTM grid cell sizes, and (iii) the influence of forest and DTM grid cell size on maximum kinetic energy and maximum passing height of rocks.
2. Materials and Methods
2.1. Study Area
Our study involved modelling at two rockfall sites in the Trenta valley, which are located in the Julian Alps in the northwestern part of Slovenia (Figure 1a,b). Geotectonically, the valley belongs to the upper part of the Southeastern Alps or thrust unit of the Julian Alps [37,38]. The thrust unit consists of sediment rocks from the Lower Triassic and Cretaceous with carbonate rocks from Upper Triassic being predominate. The dominant direction of faults is transverse dinar (NE-SW), while prominent faults are also present in the Dinaric direction (NW-SE) and are shifted along the transverse dinar faults. Less significant faults are in the direction of N-S and E-W. The prevailing bedrock in the valley, which is also present at the locations of both rockfalls, is the Dachstein limestone of the Norian–Rhaetian age. The limestone layers are 0.5 to 5 m thick. Typically, limestone is passing into dolomite in vertical and lateral direction [39]. Tectonically damaged rocks form very steep and even subvertical slopes that are covered by the Anemono trifoliate-Fagetum forest community [40], which predominantly has the protection function.
Rockfalls in the studied area mostly occur due to the freeze–thaw cycle occurring in fractured bedrock in the spring, and extreme rainfall especially in autumn [41,42,43]. As Slovenia lies on the seismically active area between the Alps and the Mediterranean Sea, where the Adriatic/Apulian Sub-plate hits the Eurasian tectonic plate, forming the Periadriatic Seam, rockfalls in this part of the Alps can also be earthquake-induced [44,45]. At the first rockfall site (RS1) (Figure 1c), the last large rockfall event occurred in April 2017, when approximately 29,000 m3 of coarse material was released. RS1 is located between 940 and 1160 m a.s.l. (Figure 1e), and its runout covers an area of approximately 19,300 m2. The source area is located at a vertical rocky cliff (fall height ≈150 m), and the majority of the material was deposited directly below the cliff, while individual fragmented blocks were deposited further down the slope. Previous rockfall activity can be distinguished by the color of the rocks (older rocks are darker), and by vegetation cover (older rocks are overgrown by low shrubs and dwarf pines). A similar situation is present at the second rockfall event (RS2) (Figure 1d), which lastly occurred in March 2020 when approximately 11,000 m3 material collapsed from a rock cliff (fall height ≈50 m) and impacted the area of 9500 m2. RS2 is located between 890 and 1010 m a.s.l. (Figure 1f). Previous activity is also noticeable by the older rockfall deposits and the rock-induced damages in the surrounding forest. Rockfalls at this site occur in minimum on average every 8.5 years [46]. Both rockfalls are located in the vicinity of housing and infrastructure; consequently, the forest on these slopes plays an important role in reducing the rockfall runout.
2.2. Rockfall Modelling
The Rockyfor3D rockfall model [21] was used for modelling rockfall propagation and runout area. Rockyfor3D is a probabilistic, process-based rockfall trajectory model of falling blocks in three dimensions [21] that can be used for regional, local, and slope-scale rockfall simulations. Rockfall trajectory is simulated as 3D vector data by calculating sequences of classical parabolic free falls through the air, rebounds on the slope surface, and also impacts against trees (optional). In the model, rolling is represented by a sequence of short-distance rebounds, while rock sliding is not modelled. The required input data include the topography and surface slope characteristics as well as a set of parameters that define the release conditions. The minimum input data required include the DTM, definition of the source area, rock density, rock size and shape, surface roughness, and soil type [21]. Rockyfor3D also enables simulation with forest, which can be done either (i) by providing a text file with locations of the trees, their DBH, and the percentage of coniferous trees, or (ii) by providing four raster maps containing the number of trees, mean DBH, standard deviation of DBH, and the percentage of coniferous trees. Namely, based on the DBH and the tree type, the model calculates the maximum amount of kinetic energy that could be absorbed and dissipated by a tree. The amount of maximum kinetic energy is determined as following:
(1)
where is the maximum amount of kinetic energy that can be dissipated by the tree (J), is the fracture of energy ratio of the tree type, and is the stem diameter at breast height. Since the broad-leaved trees are more resistant to breakage due to rockfall impact than coniferous trees, the model uses a higher for broad-leaved trees ( = 1.59) than for coniferous trees ( = 0.93).Additionally, the model enables simulation with rockfall nets as protection structures on a slope. The main outputs of the model are maximum kinetic energy (90% confidence interval of all maximum kinetic energy values), maximum bounce height, the number of block passes through each cell, rockfall propagation probability, the number of deposited blocks, maximum simulated velocity, maximum tree impact height, and the number of trees impacts per cell [21].
2.3. Used Datasets
2.3.1. Simulation Settings and Block Definition
Rockfall source areas were in the case of both rockfall areas determined by comparing the point clouds before (national LiDAR point cloud from year 2014 and survey with unmanned aerial vehicle after the event) and after the rockfall event. In the case of RS1, the photogrammetric point cloud was obtained from the UAV survey in July 2018, while in the case of RS2, the survey was done in March 2020. The area of source area at RS1 encompasses approximately 4000 m2, and at RS2, it encompasses approximately 330 m2. The rock density of both source areas was set to 2500 kg/m3 (limestone) [47], and based on the height difference between the source and runout areas, it was calculated using the LiDAR DTM1 data, and the initial fall height was set to 50 m in both cases. Block dimensions and block sizes are based on the field measurements of rocks at the maximum runout of the rockfall. At RS1, the block dimensions that were used for modelling were 1.4 m (d1), 0.9 m (d2), and 0.8 m (d3), and at RS2, they were 1.9 m (d1), 1.5 m (d2), and 1.1 m (d3). The prevailing block shape at both sites was set to rectangular. Variation in block volume was set to ± 0% in order to specifically focus on studying the effect of grid cell size and effect of forest on rockfall modelling, and neglecting variations in shape and size for a given block diameter. Throughout the paper, we use the term “block” for defining the simulation of single, individual pieces of falling rocks within the model, since the same term has been used by the author of the model [21]. This term does not necessarily reflect sedimentological nomenclature where rock blocks would, based on their dimensions, be treated i.e., in planetary and space research as large clasts, divided to boulders (<1 m) and megaclasts (blocks 1–10 m; megablocks 10–100 m; superblocks > 100 m) [48,49,50,51,52].
The input data that were of interest in this study were the DTM. In the Rockyfor3D guidelines, it is stated that both the spatial precision of the simulated maps and the accuracy of the simulated kinematics decrease with increasing cell size [21]. In order to test the preferred grid cell size of the data, we used the following set of DTM grid cell sizes: 1 × 1 m (DTM1), 2 × 2 m (DTM2), 5 × 5 m (DTM5), and 10 × 10 m (DTM10). Different grid cell sizes of DTM were created based on the LiDAR 1 × 1 m point cloud [53] obtained in 2014, using the binning interpolation type with average elevation values in ArcGIS Pro [54].
2.3.2. Surface Roughness Parameters and Soil Types
Surface roughness parameters (rg) represent rocks lying on the slope that form obstacles for the falling rocks [21]. The parameters define the surface roughness, which is expressed as the size of the material covering the slope’s surface in the downward direction of the slope. Rg70, rg20, and rg10 correspond to 70%, 20%, and 10% of the cases during a rebound on the slope, respectively, and they represent values from 0 to 100 meters (0 represents a smooth surface).
Determination of rg subareas was mapped directly on the field at each rockfall site (Figure 2). At RS1, we have determined three different rg surface slope type areas (Figure 2a): the first area represents the scree terrain, which is built from material from previous rockfall events; the second area is the area with dwarf pines and other coniferous tree species that are overgrowing rocky terrain, and the third area is the area of continuous forest area overgrowing steep terrain represented by surface rockiness, rock deposits, and laying logs. Soil type was at all areas set as type 4 (talus slope with Ø > ~10 cm, or compact soil with large rock fragments). At the RS2 (Figure 2b), four different rg subareas were determined: the first area represents scree material; the second area represents the terrain within the mixed forest that is characterized by surface rockiness, rockfall deposits, laying logs, and stumps; the third area is a meadow within the spruce forest and open-planed meadow; and the forth area is the river torrent. Correspondingly, two types of soil types were set; rg area 1, 2, and 4 were set as type 4 (talus slope with Ø > ≈10 cm, or compact soil with large rock fragments), while rg area 2 was set as type 1 (fine soil material with depth ≥ 100 cm). Based on the guidelines of the model, we defined the main rg values for each rg area at both rockfall sites (see Table 1); however, in order to calibrate the model, ranges of rg values were used in the calibration process for each rg area separately (see Section 2.4). The number of simulations was set to 1000, as recommended by the authors of Rockyfor3D.
2.3.3. Simulation with Forest
The modelling was done also with considering the scenario with the effect of forest. Simulation with forest was done by using tree locations, DBH, and percentage of coniferous trees. Tree locations and additional attributes were using the tool FINT [55]. The tool calculates the positions of dominant and co-dominant trees based on high-resolution digital surface models (DSM); namely, it evaluates for each grid cell in a normalized surface model (NSM), which has a value larger than the defined minimum tree height, whether it is a local maximum by initially using a 3 × 3 window. FINT assesses the dominance of the analyzed cell over its surrounding cells [55]. The recommended grid cell size are 0.5, 1, and 2 m, while with the increasing grid cell size, the spatial precision of the calculated tree positions decreases. The tool uses a standard function for calculating the DBH for each identified tree on the basis of its height:
(2)
where DBH (in cm) presents diameter at breast height and h the height of the tree (in m). However, the user can use a custom DBH function based on the relationship between the DBH, tree height, and tree altitude. The location of the trees in provided in a .txt format, and it contains X and Y coordinates for each tree along with the DBH.In the case of this study, we used DTM and a digital surface model (DSM) with a grid cell size of 1 m. The function for calculating DBH (in cm), based on the height, was adjusted based on the DBH and height measurements of trees present at both study sites. Fifty trees were measured at each site, 50% coniferous (Picea abies Karst. and Larix decidua Mill. at both sites, additionally Pinus sylvestris L. at the first rockfall site), and 50% deciduous tree species (Fagus sylvatica L. at both sites) of various DBH ranges. DBH functions were extracted separately for each rockfall site (R2 = 0.59 for RS1, R2 = 0.79 for RS2):
DBH (RS1) = 10.70 × e(0.0529 × h),(3)
DBH (RS2) = 7.87 × e(0.0529 × h).(4)
In addition to the location of the trees, the model also needs information about the ratio between broadleaved and coniferous trees; namely, when modelling with forest, one must provide a raster file where cell values represent the mean percentage of coniferous trees (%) within each raster cell. The mean percentage of coniferous trees for both locations was extracted from the stand map of the Slovenian Forest Service [56], and the spatial extents were additionally checked and corrected by the inspection on the field and by using orthophoto images. At both rockfall site areas, we have determined three areas with different percentage of conifers (Figure 2); at RS1, an area with 0% (no forest), an area with 30% conifers, and an area with 60% conifers, and at RS2, an area with 0% (no forest), an area with 40% conifers, and an area with 90% conifers.
2.4. Calibration of Surface Roughness Parameters
As the model is sensitive to the surface roughness parameters (rg70, rg20, and rg10) [21,57,58], the calibration of these parameters was done using a range of rg values within an individual rg area where rg70, rg20, and rg10 values were changing randomly between rg and rg subareas. For the calibration, we have used 30 different combinations of rg values at each subarea with a different surface roughness in each rockfall site, and the range of rg values was randomly chosen from the intervals that were set during the field inspection. For each subarea, we have selected a typical rg value for different slope surface types as provided by the authors of the model. Around that value, we have created an interval; in the case of RS1, ±0.4 to the initial rg value, and in the case of RS2, ±0.3. The range of rg values was larger in the case of RS1, since that rockfall impacts a larger area where subareas are more diverse and also larger when compared to RS2. Then, this interval (Table 1) represented the range of the rg values within which the random combinations of rg values for all subareas were created at the same time, and that the same pairs were later used in the model for simulation. Random values were created in R [59] using the function runif. The values within the same rg parameter (70/20/10) were unique and could not be used multiple times in the simulations. Calibration of the model was performed using DTM2, since the preferred grid cell size for the modelling lies between 2 × 2 m and 10 × 10 m [21,31], and by the experience of the author of the model, the 1 × 1 m grid cell size does not necessarily increase the quality of the results, and it increases the amount of data substantially [21]. In the calibration process, we also considered the effect of forest.
2.5. Evaluation of Calibration and Validation
Evaluation of the performance of each calibration run was done by using goodness-of-fit indices (GOF) as presented by Formetta et al. [60]. GOF indices are based on pixel-by-pixel comparison between the observed rockfall area map (OR) and predicted rockfall area map (PR). Comparison of these two maps results in binary maps with positive values corresponding to “actual rockfall area” and negative values corresponding to “not a rockfall area”. Correspondingly, four types of pixel outcomes are possible for each raster cell: (i) true positive (TP) is a pixel mapped as an “actual rockfall area” on both the OR and PR (correct prediction), (ii) true negative (TN) is a pixel mapped as “not a rockfall area” on both the OR and PR (correct detection of areas where rockfalls do not occur), (iii) false positive (FP) is a pixel that is actually “not a rockfall area” on the OR but is mapped as an “actual rockfall area” on the PR (false alarm), and (iv) false negative (FN) is a pixel that is an “actual rockfall area” on the OR but is mapped as “not a rockfall area” on the PR (missed alarm) (summarized based on Formetta et al. [60]). These indices are the basis of the concept of receiver operator characteristics (ROC) [61] that is used for assessing the model performance using the relation between benefits (TP) and costs (FP). Formetta et al. [60] incorporated eight GOF indices in the ROC system for quantification of model performance; however, four indices have been shown to be the most suitable for the evaluation of calibration runs and were used in this study (Table 2): the success index (SI), distance to perfect classification (D2PC), the average index (AI), and true skill statistics (TSS). More comprehensive and detailed descriptions of the indices are available in Formetta et al. [60].
The calibration run that achieved the most optimal value with those indices was selected as the most successful, and the rg values of that calibration run were used for modelling rockfall runout area at other DTMs considering two modelling scenarios: with and without forest. The actual rockfall area was mapped on the field by using a mobile application for collecting past rockfall deposits [62]. For the modelled rockfall area, the output propagation probability was used both in the calibration and validation step, since this raster layer represents the most realistic spatial distribution of the current rockfall event and can be used for calculating spatial occurrence probability, which is used in rockfall hazard analyses [21]. The same GOF indices as in calibration were also used for validating the performance of the modelling at DTMs with different grid cell sizes and for evaluating the effect of forest, following the same approach as for calibration.
2.6. Comparison of Additional Output Data
In order to support the validation procedure, especially focusing on the effect of forest on modelling rockfall propagation and runout zone, additional outputs of the simulation model were analyzed and argued in the results, namely: the surface area of the modelled area with and without forest, and achieved maximum kinetic energies (E_95CI—the 95% confidence interval of all maximum kinetic energy values in kJ recorded in each cell) and maximum passing heights (Ph_95CI—the 95% confidence interval of all maximum passing height values in m, measured in normal direction to the slope surface). The E_95CI and Ph_95CI are in the model calculated for the entire modelled runout area, and they can be used for dimensioning rockfall protective measures. With that reason, the results of these two outputs are shown only for the selected parts of the slope. We have placed fictive fences (screens) on the rockfall slope in order to assess the effect of forest and DTM grid cell size on the modelling. On each rockfall location, we have put one fence that crosses from the lateral part of the rockfall runout zone toward the central part (Figure 3). The length of the fence on the RS1 was 120 m, and it crosses the area of all real rockfall outline, and the fence at RS2 was 60 m long and it crosses only until the middle part of the rockfall outline. With these positions of the fences, we wanted to capture the momentum of the rockfall runout zone crossing between the forested and not forested area.
2.7. Model Sensitivity to Initial Fall Height and Variation of the Rock Volume
The model has been tested for its sensitivity on two input parameters: initial fall height and variation of the rock volume. The preference of the model in both case studies was to use 50 m as the initial fall height, since that is the real measured height difference within the rockfall source area at both sites, and a variation of volume of 0%, since we wanted to focus on the influence of the forest on the modelling. Nevertheless, beforehand, we ran the model with three initial fall heights (10, 20, and 50 m), and five variations of volume (0%, 5%, 10%, 20%, and 50%). The initial fall heights and variation of the rock volume values that were used in the sensitivity assessment are the values that one can choose from the model default values while setting up the input parameters. The evaluation of the model performances was done in a similar manner to that in the calibration and validation step.
3. Results
3.1. Model Calibration for DTM2
The results of calibration evaluated with GOF indices (AI, SI, D2PC, and TSS) are presented in Figure 4. The performance of the models was evaluated with each GOF indices separately, and in the end, they were combined into a common rank from best to worst, considering the model success of all four indices. Since RS1 had three subareas, and RS2 has four subareas, we combined them into one averaged rg70, rg20, and rg10 only for the visualization purpose on Figure 4. Averaged rg values were calculated by considering the rg70/20/10 value of each rg area and the area of which they actually covered in the modelling process.
The most successful model performance was in the case of RS1 achieved by calibration run number 17, and in the case of RS2, it was achieved by calibration run number 4. The real values of rg70, rg20, and rg10 for both RS are available in Table 1. At RS1, calibration run 17 was the most successful with three indices (not with D2PC), while at RS2, calibration run 4 was the most successful with all indices (Table 3).
Observing the three most successful model performances at each RS, it is possible to recognize that AI and SI achieve higher values that D2PC and TSS. The difference between the rockfalls could partially be due to the different combinations of rg values in subareas that could have been more or less suitable to begin with.
3.2. Validation of the Modelled Rockfall Runout Zones with and without Forest
The results of validation of models with different DTM grid cell size and scenarios (with/without the forest), evaluated with GOF indices, are shown in Table 4. The best model performance is in all cases achieved by the forest scenario. DTM1 achieves the best model performance with forest with both rockfalls, which is followed closely by DTM2. At RS1, DTM10 gets better results when forest is not considered in modelling. It is even, according to the GOF indices, performing the best of all models, which can be explained by the overestimation of the model due to the size of the DTM grid cell. At RS2, the differences between the effect of forest and without the effect of forest are small between all grid cell sizes.
The rockfall propagation and runout area decreases with smaller DTM grid cell size in both modelling scenarios at RS1, while it increases at RS2 (Figure 5). All models overestimate the actual rockfall extent at all grid cell sizes but at different locations. The largest overestimation of the propagation area at RS1 is at DTM1 and DTM2 (both scenarios) in the northern part of the rockfall propagation area, while in the southern part, it exhibits the greatest match between the modelled results and actual rockfall outline. On the other hand, overestimation of the model is not present as much at DTM5 and DTM10 (as in DTM1 and DTM2); however, both DTMs underestimate the maximum runout length. The underestimation of the propagation area is the lowest in the southwestern part with all DTMs except with DTM1, where it almost achieves a perfect fit in the forest modelling scenario. There is a difference in the modelled runout area between the northern part and southwestern part of the propagation area, which can be explained by the fact that in the northern part of the runout area, the surface slope is steeper, providing more kinetic energy to the rocks that can move further and cannot be stopped quicker due to the effect of forest. In the southwestern part, the runout distance is larger and changes in slope values are smaller; therefore, the difference in the DTM models’ results will be smaller. At DTMs with a smaller grid cell size, the changes in the terrain are captured in larger detail; therefore, the overestimation will be larger. Comparing the modelling results of RS1 to those of RS2, the differences between the shape of the modelled propagation area DTMs are not as evident, since both lateral sides of the propagation area have very similar surface slope and terrain features.
The forest has the largest impact on rockfall propagation area at RS1, namely at DTM2—reducing it by 13%, followed by DTM1 (12%) and DTM5 (4%) (Figure 5). In the case of DTM10, the rockfall propagation area did not change between the modelling scenarios. Differences between the forest and no forest scenario at RS2 are almost negligible; forest reduced the rockfall propagation area by 4% at DTM2, 3% at DTM1, and 2% at DTM5, while at DTM10, the propagation area does not change. When observing the forest and no forest modelling results at RS1, it can also be stated for RS1 that at DTM1 and DTM2, the shape of the forest scenario changes according to the no-forest scenario, and the change is only observed for the maximum runout distances. On the other hand, at DTM5 and DTM10, the shape of the rockfall runout area is modelled correctly in both modelling scenarios; however, there are no major differences between the forest and no-forest scenarios. In general, all DTMs at RS1 provide comparable results with respect to the shape of the rockfall propagation.
3.3. Comparison of the Maximum Kinetic Energy and Maximum Passing Heights between DTMs and Modelling Scenarios
Two fictive fences placed perpendicular on the slope direction were used to calculate the maximum kinetic energy (E_95CI) and maximum passing height (Ph_95CI). The results are shown in the form of a distance profile along both of the fences (Figure 6). Within the profile, it is possible to distinguish between two momentums: when the rocks in the model are in interaction with the forested area (trees) and when not. At both rockfall sites, the model is in interaction with forest within the first 40 m (RS1) and first 30 m (RS2), as marked on the distance profile.
Observing E_95CI results at both sites of the profile, the kinetic energy values experience high variations along the distance profile. The maximum kinetic energy in this part has smaller variations compared to the part of the profile that is not in interaction with the forest, and larger jumps in kinetic energy can be recognized at RS1 when observing the forested part. Focusing only on the forested part in RS1, the deviations between scenarios with/without forest can be recognized with the higher peaks being scenarios with forest. The DTM1 scenario with forest achieves the highest values of E_95CI compared to other DTMs (903 kJ forest, 650 kJ no forest scenario). The maximum kinetic energies are the lowest at DTM10. Differences between scenarios with and without forest are smaller at RS2 compared to RS1. Absolute maximum values are achieved by DTM2 scenario with forest (1180 kJ).
In the part of the slope with no forest, maximum kinetic energies are higher. At RS1, differences between scenarios with and without forest are small (in average 1.6%), except for DTM10 (26% larger values in case of forest scenario). DTM1 and DTM2 achieve the greatest values with scenarios when forest is included into the model. DTM10 stands out from others, as the maximum value of maximum kinetic energy is much lower along all the distance profile from other DTMs. More diverse curve lines are present at RS2. There are two that stand out—DTM5 and DTM10. DTM5 has completely different curves with scenarios with and without forest. While with scenarios with forest it gets the absolute highest maximum kinetic energy values, it gets the lowest in scenarios without forest in the beginning of the interval and the highest at the end, when the forest scenario part gets very low again. DTM10 experiences very drastic changes in values; scenarios with and without forest are completely different. Differences scenarios with and without forest are larger than at RS1 (for 114 kJ higher in average). The lowest are achieved by DTM1 (55 kJ).
The maximum passing heights at RS1 experience the largest deviations in the part where the runout is in the contact with forest. There is a high exchange between lower and higher passing heights (between 2 and 4 m) at DTM1, with maximum passing heights in general being smaller in scenario without forest (with some exceptions). With other DTMs, the passing heights are not fluctuating as much—DTM2 has small fluctuations in height (<0.5 m) in both modelling scenarios, while at DTM5 and DTM10, the maximum passing height at this distance from the source area does not change anymore, and is fixated at 1.7 m through the whole distance profile. When passing to the area mostly outside the forest, there are still deviations in the maximum passing height at DTM1; however, they are much smaller compared to the forested part of the slope. In addition, the maximum values are lower (up to 3 m), and in some parts, the values of the scenario without forest surpass the values of the forest scenario. From the forested part of the slope, the trend also continues for DTM2. At RS2 at this distance from the source area, there are practically no differences in maximum passing height between scenarios with and without the forest, and also with all DTMs. It is fixated between 2.4 and 2.6 m. Some peaks are at DTM1 at both scenarios where the height changes for less than 1 m, except in the middle of the slope where it drastically jumps.
3.4. Results of Model Sensitivity Assessment
Regarding the input parameters initial fall height and variation of volume variability, we have performed a sensitivity analysis (Table 5). We have concluded that the initial fall height has proven to have a more important effect on modelling than variation of volume variability; in the case of this study, the chosen initial fall height of 50 m (for both rockfall sites) has proven to be the most successful, since the smaller the fall height, the less successful the performance of the model. In the case of variation of volume variability, the following parameter does not significantly improve or reduce the performance of the model. Slightly more successful models were achieved when using a variation of volume of 50%, but the maximum difference in the modelled area of true positives (actual rockfall outline) was only 0.5% (forest scenario RS1)/0.22% (no forest RS1) and 0.9% (forest scenario RS2)/1.4% (no forest RS2). TSS and D2PC show the most sensitive behavior to the changes to input parameters. In the sensitivity assessment, we only used one DTM grid cell size, since we predicted the similar behavior with others, as well as we were only changing one parameter at the time, and it does not consider simultaneous variations between the parameters, which could result in different modelling success.
4. Discussion
In this study, we investigated the influence of DTM grid cell size and the protection role of forest on modelling rockfall propagation and runout area. For the calibration of the rg indices (at DTM2 forest scenario), we have used four indices (AI, SI, D2PC, TSS), even though the authors suggest that using one is enough for the parameter estimation. By calculating all indices, we conclude that there are differences between two groups of indices, namely AI and SI, and D2PC and TSS. In general, the first two provide the results that indicate more successful models achieving higher values, while the other two achieve lower values in general. The difference is most likely due to the formulation of AI and SI that give more weight to TN values compared to D2PC and TSS [60]. In the case of our study, TN presented a more “stable” part of the equation by not changing as much between the model calibration runs as the FP and FN values that would have a more important weight in the calculation of D2PC and TSS. In the end, the combination of all four by using ranks was the option that was the most secure to ensure that all aspects of correctly and wrongly modelled areas were considered in the calibration as equally as possible.
In order to capture the surface roughness in detail, the surface slope at both rockfall sites was divided into multiple subareas that were then given different rg values according to the mean obstacle height and percentage of coverage. Differences between values of rg70, rg20, and rg10 between subareas are small, all falling in the size of the surface roughness classes between 10 cm and 1 m. At both rockfall sites, the largest rg values were in subarea 2, and this area was the most crucial on reducing the maximum runout zone. Area 2 is located between the forested and non-forested area, where the majority of larger rocks from previous events were deposited due to the obstacles such as trees, stumps, laying logs, concavities, etc. The combination of rg values calibrated for area 3 at RS1 and especially at area 4 at RS2 is mostly a consequence of a good performance of combination of values at subarea 1 and 2, since it is located out of the reach of the majority of rocks in the model.
There could be some improvements in the calibration process. The calibration was done without the inclusion of soil types—based on the average rg values, soil type areas could also have been delineated additionally and calibrated together with the rg values [63], leading to the improved final result. When comparing the results of DTMs, we can observe that the maximum runout is underestimated at DTMs with a larger grid cell size. The same calibrated rg values were used at all DTMs; however, this indicates that perhaps the separate calibration of these parameters could be done for different grid cell sizes in order to achieve more realistic results. For example, for RS1—DTM5 and DTM10—it would have been better to use lower rg values because due to their cell sizes, the rocks in the model stop too soon, since the surface is “more rugged” according to the combination of DTM grid size and rg values. By decreasing the rg values, the surface would be smoother, allowing longer travel distances and reducing the underestimation of the model.
In general, all DTMs predict the shape of the rockfall area accurately; at RS1, DTM1 and DTM2 are overestimating the lateral extent, and DTM5 and DTM10 are underestimating the maximum runout, while at RS2, all models overestimate in the lateral direction, DTM5 and DTM10 additionally at maximum runout. The shape of the area is more complex and detailed at DTM1 and DTM2 as a result of greater interaction of the rocks and surface slope [29] in the model in comparison to DTM5 and DTM10, where it is possible to observe the impact of grid shape of raster on the actual rockfall runout shape [32,35]. All models overestimate the extent at the lateral side of the rockfall outline, but it is more expressed on the smaller grid cell sizes, which can be explained by smaller roughness at larger DTMs and lack of topographic channeling [64].
The effect of forest is the most evident at RS1, while it did not have an important role at the RS2. At RS1, forest reduces the maximum runout area at DTM1 and DTM2, and it also changes the shape of the modelled propagation area, while at the other two, it can practically be neglected. Even though we found field evidence that forest had a crucial role on stopping rocks at RS2, it did not have that effect in the simulation. The following is mostly due to the distance of the forest from the source area and the initial fall height with which the model gets high kinetic energy. The rockfall site is surrounded by forest, and in the areas where it should have reduced the maximum runout area or lateral extensions, the distance from the source is less than 30 m. When the rocks hit the area of forest within the model (Figure 6b), the maximum kinetic energy (between 1000–2000 kJ) is too high and the trees would not be able to dissipate that much energy [8,22]. At RS1, the distances to the forest are larger (>100 m), and the kinetic energy is lower (between 500 and 800 kJ) (Figure 6a) when rocks hit the forested area; therefore, forest is able to reduce their impact to some extent, as it can be observed in Figure 5.
Changes in the modelled rockfall runout area will not be as evident between different DTM grid cell sizes when the surface slope is more smooth and there are no major terrain features that stand out [30,35]. The same trend can be observed when modelling with or without forest. At the scale of one rockfall area only, DTM1 and DTM2 detect different movement of rocks within the model due to forest effect. The model used in this paper is more intended for use on the local scale, meaning that changes in the DTM grid cell size will be expressed quicker than in the case of the models that are intended for the use on the regional level. This is mainly due to the different inclusion of the terrain features into the model, where in the case of regional models, the topography is accounted for in a more general way (e.g., energy line angle principle) with all DTM grid cell sizes than by such complex 3D models such as Rockyfor3D.
Forest not only reduces the rockfall propagation and runout area at all DTM grid cell sizes but also reduces the maximum kinetic energy [30,64,65,66,67]. When forest is not included into model, maximum kinetic energies changes along the profile are quicker at DTMs with higher grid cell size (DTM5, DTM10), while with lower DTMs, the changes are more gradual. With larger grid cell sizes (DTM5 and DTM10), the slope geometry will be smoother and simplified, resulting in larger deviations in both kinetic energies, since there are less micro-topographic features that could absorb kinetic energy during the movement [29,68]. The observed differences between DTMs in the maximum kinetic energy at the area outside the forest at RS1 were smaller, indicating that in smoother areas, the influence of DTM grid cell size on the results is smaller than on more rugged surfaces [30]. If forest is included into the model, kinetic energies experience less changes in the maximum values, since the obstacles (trees) are the same for all models. Outside the forest, the situations change rapidly—kinetic energy is mostly still similar between the scenarios at DTMs with a larger grid cell size, but it changes with larger cell sizes where micro-topography is neglected (e.g., absence of channels and concavities) and the cell size can be averaged [68], leading to extreme jumps in kinetic energies as observed in Figure 6.
The maximum passing height was higher in the case of forest modelling scenarios (DTM1). This indicates that the rocks rebound at the impact with trees, and combined with well-expressed micro-topographical features, that means that there will be larger changes in the bounce height and changing the trajectories of the rocks. Details of terrain features can even be too detailed at DTM1, and along with the errors and artefacts that might be present in a DTM, this could lead to a possible overestimation of a passing height [29,34]. The changes in passing height at DTM1 can be very extreme, but once moving to smaller grid cell size (e.g., DTM2), the changes in passing height reduce drastically and do not differentiate as much between the forest and no-forest scenarios. Transitioning from the movement (high passing height to low) to the steady mode is achieved much faster in the areas with lower mean slope values and smoother slopes, which can in the case of this study be seen at DTMs with larger grid cell sizes (DTM5, DTM10), especially at RS2.
5. Conclusions
As already confirmed in the literature [24,28,69] and by the results of this study, forest should be considered as part of rockfall hazard assessment. However, forest is not always as efficient at reducing rockfall propagation area, maximum kinetic energies, or maximum passing heights, especially if the forest is too close to the rockfall source area or the rockfall event is too large to be stopped by forest. A combination of rockfall simulation tools that consider the forest protection effect can enable us to identify areas of forest that are insufficient for stopping rocks and where we consequently need additional protection measures [24,28,67,70]. These data can help with planning the location and capacity of the required protection measures (e.g., the height of the nets and their capacity) and the associated costs [24,71].
Nevertheless, the decision on the protective effect of forest and additional protection measures will be strongly correlated with the DTM grid cell size. Even though the grid cell size 1 × 1 m is not recommended by the model regarding the final results compared to larger grid cell sizes, based on our results, we can conclude that 1 × 1 m and 2 × 2 m were the most suitable grid cell sizes. They most accurately portrayed the protection effect of forest, while at other DTM grid cell sizes, the effect of forest in reducing rockfall propagation and runout area is practically neglected. If considering the computational time needed for modelling, 2 × 2 m is more efficient. Grid cell sizes of 5 × 5 m and 10 × 10 m would be too coarse in order to be used in planning of protection measures on a local scale. In the results, we can also observe that with those two DTMs, the maximum kinetic energy and bounce heights differ greatly from 1 × 1 m and 2 × 2 m grid cell sizes, and they express major differences along the rockfall slope. It is also important to state that a smaller DTM grid cell size might be very sensitive to small changes in elevation, which can be observed in the rockfall simulation model results as high fluctuation/variability of e.g., maximum kinetic energy or bounce heights in a small area. Yet, it is still more preferable to use 1 × 1 m or 2 × 2 m for the calculation of maximum kinetic energies and passing height than the other two, especially when planning protection measures, since the other two grid cell sizes provide reliable simulations on those data. Thus, the decision on DTM grid cell size must be consistent with the main goal of the final modelling results.
Author Contributions
Conceptualization, B.Ž., M.M. and M.K.; methodology, B.Ž., M.M. and M.K.; formal analysis, B.Ž.; data curation, B.Ž.; writing—original draft preparation, B.Ž.; writing—review and editing, M.M. and M.K.; visualization, B.Ž.; supervision, M.M. and M.K.; project administration, M.K.; funding acquisition, M.K. All authors have read and agreed to the published version of the manuscript.
Funding
This research is a result of the work on the Interreg Alpine Space project GreenRisk4ALPs (ASP 635) and was also supported by the Pahernik foundation. The second author acknowledges the financial support from the Slovenian Research Agency (research core funding No. P2-0180).
Institutional Review Board Statement
Not applicable.
Informed Consent Statement
Not applicable.
Data Availability Statement
The data presented in this study are available on request from corresponding author.
Conflicts of Interest
The authors declare no conflict of interest.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figures and Tables
Enlarge this image.Figure 1. (a) The study sites are located in the northwestern part of Slovenia in the Julian Alps. (b) Both rockfalls are a part of the Trenta valley. (c) Orthophoto image with marked source, and propagation and runout area of rockfall site 1. (d) Orthophoto image with marked source, and propagation and runout area of rockfall site 2. (e) Transverse profile of rockfall site 1. (f) Transverse profile of rockfall site 2.
Enlarge this image.Figure 2. Delineation into roughness parameters (rg) surface slope type areas into rg subareas at (a) rockfall site 1 (RS1) where there are three rg subareas, and (b) rockfall site 2 (RS2) where there are four types of rg subareas. Percentage (%) value on individual polygons indicates the percentage of conifers is that rg area. If one area has different percentage of conifers within, the polygon is additionally divided with the dashed line.
Enlarge this image.Figure 3. The position of the fictive fences on the slope for purpose of calculating maximum kinetic energy (E_95CI) and maximum passing height (Ph_95CI) for (a) RS1 and (b) RS2. The dashed lines indicate the position of the fictive fence according to the rockfall source area.
Enlarge this image.Figure 4. Models performance results in 30 calibration runs for (a) rockfall site 1 (RS1), and (b) rockfall site 2 (RS2) according to the rg70, rg20, and rg10. The values of rg70, rg20, and rg10 are for the demonstration purposes adjusted according to the values of rg within the individual subareas altogether, meaning that the values on the x, y, and z axis are actually showing average rg values that were calculated by considering an average rg value of each calibration run considering the subarea size. The real rg values were in the range shown in Table 1. Goodness-of-fit (GOF) indices that were used for the evaluation of calibration runs were average index (AI), success index (SI), distance to the perfect classification (D2PC), and true skill statistics (TSS), and on this plot, they are merged into a joint result, and the size of the sphere is reflecting the rank performance of each calibration run. The larger the sphere on the graph, the lower the rank, indicating a more successful model performance (rank 4 is the most successful model, rank 120 is the least successful). In the case of RS1, calibration run number 17 has the lowest rank of 8, and in the case of RS2, calibration run number 4 has the lowest rank of 4.
Enlarge this image.Figure 5. Modelling results for digital terrain model (DTM) grid cell sizes (1, 2, 5, 10 m) considering two modelling scenarios: with and without the effect of forest.
Enlarge this image.Figure 6. Calculated E_95CI and Ph_95CI values for (a) RS1 and (b) RS2 for the locations of two fictive fences located perpendicular to the surface slope direction. In the case of RS1, the location of the fence is 140 m from the initial point in the rockfall source area, and in the case of RS2, it is 90 m. The direction of the distance along the fence is marked in Figure 3.
The left side of the table is showing the subareas with different slope surface roughness with the interval values that were used for the creation of a random set of rg coefficients (30 combinations) for each subarea and parameter at both rockfall sites. The right-hand side of the table is showing the combinations of rg coefficients at each subarea that provide the most successful model performance in the calibration.
| ROCKFALL SITE 1 | ||||||
|---|---|---|---|---|---|---|
| Intervals Used for Creation of 30 Random Calibration Runs | Optimal rg Values Obtained in Calibration Run 17 | |||||
| sub area | rg70 | rg20 | rg10 | rg70 | rg20 | rg10 |
| area 1 | 0.31 |
0.6 |
0.7 |
0.14 | 0.36 | 0.65 |
| area 2 | 0.31 |
0.31 |
0.4 |
0.44 | 0.37 | 0.54 |
| area 3 | 0.31 |
0.35 |
0.4 |
0.13 | 0.44 | 0.53 |
| ROCKFALL SITE 2 | ||||||
|
Intervals Used for Creation of 30 Random
|
Optimal rg Values Obtained in Calibration Run 4 | |||||
| sub area | rg70 | rg20 | rg10 | rg70 | rg20 | rg10 |
| area 1 | 0.21 |
0.4 |
0.5 |
0.26 | 0.21 | 0.31 |
| area 2 | 0.35 |
0.4 |
0.45 |
0.5 | 0.53 | 0.26 |
| area 3 | 0 | 0 | 0 | 0 | 0 | 0 |
| area 4 | 0.3 |
0.3 |
0.35 |
0.17 | 0.32 | 0.53 |
Indices of goodness-of-fit (GOF) for comparison between actual rockfall area and predicted rockfall area, after Formetta et al. [60].
| Name | Definition | Range | Optimal Value |
|---|---|---|---|
| Success index (SI) | SI = ½ × ((TP/(TP + FN)) + (TN/(FP + TN)) | [0, 1] | 1.0 |
| Distance to the perfect classification |
D2PC = √((1 − TPR)2 + (FPR)2) | [0, 1] | 0.0 |
| Average index (AI) | AI = 1/4 × ((TP/(TP + FN)) + (TP/(TP + FP)) + (TN/(FP + TN)) + (TN/(FN + TN))) | [0, 1] | 1.0 |
| True skill statistics (TSS) | TSS = ((TP × TN) − (FP × FN))/((TP + FN) × (FP + TN)) | [−1, 1] | 1.0 |
The results of GOF indices for the three most successful model runs in the calibration step.
| ROCKFALL SITE 1 | ROCKFALL SITE 2 | ||||||||
|---|---|---|---|---|---|---|---|---|---|
| Calibration Run | AI | SI | D2PC | TSS | Calibration Run | AI | SI | D2PC | TSS |
| 17 | 0.855 | 0.837 | 0.285 | 0.673 | 4 | 0.892 | 0.892 | 0.153 | 0.785 |
| 27 | 0.820 | 0.824 | 0.250 | 0.647 | 25 | 0.886 | 0.886 | 0.162 | 0.772 |
| 23 | 0.820 | 0.822 | 0.253 | 0.644 | 3 | 0.885 | 0.886 | 0.167 | 0.771 |
Validation results for modeling scenarios considering different DTM grid cell sizes and two modeling scenarios: with and without the effect of forest.
| ROCKFALL SITE 1 | ||||||||
|---|---|---|---|---|---|---|---|---|
| scenario without forest | scenario with forest | |||||||
|
DTM grid
|
AI | SI | D2PC | TSS | AI | SI | D2PC | TSS |
| DTM1 | 0.927 | 0.957 | 0.076 | 0.914 | 0.893 | 0.934 | 0.126 | 0.868 |
| DTM2 | 0.927 | 0.947 | 0.076 | 0.894 | 0.897 | 0.931 | 0.110 | 0.862 |
| DTM5 | 0.927 | 0.930 | 0.110 | 0.859 | 0.922 | 0.931 | 0.103 | 0.862 |
| DTM10 | 0.922 | 0.932 | 0.100 | 0.864 | 0.928 | 0.939 | 0.090 | 0.877 |
| ROCKFALL SITE 2 | ||||||||
| scenario without forest | scenario with forest | |||||||
|
DTM grid
|
AI | SI | D2PC | TSS | AI | SI | D2PC | TSS |
| DTM1 | 0.920 | 0.940 | 0.089 | 0.879 | 0.918 | 0.940 | 0.092 | 0.880 |
| DTM2 | 0.924 | 0.929 | 0.106 | 0.857 | 0.914 | 0.924 | 0.109 | 0.848 |
| DTM5 | 0.908 | 0.927 | 0.104 | 0.853 | 0.903 | 0.923 | 0.110 | 0.846 |
| DTM10 | 0.908 | 0.933 | 0.107 | 0.867 | 0.902 | 0.926 | 0.114 | 0.853 |
The sensitivity analysis of initial fall height and variation of rock volume using both forest and no-forest scenarios.
| ROCKFALL SITE 1 | SCENARIO WITHOUT FOREST | SCENARIO WITH FOREST | |||||||
|---|---|---|---|---|---|---|---|---|---|
| AI | SI | D2PC | TSS | AI | SI | D2PC | TSS | ||
| initial fall height (m) | 50 | 0.925 | 0.982 | 0.032 | 0.963 | 0.945 | 0.976 | 0.046 | 0.951 |
| 20 | 0.936 | 0.925 | 0.149 | 0.849 | 0.937 | 0.914 | 0.172 | 0.828 | |
| 10 | 0.928 | 0.889 | 0.221 | 0.778 | 0.924 | 0.878 | 0.244 | 0.755 | |
| variation of rock volume (%) | 0 | 0.925 | 0.982 | 0.032 | 0.963 | 0.945 | 0.976 | 0.046 | 0.951 |
| 5 | 0.925 | 0.982 | 0.032 | 0.964 | 0.945 | 0.977 | 0.044 | 0.953 | |
| 10 | 0.925 | 0.982 | 0.031 | 0.964 | 0.946 | 0.978 | 0.041 | 0.956 | |
| 20 | 0.924 | 0.982 | 0.032 | 0.964 | 0.944 | 0.977 | 0.043 | 0.954 | |
| 50 | 0.924 | 0.983 | 0.030 | 0.966 | 0.945 | 0.979 | 0.039 | 0.959 | |
| ROCKFALL SITE 2 | SCENARIO WITHOUT FOREST | SCENARIO WITH FOREST | |||||||
| AI | SI | D2PC | TSS | AI | SI | D2PC | TSS | ||
| initial fall height (m) | 50 | 0.936 | 0.953 | 0.088 | 0.906 | 0.943 | 0.951 | 0.093 | 0.902 |
| 20 | 0.921 | 0.877 | 0.245 | 0.753 | 0.917 | 0.872 | 0.254 | 0.744 | |
| 10 | 0.909 | 0.846 | 0.308 | 0.691 | 0.905 | 0.838 | 0.323 | 0.676 | |
| variation of rock volume (%) | 0 | 0.936 | 0.953 | 0.088 | 0.906 | 0.943 | 0.951 | 0.093 | 0.902 |
| 5 | 0.940 | 0.955 | 0.086 | 0.909 | 0.939 | 0.949 | 0.098 | 0.897 | |
| 10 | 0.937 | 0.956 | 0.083 | 0.912 | 0.940 | 0.948 | 0.099 | 0.896 | |
| 20 | 0.942 | 0.958 | 0.079 | 0.916 | 0.936 | 0.942 | 0.112 | 0.884 | |
| 50 | 0.941 | 0.961 | 0.073 | 0.922 | 0.941 | 0.951 | 0.093 | 0.903 | |
© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Details
; Mikoš, Matjaž 2
; Kobal, Milan 1
1 Department of Forestry and Forest Renewable Resources, Biotechnical Faculty, University of Ljubljana, Večna pot 83, 1000 Ljubljana, Slovenia;
2 Department of Environmental Civil Engineering, Faculty of Civil and Geodetic Engineering, University of Ljubljana, Jamova cesta 2, 1000 Ljubljana, Slovenia;