Comparing Flow-R, Rockyfor3D and RAMMS to Rockfalls from the Mel de la Niva Mountain: A Benchmarking Exercise

Abstract

Rockfall simulations are often performed at various levels of detail depending on the required safety margins of rockfall-hazard-related assessments. As a pseudo benchmark, the simulation results from different models can be put side-by-side and compared with reconstructed rockfall trajectories, and mapped deposited block fragments from real events. This allows for assessing the objectivity, predictability, and sensitivity of the models. For this exercise, mapped data of past events from the Mel de la Niva site are used in this paper for a qualitative comparison with simulation results obtained from early calibration stages of the Flow-R 2.0.9, Rockyfor3D 5.2.15 and RAMMS::ROCKFALL 1.6.70 software. The large block fragments, reaching hundreds of megajoules during their fall, greatly exceed the rockfall energies of the empirical databases used for the development of most rockfall models. The comparison for this challenging site shows that the models could be improved and that combining the use of software programs with different behaviors could be a workaround in the interim. The findings also highlight the inconvenient importance of calibrating the simulations on a per-site basis from onsite observations. To complement this process, a back calculation tool is briefly described and provided. This work also emphasizes the need to better understand rockfall dynamics to help improve rebound models.

1. Introduction

The back analysis and trajectory reconstruction of rockfall events, such as those covered in the companion publication Noël et al. [1], provide valuable data for the important calibration of rockfall simulations (Jones et al., [2]; Labiouse, [3]; Berger and Dorren, [4]; Berger et al., [5]; Volkwein et al., [6]; Dorren et al., [7]; Valagussa et al., [8]; Bourrier et al., [9]; Noël et al., [10]). Simulations are often performed with various levels of detail with associated calibrations depending on the required level of certainty and safety margins. Since the required resources increase with the complexity of the required level of detail, rockfall hazards are often approached in different stages. As described in detail by Hantz et al. [11], Volkwein et al. [6], Jaboyedoff et al. [12] and Loup and Dorren [13], the typical stages where simulations are performed with an approximative increasing order of required certainty are as follows:

• Rockfall susceptibility/indicative hazard evaluations;
• Quantitative rockfall hazard zoning;
• Detailed risk assessment related to rockfalls;
• Mitigation measure designs.

First, the susceptibility stage is meant to quickly cover large regions with conservative safety margins. This ensures that only the areas, for which it is not possible to fall within any of the subsequent rockfall hazard zone classes, would be phased out at this stage. Geometrical methods are commonly used for delimiting the envelopes of areas that are potentially reachable by rockfalls at that stage, given their high efficiency to perform at a regional scale with grid-based/GIS-based approaches (Dorren, [14]; topographic–hydrologic approaches in Volkwein et al., [6]; Michoud et al., [15]). Such derived models can be expressed geometrically with the propagation defined from the reach angle, sometime called energy line or Fahrböschung, obtained from the expected horizontal and vertical distance travelled by the rocks. The alpha–beta model (Domaas, [16]), CONEFALL’s model (Jaboyedoff and Labiouse, [17]) and Flow-R’s simplified friction-limited model (Horton et al., [18]) are examples of such geometrical models. A rough estimation of the kinetic energies can also be obtained geometrically from these methods by considering the potential energy of the height difference (Δh) between the terrain elevation and the reach angle line, or the cone for a given rock mass and gravitational acceleration (g) (Jaboyedoff and Labiouse, [17]). Pseudo velocities (v_pseudo) without rotation can also be roughly estimated in the same way with Equation (1) (modified from Dorren et al. [7]) as follows:

$$v_{pseudo}=\sqrt{2g∆h},$$

(1)

However, these first geometric estimations should be distinguished from the process-based methods, sometimes called “dynamic” or “trajectographic” methods, often required in the subsequent stages (ONR 24810, [19]; OFEV, [20]; Mölk and Rieder, [21]; Direktoratet for byggkvalitet, [22]; NVE, [23]; Loup and Dorren [13]). With the constantly improving processing capabilities of computers, some have already transitioned toward using process-based methods at the susceptibility stage for potentially more realistic predictions with simulated behaviors closer to a rockfall’s modes of motion. For examples, Frattini et al. [24], Noël [25], Dupire et al. [26], Kalsnes et al. [27], Alvioli et al. [28], Alvioli et al. [29] and Dorren et al. [30] ran process-based simulations in batches with models derived from Pfeiffer and Bowen [31] to obtain billions of simulated rockfall trajectories over Norway, the Alps, and hundreds of kilometers of railways in Italy and Canada.

Given the large safety margins ensuring long runouts that would cover all subsequent hazard zone classes, the calibration process is rather simple at the susceptibility stage. As conceptualized in Figure 1, the simulated envelope should include most of the observed/mapped deposited rockfall blocks for the covered region. One should be careful at this stage to keep enough margin to account for the lack of per-site calibrations from field observations and an oversimplification of the simulated rockfall’s modes of motion depending on the simulation models used.

The quantitative rockfall hazard zoning and risk assessment of the second and third stages (E2 and E3 in OFEV, [20]) are well described by Jaboyedoff et al. [12], Volkwein et al. [6], Hantz et al. [11] and Loup and Dorren [13]. For a reference period, the spatial-temporal hazard probability is defined by the product of the probability of failure of a rock mass generating rock fragments within a given reference volume/energy range and the probabilities of the rock fragments reaching exposed objects while they propagate downslope. As covered by Jaboyedoff et al. [12] and Hantz et al. [11] and conceptualized in Figure 1d, it is important to choose the proper width of influence of the simulated trajectories combined with the exposed “reachable” objects or the proper related cell size of the gridded results since this has a great effect on the obtained reach probabilities. As the distance from a rockfall source increases, the propagation probabilities decrease not only because the rocks gradually come to rest but also because their lateral spreading gradually covers a wider area (Evans and Hungr, [32]; Jaboyedoff et al. [12]). It is therefore important that a model reproduces the mapped observed lateral deviations to obtain proper propagation and reach probabilities for quantitative use. Simulating a series of small samples of trajectories (e.g., 20–50) helps with performing the plausibility check by easing the visualization of each individual trajectory, as conceptualized with the yellow paths compared to the hypothetical mapped observed paths in red (Figure 1). A model producing overly channelized results, unless matching with the observations, should be fine-tuned further.

