Delimiting rockfall runout zones using reach probability values simulated with a Monte-Carlo based 3D trajectory model

At present, a quantitative basis for delimiting realistic rockfall runout zones on the basis of trajectory simulation data is generally missing. The objective of this study is to come up with standardized reach probability threshold values (RPTV ) to separate "realistic" from "unrealistic" simulated rockfall runouts. We therefore compared reach probability values (Preach) simulated with Rockyfor3D for 458 mapped, fresh rockfall blocks (silent witnesses SW ) on 18 different sites with a volume >= 0.05 m3 and estimated occurrence frequencies up to 300 years. We analysed which block, slope and forest characteristics 5 influenced Preach of the SW based on a linear mixed effects model. The results indicate that the limit of a realistic runout zone lies in the range where simulated Preach values are between >1% and approximately 3%. We conclude that RPTV can be defined to values lying in the range from 1.2% to 2.5% depending on the defined block volume and the encountered cumulative basal area in a forested transit zone. Where possible, the defined RPTV should be compared and validated by field recordings of SW . 10 keywords: rockfall hazard modelling; simulation, Rockyfor3D, silent witness, statistical analysis


Introduction
Throughout the world, at places where steep rocky slopes occur, the sudden and rapid downslope movement of single rock fragments, blocks or fragmented rock masses threatens human life and poses problems for infrastructure and residential areas.
An important basis for the implementation of above-mentioned types of measures is the rockfall hazard and risk assessment.
Although such assessments have been improved considerably over the last decades, it remains a challenge to fully grasp the unknowns of rockfall processes with respect to space and time . These are related to a range of factors which 25 include: precise release locations in a rock cliff (cf. Loye et al., 2009), the frequency and initial volume of the released rock mass (cf. Francioni et al., 2020;Farvacque et al., 2021;Hantz et al., 2021), the fragmentation, i.e. disaggregation of the initial volume during the failure process and breakage after the first and subsequent impacts (cf. Giacomini et al., 2009;Ruiz-Carulla et al., 2015;Matas et al., 2020;, the size and shape and change thereof of the individual rockfall fragments that propagate down the slope (c.f., Melzner et al., 2020), the transformation and dissipation of energy dur-30 ing rebounds on the slope surface (Caviezel et al., 2021) and impacts on standing or lying tree stems Lundström et al., 2009;Lu et al., 2020;Noël et al., 2021;Ringenbach et al., 2021), penetration depth in the soil and alteration of the terrain during impact (Pichler et al., 2005;Wang and Cavers, 2008;Guangcheng et al., 2015;Lu et al., 2019). To come up with realistic predictions of the areas that are potentially endangered by rockfall processes, modelling rockfall trajectories is one of the methods that provide an important information (Volkwein et al., 2011;Yan et al., 2020). To 35 account for some of the above mentioned uncertainties and to the intrinsic variability of rockfall as function of block shape, volume, exact initial position, etc., trajectory simulation models generally use stochastic variables in their algorithms (e.g., Bourrier et al., 2009). Since such probabilistic computational algorithms rely on repeated random sampling, also referred to as a Monte Carlo method (Von Neumann and Ulam, 1951), the numerical results are presented as probability distributions (Farvacque et al., 2020). In case of rockfall trajectory simulations, these generally include data on runout zones, kinetic energies, 40 and passing heights for all simulated individual rockfall blocks sometimes converted into individual risk estimates (Farvacque et al., 2019a).
Although attempts have been made to automatize the delineation of hazard zones on the basis of trajectory modelling (e.g., Jaboyedoff et al., 2005;Abbruzzese et al., 2009;Abbruzzese and Labiouse, 2020;Farvacque et al., 2020), in the daily practice this delineation is mostly based on human interpretation of the simulation results and subsequent definition of the realistic 45 runout zone for a given simulated scenario. Thereby, extreme long, low probability trajectories are separated from all other trajectories in the modelled rockfall runout distribution based on expert judgement, eventually supported by historical records, mapped deposited rocks (silent witnesses SW ), and other information recorded in the field (e.g., tree impacts). Legal considerations varying from one country to another also play a role, with often no clear correspondence between sound statistical and probabilistic concepts (Eckert et al., 2018). The basis provided by rockfall trajectory models for doing so is data on the number of passages per cell normalised by the total number of blocks potentially passing through a cell, which depends on the number of simulations per source cell and the number of feeding source cells. This normalised model output data can be referred to as reach probability data. The reach probability is a conditional probability depending on whether a block is released or not.
At present, no standardized reach probability threshold values (RP T V ) are defined for delimiting rockfall runout zones based on trajectory simulation data. This is usually done in a very subjective manner. Therefore, expert knowledge still represents a 55 great deal in rockfall hazard assessment.
A sound statistical comparison of simulation-based reach probability values with field mapped stopping locations of blocks from recent rockfall events could improve the quantitative basis for the delimiting rockfall runout zones. Consequently, the objective of this study is to come up with a quantitative basis for defining RP T V by comparing reach probabilities simulated with a three-dimensional rockfall model for mapped, fresh SW of 18 sites in Austria, France, Greece, Italy and Switzerland 60 for which the runout distance return period could be estimated as being maximum 300 years. In this paper we analyse which reach probability values are typically expected in the outer range of the rockfall runout zones and quantify which block, slope and forest characteristics influence the simulated reach probabilities of the mapped SW .