For rockfall hazard zoning, risk assessments or the design of mitigation measures, kinetic energies (intensities) are usually required to determine the vulnerability of the exposed infrastructure/building objects, and the population and capacity of the mitigation measures (Raetzo et al. [33]; Jaboyedoff et al. [12]; Volkwein et al. [6]; Lambert and Bourrier, [34]; Lambert et al. [35]). The simulations must then be calibrated at those stages to reproduce the observed kinetic energies. Additionally, since results can vary greatly from user to user (Berger et al. [36]; Berger and Dorren, [4]; Berger et al. [5]; Garcia, [37]; C2ROP, [38]), having real observations as a reference for calibration reduces potential user biases. Domaas [39], Descoeudres [40], Volkwein et al. [6], Wyllie [41] and Gerber [42] detail how analytical methods can be used to estimate the translational kinetic energies and bounce heights from the positions of a succession of impact marks and/or combined with the incident angle imprint in the ground at the impacts. These positions can also be determined on 3D detailed terrain models from the imprint impact marks, as in Caviezel et al. [43] and Noël et al. [1]. If video footage from previous rockfall events is available, these values can be reconstructed as well as the rock trajectory segments, as in Noël et al. [44] and Noël et al. [1], for example.

No matter what rockfall hazard stage is involved, it is important to be knowledgeable about the simulation models used, their predictability, sensitivity to different parameters, and objectivity (e.g., Lambert et al. [45]; Lambert and Bourrier, [34]). To evaluate these attributes and eventual limitations, the model’s results can be put side-by-side and compared with observed reconstructed rockfall paths and mapped deposited blocks from real events or rockfall experiments. For a relevant comparison in a practice context, equivalent time resources should ideally be spent per model and close to what is usually spent when delivering projects. With finite resources, the more time reserved to fine-tune the simulations from less predictable and slow software, the less available time there is for the preceding important data collection from field observations. Therefore, one must be careful when sharing perfectly matching simulation results to the rockfall community. These potentially site-dependent results can be misleading if the time spent for the precise calibration or how far the model needs to be pushed beyond its anticipated parameters to reproduce reality is not mentioned. Instead, one can better communicate the predictability, sensitivity, and objectivity of the simulation models by comparing their simulation results at an early calibration stage.

For this exercise, the historical events from Mel de la Niva Mountain, whose 3D trajectories and 2D paths were reconstructed in detail in the companion publication by Noël et al. [1], are used as references in this paper. The 200 m³ blocks, reaching hundreds of megajoules during their fall, greatly exceeded the common rockfall energies of the empirical databases used for the development of most rockfall models (Noël et al. [1]). Such a site can thus be challenging for the simulation models tested in this paper. For a pseudo benchmark with the historical events, preliminary simulation results from the Flow-R 2.0.9 simplified friction-limited geometrical model (Horton et al., [18]), Rockyfor3D 5.2.15 process-based model (Dorren, [46]) and RAMMS::ROCKFALL 1.6.70 process-based model (Leine et al. [47]) are qualitatively compared here. Similarities and discrepancies between the results and a comparison with the reconstructed trajectories are then discussed. As the simulations can be further fine-tuned, and the models can evolve and improve over time, the reader is encouraged to repeat and extend such back analysis calibrations and comparisons. The reconstructed paths and trajectories for this site are shared in CC BY 4.0 in Noël et al. [48], and swissAlti3D elevation models are openly and freely available from ©swisstopo. A back calculation tool is also briefly described following the benchmarking simulations and provided as Supplementary Material to help the reader repeat this comparison analysis on other sites.

2. Benchmarking Approach

For the simulations, the elevation terrain models (DTMs) were prepared with the 2013 DTM from the swissAlti3D product by ©swisstopo. It was resampled by average to the recommended cell size of 10 m for the simulations with Flow-R, and 2 m for Rockyfor3D and RAMMS. The parameters attributed to the terrain (Figure 2) for the three rounds of simulations are summarized in

Table 1. Parameters used for the three scenarios with each simulation software.

Flow-R 2.0.9
	1: Michoud et al. [15]				2: Scenario 2				3: Scenario 3
Reach Angle	33°				31°				25°
Pseudo velocity cutoff	∞				∞				30 m s−1
Rockyfor3D 5.2.15
Slope	1: Rapid automatic sim.				2: Rough				3: Rougher
[°]	Soil type	Rg70	Rg20	Rg10	Soil type	Rg70	Rg20	Rg10	Soil type	Rg70	Rg20	Rg10
0–15	1	0.01	0.01	0.01	1	0.01	0.05	0.1	1	0.01	0.05	0.1
15–25	2	0.01	0.01	0.01	3	0.3	0.5	2	3	0.3	0.5	2
25–35	3	0.05	0.05	0.05	3	0.3	0.5	2	3	0.5	1	3
35–45	4	0.05	0.05	0.1	4	1	2	4	4	1	2	4
45–52	4	0.05	0.05	0.1	5	0	0.05	0.1	5	0	0.05	0.1
52–90	6	0	0	0	5	0	0.05	0.1	5	0	0.05	0.1
RAMMS::ROCKFALL 1.6.70
Slope	1: Lu et al. [49]				2: Harder				3: Harder (ellipsoid)
[°]	Block 2 shape, 100 m3				Block 2 shape, 100 m3				Block 2 ellipsoid, 100 m3
0–15	Medium-hard				Medium				Medium
15–25	Medium-hard				Medium-hard				Medium-hard
25–35	Medium-hard				Hard				Hard
35–90	Medium-hard				Extra-hard				Extra-hard

. They are ordered from left to right following the order in which the simulations were performed. This corresponds to increasing induced runout distances with Flow-R and RAMMS, which both started with too short runouts. For Rockyfor3D, the first simulation was performed with the rapid automatic simulation parameters and was followed by rougher surfaces to reduce the induced runouts.

2.1. Flow-R 2.0.9 Methodology

The parameters of the first scenario with Flow-R correspond to the finalized fine-tuned ones used for a regional analysis in a nearby valley sharing similar terrain features (Michoud et al. [15]). The second scenario uses the same parameters but with a desired reach angle (“basal flow friction”) of 31° that corresponds to the H/L value of 0.6 in Corominas [50], and to reach angles of rockfalls from mountainsides with heights from 100 m and up to 350 m in Domaas [16]. The last scenario uses a reach angle of 25° since the reconstructed trajectories show that the blocks had relatively stable velocities for slopes around that angle (Noël et al. [1]). A pseudo velocity cutoff of 30 m s⁻¹ was added to the last scenario to restrict its runouts that otherwise tend to follow the downstream path of the river at the bottom of the valley. All 2D simulations with this software were performed as in Michoud et al. [15] using the original diffusive flow direction by the gridded groundwater hydrological runoff model from Holmgren [51] with an exponent of 1 that produces highly diffusive flow paths equivalent to that of Quinn et al. [52]. The clever addition of the height factor to emulate debris flow thickness for a diffusion that can pass over small obstacles by Horton et al. [18] was kept to 0 m. The diffuse flow pathways were refined by using the gridded debris flow model from Gamma [53] implemented in Flow-R. As this model is not based on rockfall dynamic processes, there are no parameters related to the rock shapes, sizes or volumes.