65
For this study, we used data from 18 different rockfall sites in Europe ( Fig. 1) with a total of 769 mapped SW (see table 1).
To analyse the effect of the complexity of the topography at the study site, we classified each study site in topography type 1 The mapped SW correspond to stopping locations of blocks mapped in the field. Only fresh blocks were taken into 70 account, meaning that blocks with weathered surfaces as well as those partly buried by material covering the surrounding slope surface were not recorded (Fig. 3). Preferable, we also detected additional rockfall marks such tree wounds and rockfall impact signs (craters) upslope from the mapped SW . The block volume (V ol) of each SW was calculated on the basis of the measured width, depth and height, as well as an estimated rounding factor varying from 0.52 (a perfect sphere) to 1 (a perfect rectangle). All data on the SW was gathered by the authors, except for the Gurtnellen site in Switzerland (CH), which was 75 provided by (Thali, 2009). At every site, we focused at mapping SW which, based on an expert interpretation (i.e., comparison with deposits in the direct surroundings, block volumes recorded in historical events databases of the region) corresponded to blocks that resulted from recent rockfall events with return periods of max. 300 years in terms of runout distance. Similarly, very frequent events were excluded from the analysis (cf., 2.3), so as to work with a sample of rockfall runout events that may well separate safe from unsafe locations in terms of hazard mapping for buildings and infrastructure.

80
Generally, we can distinguish 4 main types of mapped SW for the study sites (cf. table 1). At the Taesch (CH) study site, all deposited blocks of the rockfall event of August 2013 were mapped. At the Claro (CH), Flaesch (CH), Orvin (CH, c.f., Moos et al., 2018), Schmitten (CH) and Tithorea in Greece (GR) (, c.f., Saroglou et al., 2015) study sites, we mapped all de-posited rockfall blocks of multiple recent rockfall events. At six other study sites, we mapped selected freshly deposited blocks which had either a large V ol in comparison to surrounding blocks or a longer runout distance compared to the majority of the  Dorren et al., 2005b). Finally, at another six sites, we mapped selected blocks present at the study site which resulting from multiple rockfall events (e.g., Crolles in FR; cf., Farvacque et al., 2019b). We also mapped and recorded 90 field data required for the modelling (terrain roughness as well as soil types), following an intensive exchange on a common field mapping/recording method, corresponding to the Rockyfor3D manual (Dorren, 2016).  To characterise the forests at the study sites, we used an approach combining 1) mapping of forest stands based on orthophotos or, if available, high-resolution vegetation height models derived from airborne laser scanning data or drone flights and 2) 95 forest inventory plots in the field. At least one representative plot of 10 by 10 m or 20 x 20 m (depending on the tree density) was inventoried for each homogeneous forest stand type in the different study areas. From those inventory plots, we derived the required forest data for the used rockfall trajectory simulation model (cf. section 2.2).