For a visual comparison of simulated paths from Flow-R with the reconstructed trajectories from the Mel de la Niva events (Noël et al. [1]), we produced overlaying pseudo paths and deposited blocks based on Horton et al. [18]. Five paths per scenario were simulated using a random walk approach following the probabilistic flow direction distributions used by Flow-R. With this approach, a high number of random walk paths would converge to the same results as those produced with Flow-R. The same corresponding parameters as the original Flow-R simulations were used to obtain the same runout area envelopes. Preparing the input layers, adjusting the parameters and running the different Flow-R simulations were performed in less than an hour thanks to the simplicity and efficiency of the model.

2.2. Rockyfor3D 5.2.15 Methodology

For Rockyfor3D, preparing the terrain model and the five corresponding source layers took approximately 15 min. The four other required layers were generated based on the slope thresholds used by the rapid automatic simulations and the roughness of the site detailed in Noël et al. [10]. The open-source automation tools developed for assessing the rockfall potential along 260 km of railway during the ParaChute project (Noël, [25]; Cloutier et al. [54]; Cloutier et al. [55]) were used to automatically generate the terrain material layers for the three scenarios in approximately 15 min. The d₁-d₂-d₃ dimensions of block 2, 11.1-8.5-4.0 m, were used with an ellipsoid shape by Rockyfor3D, which estimates the inertia value for the simulations. Each scenario was run two times, once with 10,000 trajectories to produce the runout area envelopes and once again to randomly extract five trajectories to be compared to the observed paths for an insight into the simulated propagation and reach probabilities. The whole process was completed in less than an hour thanks to the automation tools and the efficiency of the model.

2.3. RAMMS::ROCKFALL 1.6.70 Methodology

For RAMMS, block 2 was first converted to a .pts file to be used as a real input shape by the software. The 3D model of the block, aligned on its principal axes of inertia, was subsampled by space to a total of 72 points, which is slightly higher than the included shapes with the software. It was then built as a custom block in RAMMS while limited to 263 metric tons due to software restrictions of 100 m³. This block was used for the two first scenarios and simplified to a corresponding ellipsoid for the last scenario to see the effect of the simplified shape. The block was converted to an ellipsoid with the same d₁-d₂-d₃ dimensions, subsampled by space to a total of 82 points and converted to a .pts file for the last scenario. These conversions and building custom blocks in RAMMS took approximately half an hour.

For the first scenario, the same terrain parameters previously used by the WSL Institute for Snow and Avalanche Research SLF to roughly simulate the 2015 Mel de la Niva rockfall event (Lu et al. [49]) were used. They were combined with no initial velocity and an automatic offset to the source cliff. Then, because most simulated blocks stopped a few hundred meters from the source, the terrain parameters were changed to values based on the terrain types described in the user manual: matching scree slopes and mountain meadows with their corresponding parameter. This produced short runout results similar to those of the first scenario, so this was not kept as a second scenario. Then, more realistic runouts were obtained by giving an initial horizontal velocity of 10 m s⁻¹. This could be partly justified by the toppling behavior of the failure but was discarded, as the other software did not require these additional adjustments of their default initial conditions. The hardness of the terrain material was raised instead, setting the scree to the parameter for hard bedrock and the mountain meadow to the parameter for hard rock scree. This helped increase the runouts for the last two scenarios. Surprisingly, a 100 m³ non-flat block (real_equant_1.2.pts) with the same parameters as the two last scenarios significantly increases the runouts (see the simulation results dataset in Noël et al. [48]). For the last scenario, only the shape of the block was changed to the corresponding ellipsoid based on the dimensions of block 2, keeping the same hard material surfaces.

Each scenario was run two times, once with 1000 trajectories to produce the runout area envelopes and once again to randomly extract five trajectories to be compared to the observed paths for an insight into the simulated propagation and to reach probabilities. The maximum number of simulations per scenario had to be reduced from 10,000 to 1000 due to excessive processing time, as more than one night was needed to simulate 10,000 trajectories for this large site while keeping the time step and “dump” step to 0.01 and 0.02 s, respectively. The RAMMS simulations were thus completed in approximately one day, including the fine-tuning work but not counting the first attempt for 10,000 simulations. Without considering license costs, the reader should bear in mind that for a fair comparison in terms of resources, the time spent with RAMMS should have been limited to the same time spent as the others, or the other way around. Also, the simulations should have been performed with the same computing hardware. However, this was not feasible given the license availability.

2.4. Localized Reach Probability and Hazard

To complement the raw simulated results for the comparison, they are also expressed in terms of quantified localized reach probability and guiding hazard zones. The results of the second scenarios are used for that purpose. Note that being expressed in those terms does not improve the preliminary nature of the simulated results limited by their early calibration stage. This is only done here for comparison purposes to help highlight by how much similar simulations would miss the natural phenomena if one overlooks the important fine-tunning steps from field observations. For a relative comparison from one to another, the temporal aspect must be the same for the three software programs. A finalized temporal estimation for the site is, however, not required for this comparison exercise from early calibrated simulations. Given the mass of the reference block scenario of 500 metric tons, the energies involved when the blocks are in movement (at >1.1 m s⁻¹) fall in the highest intensity level of the Swiss danger zonation matrix (combining hazard on the x axis to the corresponding intensity on the y axis). Thus, no emphasis is made related to the intensities or protective effect of the forest for this relative comparison.

For a punctual source, the probability of an exposed object to be reached by a rockfall can be evaluated by summing up the number of unique trajectories that intercept the object and exceed the intensity level at which the object is vulnerable. As mentioned in the introduction (Figure 1d), the dimensions of the object and those of the simulated rock fragment must be considered (Jaboyedoff et al. [12]; Corominas et al. [56]; Hantz et al. [11]). The local reach probability of the object (P_r) is then obtained by dividing the obtained number of intercepting trajectories by the total number of simulated trajectories from the given source.

For simulation models considering the fragmental aspect, a unique trajectory must include the branching segments of its fragments (e.g., HY-STONE from Agliardi and Crosta, [57]; Crosta and Agliardi, [58]; Frattini et al. [24]; EG4 Risk, [59] or RockGIS from Matas et al., [60]). This way, numerous fragments intercepting the exposed object only count as one occurrence if they are related to a common unique trajectory. In that case, the temporal aspect (λ) should be attributed to an expected yearly frequency of failure events from the source (before fragmentation).

The temporal aspect can also be expressed as a yearly frequency of block fragments, for example, estimated from deposited block inventories (Evans and Hungr, [32]) or mapped past trajectory segments (after fragmentation). In that case, the number of trajectories intercepting the exposed object should be estimated from unique trajectories individualized for each simulated fragment. The parent segments starting from the source should be duplicated until their number matches with the number of fragments.