Monte-Carlo based rockfall trajectory simulation model
The Monte-Carlo based rockfall simulation model used in this study is Rockyfor3D (for details see (Dorren, 2016)), which 100 is a probabilistic, process-based rockfall model simulating trajectories of individual blocks in three dimensions. Rockyfor3D was developed on the basis of full-scale rockfall experiments in the field and uses raster maps describing topography (Digital Elevation Model -DEM), rockfall source cells, the mechanical properties of the surface material and the slope surface roughness. To represent the forest, the model requires the number of trees, the mean and standard deviation of stem diameters at breast height (DBH), as well as tree types (coniferous or broadleaved) per cell as input data (Dorren et al., 2004.

105
For each rockfall source cell, the trajectories of a given number of rocks are simulated by considering flying and bouncing (rebounds on the surface). The density as well as the indentation resistance (cf. (Pichler et al., 2005), which have an effect on the penetration depth of the block during a impact before rebounding (which in turn affects the tangential energy loss of the block ) and dampening effect (represented by the normal coefficient of restitution Rn) of the impacted material is defined by the following eight soil types: Surface roughness is represented by three raster maps. These rasters represent deposited rocks and rock fragments covering 120 the slope surface, which form "obstacles" for the falling block. Micro topography (e.g., steps in the terrain) should not be taken into account here. The surface roughness was recorded in the field by identifying homogeneous zones in the study areas. These were represented as polygons on a map (in most cases printed hillshade maps) and later digitised and rasterised. Each polygon defines the surface roughness, expressed in the size of the material covering the slope's surface, looking in the downward direction of the slope. Three roughness probability classes need to be represented by a raster map and correspond to the height 125 of a representative obstacle in m that a falling block encounters in resp. 70%, 20%, and 10% of the cases during a rebound in the defined polygon. If the slope surface is smooth, a roughness value of 0 m was used. The choice of the roughness values needs careful attention, since Rockyfor3D is sensitive to this parameter.
During each simulated rebound on slopes with a gradient less than 30°, Rockyfor3D decreases the slope angle at the position of the rebound at random (following a uniform distribution), similar to (Pfeiffer and Bowen, 1989), up to a maximum decrease 130 of 4°. Rolling is represented by a sequence of short-distance rebounds with a distance in between that is equal to the radius of the block and an absolute minimum distance of 0.2 m. Rockyfor3D explicitly calculates the deviation and energy loss after impacts with trees dependent on the DBH, impact position, and the kinetic energy of the rock before the impact.
The main output of RockyFor3D used for this study consists of a raster map containing information on the reach probability (cf. section 2.3. Additional outputs generated by Rockyfor3D are, amongst others, raster maps with information on the kinetic 135 energies, the passing heights, the number of deposited rocks, the number of tree impacts per cell and the angle of the straight line between source and stopping cell (energy line angle ELA; see Dorren, 2016).
We simulated the exact V ol mapped for each site with Rockyfor3D using a rectangular block shape and extracted the simulated reach probabilities for the position of the deposited blocks. We simulated 100 blocks per start cell. At all sites, we used elevation models with a resolution of 2 x 2 m, except for the Greece study site, where only a 5 x 5 m resolution DEM was 140 available.

Reach probability
The reach probability (P reach) value in a given cell x indicates the probability (given in%) that cell x is reached by a block that has detached from the cliff. It is calculated by:

145
Where N B (x) is the number of blocks passed through cell x, N sims is the number of individual blocks simulated from each source cell and N sources (x) is the number of source cells "feeding" cell x. In other words, the reach probability is a measure of the number of simulated blocks passed through a given cell relative to the number of blocks potentially "feeding" the cell. It is typically used as indicator to determine the rockfall runout zone of a given release area. However, a quantitative basis for differentiating the realistic P reach values from the extreme ones (in terms of runout distance) which can therefore 150 be neglected for hazard mapping, is missing. In this study, we extracted for each SW the simulated P reach as the mean of the P reach values > 0 in the eight neighboring cells as well as in the center of a mapped SW . This value is hereafter referred to as P reach SW . We excluded SW with a V ol < 0.05 m 3 and a P reach SW > 5% from the originally 769 SW , since they were regarded as irrelevant for hazard analysis or are not in the outer reach, respectively. The threshold of 5% was determined based on a two-step outlier detection according to (Yang et al., 2019) using the median absolute deviation as score. In total, 155 202 SW from 9 sites were not reached by any simulated trajectory, i.e. their deposit location could not be reproduced by the rockfall simulation. All these SW had a significantly smaller V ol than the neighbouring SW which were perfectly reached by the simulations (Fig. 4). We therefore assumed that these SW are most likely fragments of larger blocks that broke off while falling down slope. On the basis of this assumption, we excluded those SW from the analysis. This resulted in 458 SW with a V ol > 0.05 m 3 and with a P reach SW > 0 and <= 5%.

Statistical analyses of the simulation results
We analysed which block, slope and forest characteristics influenced P reach SW of the SW that were reached by the simulations (n=458). We first tested whether there are significant differences in P reach SW between sites as well as between V ol classes based on a one-way Anova with a logarithmic transformation of P reach SW using the aov function of the stats package in the statistical software R. We then fitted a linear mixed effects model (lmm) for P reach SW with the variables reported in 165 was fitted to the log-transformed P reach SW values with stepwise backward variable selection with the aim to minimize the Akaike Information Criterion (AIC) using the step function (lmerTest package; (Kuznetsova et al., 2017)). We only included explanatory variables that were not substantially correlated (Spearman correlation coefficient < 0.4) to avoid colinearity. We 170 also tested for possible variable interactions, which did not substantially improve the model. The lmm was implemented with the lmer function of the lme4 package in the statistical software R (Bates et al., 2015). P-values were obtained by Wald-Chisq tests as well as likelihood ratio tests of the full model against the model without the variable in question. The model performance was assessed based on customary residual plots (Stahel, 2017) and the marginal and conditional R 2 (Nakagawa and Schielzeth, 2013). We further fitted a random forest model (rfm) (Breiman et al., 1984) and compared predicted values of the 175 lmm to the predicted values from the rfm. The latter was implemented using conditional inference trees as base learners in the cforest function of the party package in the software R (Strobl et al., 2009;Hothorn et al., 2010). The random forest model was fitted using three times repeated 10-fold cross-validation to reduce the risk of overfitting (Kohavi et al., 1995)).

Results
The mean P reach SW of the considered 458 SW was 1.79% and the median P reach SW was 1.41% (Fig. 5)  The Anova revealed a significant difference between sites, whereby the Tukey post-hoc test showed that only Claro and Taesch significantly differ from the others (Fig. 6), meaning that only the SW of these two sites have significantly different P reach values. Furthermore, the analyses showed that P reach SW is significantly higher for blocks with a V ol of 5-10 m 3 185 (Fig. 7).  According to both the lmm and the rfm, P reach SW is significantly influenced by the block volume class (V olClass), the normalized cumulative basal area (cbA), the mean slope (Slope mean ) and the soil roughness (rg70 maj ; cf. Table 3). The different soil types (soiltype maj ) were not significant in the lmm, but remained in the final model. P reach SW for a given SW increases with increasing V ol, whereas it decreases again for the largest V olClass (Figure 7), for increasing cumulative 190 basal area of the forest (Figure 8), increasing slope angle, decreasing soil roughness and decreasing dampening capacity of the soil. The lmm explains 45% of the variance (conditional R 2 ). There is a relatively high correspondence between the predicted values of the lmm and the rfm with slightly smaller values predicted by the lmm model (Fig. 9).    x-axis).

Discussion
The results of our analysis of a large number of deposited blocks from 18 sites in Europe imply that the lower limit of rockfall 195 runout zones typically have P reach values, simulated with RockyFor3D, between 1 and 3%. The statistical models showed that the P reach values, simulated with Rockyfor3D, of deposited blocks strongly depend on block, slope and forest characteristics.
Smaller P reach values are in particular achieved for blocks with a volume between 5 and 10 m 3 and densely wooded forest (cbA > 20 m 3 ). The analysis provides important information for the interpretation of simulation results for rockfall hazard mapping.

200
This study is conceptually based on the assumptions that the SW analyzed in this study were all registered in the typical outer range of a rockfall propagation zone, as it is typically done for hazard and risk assessment. However, we used a collection of study sites where different field protocols, especially with regards to the used criteria for SW recording, were applied. There are sites where only blocks with an "extreme" runout were recorded (e.g. Evolène) and others, where deposits also in the upper this study. We attempted to offset this by excluding all SW with a P reach > 5%.As already mentioned, an important bias is very probably the result of the recording of SW corresponding to block fragments (cf. fig. 4). Fragmentation is extremely common during rockfall processes and leads to difficulties when deciding on representative block shape and block size in the hazard analysis process (Jaboyedoff et al., 2005;Matas et al., 2020). In the field, it is rarely possible to differentiate which deposited blocks resulted from fragmentation in the rock cliff, or upon first impact below the cliff, or in the lower parts of the 210 transit area. This has been one of the arguments to keep all mapped SW in the original data set. It allowed us to profit from the large data source on different sites and slope conditions. The analysis of the volume distribution of SW that were not reached by the simulations finally clearly indicated that these SW were block fragments and, thus, we excluded them in the final data set.
The developed statistical model shows that P reach values increase with increasing V ol, increasing slope angle,decreasing 215 forest cover, and decreasing soil roughness. These factors promote an acceleration of the falling blocks leading to longer runout distances. Thus, we can conclude that on slopes with favorable characteristics for rockfall propagation, the probability that rocks in reality end up in the extreme range of the propagation zone increases. The fitted regression model yields an R 2 of 48% and the residual analysis, too, indicating that the model fits the data relatively well. The predicted P reach values of the regression model and the random forest model correspond well, suggesting that the two models are relatively robust.

220
The results indicate that the limit of a realistic runout zone lies in the range where simulated P reach values are between >1% and approx. 3% (= 70%-ile of all considered SW ). As a basic rule, based on the median and mean P reach values presented here as well as on long-term practical experience, we would recommend choosing a RP T V larger than 1% and smaller than 2% in a first step. If abundant data on field mapped SW allows for a detailed comparison with the simulated data, the RP T V can be eventually increased when supported by the field observations. Here, one must be sure, that blocks with extreme runout 225 distances were not removed in the field, as is often done by farmers (in agricultural fields) or by infrastructure managers (on road and railways). To eventually come with actual probabilities of a block having "extreme" runout, propagation probabilities have to be combined with the release frequencies of expert defined homogeneous release areas in a rock face.
When interpreting simulated P reach values, it is important to execute a sufficiently large number of simulations per start cell, to make sure that the P reach in the lower ranges of the runout zone converge. To do so, experience shows that at least 230 100 simulations per start cell are required. An analysis in the study site St. Paul de Varces (FR) showed that there is, however, no significant difference, even in the raster cells with low P reach values, between 100 or 1000 simulations per source cell.
For some areas, due to a topographic configuration leading to strongly diverging trajectories, it might be necessary to simulate 200 or in extreme cases 500 trajectories per start cell. Though for most model simulation based rockfall hazard analyses in the practice, 100 simulations per start cell will be sufficient. Furthermore, P reach depends on the definition of the source area and 235 can possibly be affected by the size, form and position of the source area. A known problem occurs when two or more source areas are superposed on a slope with a transit area between them. If some blocks from the upper source area fall over the rock face below this locally results in a higher number of sources underneath the rock face below, whereas the number of passages does not increase significantly. Finally this leads to lower P reach values compared to neighbouring cells which were only fed by the lower rock face. In such a case, it might be useful to simulate each of the superposed source areas separately.

240
To come up with more detailed categories of RP T V , we applied the developed lmm model to predict P reach values for diverse V ol as well as varying slope and forest characteristics. The derived RP T V categories are presented in Fig 10. As shown by table 3, the key data required for defining RP T V is finally the V ol and the cbA. For the latter, one can roughly differentiate between low, medium and high cbA.

Conclusions
On the basis of the presented results, we conclude that simulated P reach data are a valuable basis for delimiting a rockfall runout zone downslope from a release area. The limit of a realistic rockfall runout zone for rectangular rocks simulated with Rockyfor3D lies in the range where P reach values are between >1% and, depending on the V ol, slope roughness and the type and area of forest cover, smaller than 3%. As a basic rule, we therefore recommend choosing a RP T V larger than 1% 255 and depending on the before mentioned variables, RP T V can be defined in detail and fixed to values lying in the range from 1.2% to 2.5%. We recommend to delimit runout zones based on rockfall simulations only in combination with expert-based Table 1. Characteristics of the 18 study sites with mapped SW used in this study. Site provides the site name, which are abbreviated with the first 5 letters throughout the article, as well as the country (AT = Austria, CH = Switzerland, FR = France, GR = Greece, IT = Italy) and a sequential site number per country; Hc refers to the total cliff height (rockfall release or source area); φT r is the mean slope angle in the transit area; ST ype refers to the complexity of the topography explained in Fig. 2; LT r is the average length of the transit slope; G is the mean basal area in the forest between the SW and the foot of the cliff; n represents the number of mapped SW ; SW vol is the volume range of the mapped SW ; SW des described which SW we mapped in the field (1 = all blocks of one specific recent rockfall event, 2 = all blocks from multiple recent events, 3 = selected blocks (for explanation see text) of one specific recent rockfall event, 4 = selected blocks from multiple recent events). Site