The latter situation can be approximated with simulation models that do not account for the fragmental aspect, like for the software tested. Each simulated trajectory corresponds to a simulated fragment. The effect of this simplification is marginal when the fragmentation is estimated to happen in the vicinity of the source. However, one must remember that the fragmentation induces changes in shape, size and volume of the propagating blocks. Thus, their initial propagating characteristics and related intensities are different than those of the final fragments, which is not accounted for by the simplified simulations.

For an object exposed to numerous (diffuse) sources, rather than evaluating its local reach probability (P_r), it is more appropriate to estimate its local average expected frequency (λ_object) of being intercepted by trajectories with exceeding intensities. For that purpose, a fraction of λ must be attributed to each unique trajectory, ideally expressed as a ratio of the corresponding source surface area covered (not the 2D planimetric surface). Additionally, this ratio has the positive side effect of compensating for the strong biases associated with planimetric gridded source distributions (Figure 3 in Noël et al. [10]). For example, if 10,000 trajectories are simulated from a 1000 m² cliff from 100 sources evenly spread (without grid bias for simplicity), each trajectory corresponds to 0.1 m², and each source covers 10 m². Assuming a rockfall event produces 20 fragments in average for that site, an expected yearly λ of 0.2 failure events for that cliff would correspond to a yearly λ of 4 fragments per 1000 m² cliff in average. Therefore, each source would carry a yearly λ of 0.04 fragments for its 10 m², which would be subdivided to its 100 related trajectories. An exposed object intercepting 400 unique trajectories with exceeding intensities would thus be related to 40 m² of cliff and have an expected average λ_object of 0.16 fragments (or 0.008 failure events if the trajectories were obtained with a fragmental simulation model). One could develop this further by considering variating frequencies of release function of the likelihood of failure attributed to different segments of the cliff (e.g., overhangs, unfavorable orientation related to joint sets, unfavorable presence of roots and trees with jacking effect, unfavorable water table conditions, variating rock quality, etc.).

Finally, the local hazard probability for the object to be reached at least once (P_Hazard) during a given period (t) can be obtained with Equation (2) for a punctual source or Equation (3) for multiple sources, modified after Jaboyedoff and Labiouse [17], Abbruzzese and Labiouse [61], and Hantz et al. [11]:

$$P_{Hazard}=1-{\left(1-P_r\right)}^{\lambda t},$$

(2)

$$P_{Hazard}\approx 1-e^{-{\lambda }_{object}t},$$

(3)

If evaluated from finalized fine-tuned simulations, it is recommended to express the results (e.g., guiding zonation limits based on regulated hazard thresholds) for the boundaries of exposed objects rather than their centers. A practitioner could then draw their desired zones by expert judgment while being guided by the simulated guiding zonation limits and field observations.

With RAMMS’s results being exported in vectorial format (in addition to the gridded results), no workaround was required for producing the preliminary reach probability and guiding hazard zones (from the early calibrated scenarios). The methodology for expressing RAMMS’s results in those terms is thus described first, followed by the one of the workarounds required with Rockyfor3D and Flow-R. Following the recommendations of Jaboyedoff et al. [12], a 30 m wide hypothetical exposed object was used in combination with the maximal diameter (d₁) of the simulated block. It was represented by a circular moving window with a diameter of 30 m + d₁ used to evaluate the number of intercepting trajectories in every location. The obtained values were then divided by the total number of simulated trajectories to obtain the reach probabilities from RAMMS’s results. A yearly frequency of 13/40 block fragments estimated from the 13 mapped rockfall paths observed over a period of about 40 years (Noël et al. [1]) was used for the local hazard probability with Equation (2). The values of the moving window being saved in its center, the results are expressed for the center of exposed objects.

Counting the number of intercepting trajectories is complexified with Flow-R and Rockyfor3D because of their rasterized gridded results. Moreover, susceptibilities are obtained with Flow-R from its hydrologically diffused “rockfall” quantities. For a comparison with the process-based results of the other software programs, the susceptibilities were considered as pseudo propagation probabilities. Normally, propagation probabilities are obtained as gridded results by dividing the number of passing trajectories with the total number of simulated trajectories. This relation was used to convert the pseudo propagation probabilities to pseudo passing trajectories assuming a pseudo total number of simulated trajectories of 10,000.

The output width (d_{traj_out}) of the 10 m wide pseudo passing trajectories from Flow-R and 2 m wide passing trajectories from Rockyfor3D are defined by their respective cell size. As for RAMMS, the hypothetical exposed object was represented by a circular moving window with a diameter of 30 m + d₁. Inspired from previous studies (Jaboyedoff et al. [12]; Lan et al., [62]), the passing trajectories inside the window, including the cells with zero passing trajectories, were averaged (n_{pass_avg}) using the SAGA Focal Statistics or GRASS r.neighbors tools in QGIS [63]. For a better average with Flow-R, its pseudo passing trajectories were temporarily upscaled to a cell size of 2 m via bilinear interpolation for this averaging step. Assuming the paths were mostly parallel, the averaged values were then converted into an estimated number of intercepting trajectories (n_reach) for a pseudo width of 30 m + d₁ using Equation (4) as follows:

$$n_{reach}\approx n_{pass\_avg}\frac{\left(d_{object}+d_1\right)}{d_{traj\_out}},$$

(4)

The approximated corresponding local reach probabilities (P_r) from a single source and local hazard probabilities (P_Hazard) are finally obtained as before by dividing the number of intercepting trajectories by the total number of simulated trajectories and with Equation (2). Complementary details about this method for quantitative guiding zonation and a related freeware tool are described in Appendix A.

3. Simulation Results and Analysis

In this section, the simulated results are presented in the form of four figures (Figure 3, Figure 4, Figure 5 and Figure 6) where the models are put side-by-side for an easy qualitative visual comparison. They are overlayed by the mapped rockfall paths, with deposited block fragments shown in red to follow the style and principles of the conceptual Figure 1. As a reminder, the reader can find detailed information about the compared mapped observations in the companion publication (Noël et al. [1]). As for any mapped datasets, one should bear in mind that mapped observations are never fully complete. Deposits and paths from past events are slowly “erased” by the fines and morphological modifications from other processes. Also, the origin of the mapped deposited block fragments does not necessarily correspond to the same source as the one of the simulations and some may be glacial erratics erroneously mapped as rockfall blocks. However, 11 of the 13 mapped historical paths converge toward an origin in the vicinity of the simulated source. The closer the simulation results are to the observations, the more confident one can be about the model for that site (Dorren et al. [7]). The four figures are designed based on Figure 1 to help answer the following related guiding questions that should be explored when validating simulation models. The reasoning steps for validation and their associated questions are interspersed with the following result figures.

Figure 3 and Figure 4 are designed for evaluating the sensitivity and predictability of the models from the simulations performed within an anticipated realistic range of parameter values and for observing how far the preliminary results vary from reality. For this validation step, the following questions can be explored:

• By how much do the simulations fail/succeed at reproducing the real events and at covering the mapped deposited rocks?
• How far/close do the parameter values need to be pushed beyond the previously anticipated realistic range of values for reproducing real events, and how does that translate when back analyses are not possible for a given site? In other words, is the model predictable, or should the previous range of parameter values be extended by precaution when back analyses are not possible?
• Do the results exaggerate the runout distances, involved energies and bounce heights leading to safe but costly/constraining susceptibility and hazard zoning, risk assessments and mitigation designs?
• Contrarily, can the results potentially lead to consequences in case of any eventuality due to underestimated simulated values?

Combined with Figure 5, Figure 3 and Figure 4 are also designed to evaluate the simulated lateral deviations, overall distribution, and values to the mapped and reconstructed rockfall data from a small sample set of paths (e.g., 15–50). For this validation step, the following questions can be explored:

• Do the simulated paths unrealistically oscillate laterally compared to the observations, potentially leading to unrealistically undertaken paths and erroneous associated reach probabilities (e.g., the yellow vs. red paths in Figure 1a,b)?
• Do the small sets of simulated paths and their distribution fall close to the observations for similar paths so that the simulations are statistically representative (e.g., the yellow vs. red paths in Figure 1a,c)?
• Do the vaults or curvature of the simulated freefalling parabola segments and separating distances in between the impacts compare well with the observations so that the bounce heights are realistic for use of the simulation with mitigation design tasks?
• Do the simulated velocities or related energies compare well with those from the reconstructed trajectories so that intensities based on preset building vulnerabilities can be estimated for land planning zonation tasks, risk assessment studies or mitigation design tasks?
• Are the simulated velocities, bounce heights and deposited rocks matching with the reconstructed data? In other words, is it possible to obtain all those characteristics at once, or must simulations with different parameter values be performed in parallel for each desired piece of matching information?

Combined with Figure 6, Figure 5 is also designed to qualitatively evaluate if the area covered by the simulations is statistically representative of a high number of simulated paths (and/or a number equivalent to what may be used when delivering projects). For this validation step, the following questions can be explored:

• Does the simulated envelope cover most of the mapped deposited rocks suspected to come from the same source(s) so that the results could be used for guiding susceptibility/indicative hazard zoning (e.g., the yellow vs. red dotted envelopes in Figure 1a–c)?
• Would an additional buffer around the simulated paths and/or different parameter values be needed to increase the area covered by the envelope to catch most of the mapped deposited blocks?
• Are the simulated paths and deposited rocks overly channelized compared to the mapped observations, leading to unrealistic statistics for the reach probabilities (e.g., the yellow vs. red paths and blocks in Figure 1a,b)?
• Contrarily, are the simulated rocks deposited near mapped deposited ones in proportional numbers so that the simulations are statistically representative for quantitatively guiding hazard land use zonation and risk-assessment-related tasks (e.g., the yellow vs. red paths and blocks in Figure 1a,c)?
• Are the path widths properly considered to obtain statistically representative reach probabilities from the propagation probabilities (e.g., see the width influence on the reach probabilities in Figure 1d)?

The results of each software are analyzed individually in the following subsections for objectivity. As a reminder, the volumes and energies involved for that site greatly exceed those of the empirical data used for the development of the process-based simulation software used in this comparison (companion publication Noël et al. [1]). It can thus be challenging for the models to reproduce the rockfall events of that site.

3.1. Flow-R 2.0.9 Analysis

Concerning the first scenario with Flow-R (Figure 3a and Figure 4a) reusing the fine-tuned parameters from Michoud et al. [15] in yellow, many observed paths and deposited blocks (in red) are surprisingly not covered by the 2D gridded simulations despite the similarity of the neighboring valleys covered in Michoud et al. [15]. In fact, most observed blocks that deviate toward their left side when looking downslope are not reached by the simulated envelopes of the three scenarios despite using an exponent of 1 for a high lateral diffusion of the flow equivalent to the one from Quinn et al. [52]. Indeed, the maximal induced lateral diffusion of the simulated flow is not high enough to reproduce the rockfall paths observed for this site. However, increasing the diffusivity by using a smaller exponent value outside of its designed range [1 to ∞] (Horton et al. [18]), which goes against the original intended use of the exponent (Holmgren, [51]), would have the adverse effect of lowering the runouts. Indeed, Flow-R is designed to stop the flow locally when its pseudo quantity, which becomes “diluted” with increasing diffusion, is not sufficient for sustaining the flow propagation, i.e., when the local susceptibility falls below a certain threshold (Horton et al. [18]). When experimenting with a smaller exponent, one should thus lower this threshold.

The lack of lateral spreading is also visible for the 15 extracted paths from the random walk that mostly follow the steepest downslope gradient with local sharp turns due to their eight allowed directions toward the neighboring cells in the raster grid (N, NE, E, SE, S, SW, W and NW) (Figure 3a and Figure 5a). The 2D gridded simulated diffuse flows follow the topography instead of reproducing the observed rockfall paths since 2D hydrological runoff models generally follow the steepest path (Volkwein et al. [6]).

The runouts in directions other than the eight neighboring cells are shortened due to an increased number of sharp turns and small oscillating detours via the neighboring cells (Dorren, [14]). As a result of this grid bias, the simulated envelopes from Flow-R barely managed to reach the CONEFALL reach angle contours of the corresponding parameters used (33°, 31° and 25°), even without a pseudo velocity cutoff (Figure 3a, cumulative curves in Figure 4a and Figure 7). When drawing maximum runout limits based on reach angle thresholds with Flow-R, one should interpret its results with caution because this grid bias affects the predictability of the model. Additionally, one should account for such bias by using high pseudo velocity cutoff values and lower reach angles than the desired ones from empirical rockfall databases. For rockfall modeling, lowering the reach angle by approximately two degrees ensures all areas that would be reached without the grid bias are reached (Figure 7c). The reach angle can be precisely adjusted geometrically (Figure 7b) with Equation (5) as follows:

$$\tan {\phi }_{Adjusted}=\cos \frac{45°}{2} \tan \phi ,$$

(5)

Future versions of the model could improve on this aspect by considering the second and third 16 to 24 cell neighbors, higher resolution grids and recalibrated the susceptibility threshold as a workaround.

Nevertheless, the reachable envelopes from Flow-R are the easiest to read with a clear distinction from one scenario to the others compared to those of the other software (Figure 3). This can help draw susceptibility/indicative hazard limits. The simulated results were also easy to quickly produce thanks to the efficiency of Flow-R. Therefore, for sites where the expected rockfalls are channelized and do not deviate laterally as much, Flow-R can be an option for the estimation of the area susceptible to rockfalls if low reach angle and high pseudo velocity cutoff parameters are used.

In the case of the pseudo velocities (Figure 4), they differ from those of the reconstructed trajectories with their smooth velocity transitions. Indeed, they lack the typical rockfall sawtooth shapes of the reconstructed trajectory velocities. This is due to the 2D gridded geometrical reach angle approach of Flow-R based on diffusive ground water hydrological runoff and debris flow models (Holmgren, [51]; Gamma, [53]). The applied pseudo velocity cutoff of the third scenario in blue is clearly visible in Figure 4a with the horizontal velocity segments at 30 m s⁻¹. When ignoring the lack of sawtooth shapes from Flow-R’s results, the decreasing velocities of the second scenario match well with the last 1000 m of block 2.

With the questionable hypothetical situation where quantitative rockfall hazard zoning would not require complementary process-based rockfall simulations, one should interpret the simulated results from Flow-R with caution if solely used for that purpose (Figure 5a,b and Figure 6a). Indeed, the 15 extracted paths differ greatly from the historical paths (yellow, green, and blue paths vs. red paths in Figure 3a). However, where in reach by the simulations, the susceptibilities from scenario 2 directly outputted by Flow-R are generally higher than the propagation probabilities from Rockyfor3D and RAMMS, which are safe (Figure 5a,c,e). However, this difference disappears once converted to “reach probabilities” as an intermediate step toward hazard zoning (Figure 5b,d,f). Also, this step does not increase the footprint of lower probability areas in that case (Figure 6a) in comparison to the other process-based results (Figure 6b,c).

3.2. Rockyfor3D 5.2.15 Analysis

Concerning the first scenario with Rockyfor3D (Figure 3b), its blue envelope safely manages to catch most historical paths. It is interesting that this first simulated scenario covers all the observed historical paths using the rapid automatic simulation parameters. Indeed, this preliminary result lies on the safe side, while subsequent results can become increasingly less conservative by reducing the safety margins as the knowledge of the site increases is a safe and effective approach. The envelopes of the first and second scenarios also reproduce the curved end of the longest runout in the lower right corner (the blue envelope plus one thin green trajectory envelope from the second scenario). The envelope of the second scenario is reasonably pulled back, using roughness dissipating (damping) parameters close to the observed roughness in the field. The addition by Dorren [65] of such measurable dissipating parameters over the “traditional” dissipating coefficients (R_N and R_T) from Pfeiffer and Bowen [31] helps improve the objectivity of the rebound model.

When looking at the cumulative profiles with the previously stated simplifying assumption that the mapped blocks share a common origin (Figure 4b), the rapid automatic simulations in blue safely overestimate the runouts by slightly more than one degree beyond a reach angle of 32.5°. This is the closest to the cumulative curves of the mapped blocks in the most critical part from 80% to 100% of the cumulated blocks of the three software programs tested. Above 20% of the cumulated blocks, the second scenario in green succeeds at closely following the curve of the mapped blocks (all n = 449) and of the 11 blocks with mapped paths that converge to a similar source origin. However, it is slightly shifted to the left around the most critical area, marginally underestimating the maximum runouts as a result. Further fine-tuning of the parameters toward obtaining final simulation results for that case should thus focus on finding values in between the first and second scenarios. From the cumulative profiles, it is interesting to note the large effect of adding 0.2 m, 0.5 m, and 1.0 m to the Rg70, Rg20 and Rg10 roughness of the 25–35° slopes of the third scenario in comparison to the second scenario. These are the only differences between the parameters of the two scenarios and seem to significantly affect the runouts of the 11.1 m block (d₁).

Concerning the lateral deviations with Rockyfor3D, the simulated trajectories are marginally less channelized than with Flow-R and form curved paths slightly closer to those observed (Figure 3 and Figure 5). Nevertheless, slight lateral oscillation persists due to the probabilistic deviations applied by Rockyfor3D’s rebound model at impact. However, the deviations are not applied during the freefalling phases. Indeed, the simulated blocks pursue their strait parabolic course without being affected by the obstacles they fly over until their next impact. This allows the blocks to resist deviating toward the steepest downslope direction during freefall, reducing the channelization as a result. However, the 15 extracted simulated trajectories are too channelized to reproduce the 13 historical mapped paths at that site.

Drawing the susceptibility/indicative hazard limits would be facilitated by adding a buffer around the 2 m thin simulated paths due to their small footprint limited to the width of the cell size. Alternatively, this could also be done from a quantitative hazard approach as described in Section 2.4 (Figure 6) using conservative simulated runout distances or a reach probability threshold for the susceptibility contours. Indeed, a pseudo buffer of the width of the block combined with the exposed object is added to pass from the propagation probabilities (Figure 5c) to the reach probabilities (Figure 5d). However, the second scenario is overly channelized, and many mapped deposited block fragments would not be covered unless a very low reach or hazard probability were to be chosen in that case. Such very low probabilities can, however, be obtained easily by simulating a high number of trajectories thanks to the high efficiency of Rockyfor3D.

For the process-based translational velocities (Figure 4b), it is possible to recognize typical rockfall sawtooth velocity profiles on the simulated trajectories. Additionally, most of them stay under 43 m s⁻¹ with what strangely looks like an undocumented velocity cutoff. A few points go beyond this velocity, however. Further investigations on the velocity cutoff can be found in Appendix B. While finalizing this manuscript, the authors received confirmation that there is indeed a velocity cutoff implemented in Rockyfor3D and that its description will be added to the transparent description of the model (Dorren, [46]). Despite this strange behavior, the simulated velocities from the first scenario with the rapid automatic simulation parameters oscillate in a sawtooth pattern near the velocities from the reconstructed trajectories. Those from the second scenario with a roughness similar to that observed onsite do not increase sufficiently after 800 m to match with those from block 2.

Similar to Flow-R, Rockyfor3D’s results were quickly produced thanks to the efficiency of the software. This leaves more resources for the important on-field data collection aspect and facilitate the fine-tuning and calibration processes. Rockyfor3D’s process-based results should also be interpreted with caution with the tested version if used for quantitative rockfall hazard zoning of such challenging sites with significant lateral spreading. If Rockyfor3D’s good predictability and objectivity remain valid when applied to sites involving more channelized rockfall paths, Rockyfor3D should be a good option for safe and objective quantitative hazard zoning. With limited historical paths available or limited resources for extensive back analyses from field observations, Rockyfor3D is a good option for susceptibility mapping out of the three tested software programs with its rapid automatic simulation parameters. Indeed, it has a similar efficiency to Flow-R but with simulated behaviors closer to rockfall dynamics. Additionally, based on the simulation results limited by their early calibration stage for this site, Rockyfor3D seems to have a more objective and safer approach related to the prediction of the runout distances compared to the tested version of RAMMS::ROCKFALL analyzed hereunder.

3.3. RAMMS::ROCKFALL 1.6.70 Analysis

For RAMMS’s results, the envelopes obtained are 10× sparser than those from Rockyfor3D due to the lower number of simulated trajectories despite the considerable time spent on this software compared to the others. As with Rockyfor3D, a buffer should be added to increase the width of the rasterized trajectories because of their small footprint limited to the width of the cell size (2 m). Surprisingly, the parameters of the first scenario, which were also previously used for that same rockfall event (Lu et al. [49]), did not produce the expected corresponding runouts (Figure 3c and Figure 4c). Indeed, as shown by the five yellow trajectories, such simulations struggle to go beyond the 37° reach angle contour. Four out of five trajectories did not go further than 350 m from the source, and two stopped at the foot of the cliff. The cumulative curve also shows that nearly 60% of the simulated blocks for the first scenario stopped with reach angles steeper than 42°, considerably underestimating the runout distances (Figure 4c). Similar results were also obtained using parameters corresponding to the terrain materials encountered on the site. Instead, the two other scenarios with more realistic boosted runouts were obtained by setting the material properties to values further away from what the authors would have anticipated at that site or from what is recommended in the manual of RAMMS for large blocks. The Extra-hard material was applied for most of the first 600 m and to set to Hard for most of the rest of the slopes (Figure 2).

In this fine-tuning calibration process, the results with underestimated runouts were first obtained and gradually moved toward realistic runout values while moving away from parameters that were first thought to correspond with the terrain material properties. Such an approach for the runout distances follows the opposite direction than the previous safe approach from Rockyfor3D. Indeed, time resources must be spent during this fine-tuning process to increase the area susceptible to rockfalls toward realistic values. Overlooking this back analysis process from field observations in that case would have the serious consequence of predicting too short runouts and small areas susceptible to rockfalls. The objectivity and predictability of the model thus seems lower than that of Rockyfor3D for the runout distances given that significantly underestimated runouts were obtained when initially using parameters close to the site properties or those used by others for that same site. Further fine-tuning of the parameters toward obtaining final simulation results for that challenging case should use parameters increasing the runout distances, thus moving away from the anticipated values. The software limitation on the block volume to 100 m³ may have played a role in underestimating the runout distances.

This difficulty in simulating longer runout distances with moderately to very flat blocks might also be caused by a known issue related to the rotation stability of the block shapes with RAMMS. Indeed, the block’s rotations are erroneously not stable around their d₁ or d₃ axes with the current model, as recently highlighted by Leine et al. [66], and thus struggle to reproduce a “wheel-like” behavior such as the one observed for the block fragments of the 2015 event of the Mel de la Niva site (Noël et al. [1]). Knowing this issue, a correction for the future released versions of RAMMS::ROCKFALL is being implemented. Finalizing this manuscript, the authors received a preview given by the SLF of a prototype of the future corrected version of RAMMS. With the prototype, they tested the first scenario applied in the same way as in this manuscript but with a 200 m³ block 2. The results of the preview seem promising both in terms of predictability and performance (processing time). This will of course remain to be confirmed once implemented in a released version of RAMMS. Meanwhile, the block shapes should be used with caution with the tested version of RAMMS::ROCKFALL when applied for similar challenging sites and volumes far from what was used for the development of the model. At present, the results produced with the latest version available (v.1.7.60), released in October 2022, are similar to those produced with version v.1.6.70 and were produced in a similar amount of time (see the simulation results in the dataset from Noël et al. [48]).

Nevertheless, with the boosted runouts of the second and third scenarios in green and blue (Figure 3c), their extracted trajectories and those visible on the envelopes in the background have the most realistic paths of the three software programs tested. Indeed, they are often parallel to the observed historical red paths, they reproduce similar curvature and lateral deviations, and they are not overly channelized. The envelope of the second scenario with the real shape of block 2 shows slightly shorter runouts compared to the third scenario with the simplified ellipsoid shape despite the blocks sharing the same moments of inertia. They are otherwise similar, but these shape simplification tests should be performed once the rotational stability of the model is corrected. With the most realistic paths and a wider lateral spreading of the three software programs, RAMMS’s results are promising for quantitative hazard zoning purposes on such challenging sites once carefully calibrated to reproduce proper runouts. The results of the second scenario, despite their underestimated runouts, show a wider covered area that matches better with the observations (Figure 5 and Figure 6). With the simulated trajectories also being exported in vector format, the proper addition of path width to pass from the propagation probabilities to the reach probabilities is facilitated. For low hazard probabilities, the number of simulated trajectories should however be increased for the hazard contours to stabilize. This is manifested by the irregular contours mostly visible for the 1/1000 and 1/5000 yearly hazard probabilities (Figure 6). From Equation (2) with the frequency used, an exposed object being reached by 1 trajectory out of 1625 would have a yearly hazard probability of 1/5000. So, 1000 simulated trajectories are clearly not sufficient for defining 1/5000 hazard zones quantitatively for the frequency used.

As expected, the process-based simulated velocities (Figure 4c) do not show any cutoff and show the typical rockfall sawtooth profiles. Interestingly, despite being underestimated by approximately half of the reconstructed velocities of block 2 from 900 to 1400 m, the simulated velocities slowly increase, stabilize, and decelerate, as observed for block 2. The simulated freefalling parabolas are in general more vaulted, as indirectly shown by the gentle decrease in velocity at the beginning of most sawtooth profiles followed by a steeper increase in velocity compared to the reconstructed velocities. This also indirectly indicates that the simulated bounce heights would have been overestimated if the velocities did correspond, which is safe. For the final fine-tuned simulations, one must remember in that case to also double the translational energies since only half of the volume of block 2 could be simulated.

Back analyses from field observations are especially important with the tested version of RAMMS because the parameters might have to be fine-tuned far from their anticipated values to produce safe realistic runouts. Many historical paths and/or deposited block fragments should then be mapped for the necessary meticulous time-consuming iterative fine-tuning of the simulation’s parameters. With the most realistic simulated paths of the three software programs tested, RAMMS’s results are promising for quantitative rockfall hazard zoning and related risk assessment tasks if simulated in sufficient numbers and if properly tuned from extensive back analyses. This is in line with the conclusions of Jarsve [67].

4. Benchmark Summary and Recommendations

The tested versions of the three software programs showed their strengths and weaknesses for reproducing such challenging events involving large 200 m³ block fragments, which is far beyond what the models were originally calibrated to. Future versions of the software programs may lead to different conclusions as models improve. Additionally, the comparisons were performed from simulation results taken at early calibration stages and could have been refined to a certain extent. On that topic, adjusting the runout distances is possible with each software program. Rockyfor3D stands out in this case by first quickly predicting safe runout distances close to the cumulative distribution of the observed blocks with its objective rapid automatic parameters. Any user would thus predictably reproduce the same results for this site with this objective set of parameters. Rockyfor3D also manage to reproduce relatively well the observed runout distances from parameters that objectively correspond with the terrain properties of the site. Based on this, it was shown to be a predictable model on the aspect of the runout distances. As for any models, one should reproduce such a comparison approach on other sites to evaluate the extent of its predictability. In certain conditions, its predictability may be impacted by unexpected erroneous ballistic freefall behaviors as those encountered in Appendix B.

The two other software programs produced results that followed the opposite direction than the previous safe approach from Rockyfor3D. For Flow-R, the runouts were always slightly shorter than what was anticipated based on the reach angles used or the first results using parameters for similar sites. This could be easily corrected using Equation (5) with subsequent scenarios thanks to the efficiency of the model, but one must not neglect the per-site calibration aspect from field observations. With RAMMS, the parameters had to be tuned beyond what seems to correspond to the terrain properties of the site to approach the desired runouts, and they should still be tuned further for subsequent scenarios.

However, adjusting the lateral deviation of Rockyfor3D’s simulations through its parameters is not trivial. Therefore, regardless of subsequent fine-tuned simulations, one should not expect the lateral deviations to differ much from those obtained. In this aspect, Flow-R offers more control on its geometrical hydrologic “rockfall” diffusion if one is willing to experiment beyond the intended exponent limits while carefully adjusting the susceptibility threshold. Of the three tested software programs, RAMMS::ROCKFALL, stands out here by better reproducing the lateral deviations objectively associated with the shapes of the blocks. The predictability issues related to its underestimated runout distances with moderately to very flat blocks are likely to be solved once its gyroscopic stability is fixed (Leine et al. [66]).

This comparison of simulation results taken at early calibration stages shows that there is room for improvements of the current software tested. Depending on the requirements of the rockfall hazard stage involved, combining the use of software programs with different behaviors could be a workaround in the meantime as follows:

• One could quickly evaluate the normality/abnormality of a site and preliminary reachable zones by predictably tracing desired reach angle contours from effective and nongrid-biased geometric approaches with, for examples, CONEFALL or QPROTO GIS tools (Jaboyedoff and Labiouse, [17]; Žabota and Kobal, [68]; Scavia et al. [69]; Castelli et al. [70]). Lateral constraints can be roughly estimated from common empirically observed lateral deviations in the literature (

  Table 2. Nonexhaustive compilation of planimetric lateral deviation observed for different sites.

  Observation Type	Number	Shape	Volume or Mass	Lateral Deviation	Source
  2015 Mel de la Niva rockfall event	2 blocks (largest)	Moderately to very flat and slightly elongated	~200 m3, ~500 t	±19°	Noël et al. [1]
  Rockfall event	2 rocks	Not flat, not elongated, ~ellipsoidal, d1 of 1.6 m and 2.5 m	~1.4 m3 for the first rock, ~3750 kg and ~5600 kg	±3.5°	Wyllie [41]
  Full-scale forested experiment	69 rocks	Natural and artificial rocks	69 rocks, 16–200 kg	±30°; 93% in ±15°	Ushiro and Tsutsui [71]
  Full-scale forested and non-forested experiment	218 rocks	Not flat, not elongated, ~spherical, mean d1 of 0.95 m	0.1 to 1.5 m3 280 to 4200 kg	±15°	Dorren and Berger [72]
  Full-scale experiment with reinforced concrete rocks	178 rocks	Not or moderately flat, not elongated	45 to 2670 kg (All rocks)	±19°; 94% in ±15°	Caviezel et al. [73]
  	22 rocks	Not flat, not elongated	45 kg	±16°
  	28 rocks	Not flat, not elongated	200 kg	±12°
  	43 rocks	Not flat, not elongated	800 kg	±14°
  	16 rocks	Not flat, not elongated	2670 kg	±7°
  	4 rocks	Moderately flat, not elongated	45 kg	±3°
  	23 rocks	Moderately flat, not elongated	200 kg	±15°
  	23 rocks	Moderately flat, not elongated	800 kg	±16°
  	19 rocks	Moderately flat, not elongated	2670 kg	±18°

  ).
• Then, the use of relatively predictable and time efficient models, such as Rockyfor3D, could be used to estimate the expected runout distances from parameters that can be measured in the field and/or with the conservative parameters with slope thresholds of the rapid automatic simulations.
• Finally, the realistic lateral extent and energies could be estimated from models, such as RAMMS::ROCKFALL, once its runouts are tuned to match those of the previous conservative/predictable simulations and the onsite observations. One must remember to adjust the energies if the mass of the simulated blocks is constrained by software limitations.

The comparison performed at the challenging site of Mel de la Niva highlighted significant limitations inherent to each model. This reminds us of the inconvenient importance of calibrating the simulations of the current software on a per-site basis from onsite observations, back analyses and back calculations. Even if one may be tempted to trust simulation results blindly to save resources, it is important to be aware that no model is perfect at reproducing reality. Thus, as described in Appendix C, a simple validating back calculation based on rockfall ballistics (e.g., Appendix B) and on field observations should complement simulations. On the simulation software side, providing a back calculation tool would remind and encourage their users to validate their simulation results from field observations and fine-tune the simulations if needed. This is particularly important if the simulations are used for the design of mitigation measures.

5. Conclusions

Simulations from Flow-R 2.0.9, Rockyfor3D 5.2.15 and RAMMS::ROCKFALL 1.6.70 were qualitatively compared to geometrically drawn reach angle contours based on the CONEFALL method and to the mapped rockfall trajectories and deposited block fragments of the Mel de la Niva site from Noël et al. [1]. A focus was placed on evaluating the objectivity, predictability, and sensitivity of the models by comparing simulation results taken at early calibration stages from back analysis for this challenging site. In this way, one can estimate how much the observed runouts may be missed by the simulations when using a reasonable set of parameters for similar sites and how good/poor the predictions may be if overlooking the important back analysis and calibration processes from onsite observations.

The comparison of three scenarios also allowed us to qualitatively evaluate how the models react to changes in their parameters and how far they need to be adjusted beyond what one would have anticipated. The comparison of simulation results taken at early calibration stages showed that there is room for improvements in the current software tested and that combining the use of software with different behaviors could be a workaround in the meantime.

The findings with the current evaluated software also highlight the significant limitation caused by the inconvenient required calibration of the simulations on a per-site basis from onsite observations, back calculations and back analyses. Calibration and validation reasoning steps were given and demonstrated with the presented comparison. Finally, a back calculation tool is provided as Supplementary Material to help and encourage the reader to repeat such comparison and back analyses on other sites.

From a longer-term perspective, this work emphasizes the need for resources to be spent to better understand rockfall dynamics, such as the novel rockfall experiments performed by the SLF and Geobrugg (Caviezel et al. [43,73,74]; Rigenbach et al. [75]; Noël et al. [44]). Hopefully, a better understanding of rockfall phenomena will help improve the objectivity and predictability of rebound models.