Accounting for the effect of forest and fragmentation in probabilistic rockfall hazard

10 The presence of trees along the slope and block fragmentation at impact strongly affect rockfall dynamics, and hazard as a consequence. However, these phenomena are rarely simulated explicitly in rockfall studies. We performed rockfall simulations by using the 3D rockfall simulator HY-STONEHy-Stone, modelling both the presence of trees and fragmentation through specific algorithms implemented in the code. By comparing these simulations with a more classical approach that attempt to account implicitly for such phenomena in the model parameters, and by using a new probabilistic rockfall hazard 15 analysisProbabilistic Rockfall Hazard Analysis (PRHA) method, we were able to quantify the impact of these phenomena on the design of countermeasures and on hazard. We demonstrate that hazard changes significantly when accounting explicitly for these phenomena, and that a classical implicit approach usually overestimates both the hazard level and the 95 th percentile of kinetic energy, leading to an oversizing of mitigation measures.


Introduction
Rockfalls are widespread in mountain ranges, coastal cliffs, volcanos, riverbanks, and slope cuts, and threat people, structures and infrastructures, and lifelines (Crosta et al., 2015).Although rockfalls generally have a limited size, they are extremely rapid processes that exhibit high kinetic energies, long runout and damaging capability (Corominas et al., 2017).Rockfall hazard and risk assessment (Corominas et al., 2005;Agliardi et al, 2009;Lari et al, 2014;Wang et al., 2014;De Biagi et al, 2017;Farvacque et al, 2019Farvacque et al, , 2021;;Hantz et al., 2021) and the design of defensive works (Volkwein et al., 2009) require numerical modelling of rockfalls to assess the dynamics of the blocks (i.e., velocity, kinetic energy and bouncing height) and the lateral and longitudinal spreading (Agliardi & Crosta 2003).In Italy, for example, the design of rockfall barriers is based on the use of the 95th percentiles of the blocks' height in flight and their kinetic energy, obtained from numerical models (UNI 11211; Volkwein et al., 2011).Since rockfall dynamics depends on block geometry, slope topography, surficial geology, vegetation, and some peculiar rockfall behaviours (e.g.dynamic fragmentation), the reliability of analyses and the efficiency of rockfall protections depend on the accuracy of modeling predictions and on the correct account for all relevant methodological issuesthese variables (Crosta et al., 2015).Both the characteristics of the slope (e.g.topography, material properties and presence of forests) and the type of rockfall (e.g.whether it is fragmental) must be taken into account during modelling, because they contribute to the overall extent of rockfall potential and hazard zonation in mountain areas (Frattini et al, 2012).Both these characteristics can modify the trajectories, the extent and the dynamics of the rockfall events, the frequency, and the probability of impact.
Forests provide important protection against rockfall in steep mountain terrain, defending structures and infrastructures (Berger et al. 2002;Dorren et al. 2004a;Perret et al. 2004).Thanks to this nature-based solution, maintenance and installation costs of technical protection measures, such as embankments or nets, are financially bearable or can even be avoided at many places due to the reduction of rockfall rebound heights and impact energies by previous impacts on trees (Grêt-Regamey et al., 2008;Häyhä et al., 2015;Getzner et al., 2017;Moos & Dorren, 2021).Although this protective effect is evident in hazard assessment processes because it supports decisions on risk prevention measures, it is often accounted only in implicit terms, by adopting a set of modified restitution coefficients (Pfeiffer and Bowen 1989;Azzoni et al, 1995).More rarelyLess frequently, the presence of trees is simulated explicitely by using numerical modelling approaches (Dorren et al. 2006;Stoffel et al. 2006;Berger and Dorren 2007;Bigot et al. 2009;Jancke et al. 2009;Rammer et al. 2010;Leine et al. 2014;Radtke et al. 2014;Kajdiž et al. 2015;Dupire et al. 2016;Moos et al. 2017;Toe et al. 2018).
When stiff and strong rock blocks hit a hard impact substratum or other blocks of comparable size like a talus deposit, they may fragment and explode (Crosta et al., 2015).The rockfall fragmentation process is defined as the separation of the initial rock mass into smaller pieces generally upon the first impact on the ground (Evans & Hungr, 1993), and the resultant fragments propagate downslope following independent trajectories and new dynamics (especially in terms of kinetic energy and height) compared to the source block.This definition covers both the disaggregation of the block fragments delimited by pre-existing fracturesdiscontinuities in the initial mass and the generation of new fragments due to the breakage of intact rock (Corominas et al. 2012;Ruiz-Carulla, 2018).Block fragmentation is generally at the origin of extreme behaviors, major damages and accidents, and can interact strongly with protection structures (Nocilla et al., 2009;Wang et al., 2010;Corominas et al., 2019).
Even if fragmentation during rockfall is recognized as fundamental in risk analysis (Corominas et al. 2012), a complete understanding of the process during rockfall has not been achieved so far, remaining a phenomenon largely neglected during numerical modelling.Only a few numerical codes allow modelling propagation that explicitly takes into account fragmentation (Crosta et al., 2003;Frattini et al, 2012;Matas et al., 2017;Ruiz-Carulla, 2018).When missing an explicit algorithm, the modeling of rockfalls with fragmentation can be done with two alternative approaches: either the model is calibrated to replicate the spreading of the event, including the most distal fragments, or the model is calibrated to replicate only the main deposit, neglecting the most distal blocks.The first approach leads to hazard overestimation, the second to hazard underestimation.
The aim of this paper is to quantify rockfall hazard when accounting for the presence of trees and fragmentation with an explicit simulation approach (i.e. using specific algorithmsmodel), and to evaluate the differences with a classical approach that does not simulate explicitly such phenomena.The simulator Hy-Stone (Crosta et al. 2003;Crosta et al. 2004), which allows to model both the presence of forest and fragmentation, and a new revised Probabilistic Rockfall Hazard Analysis (PRHA) are adopted to quantify the impact of these phenomena on the design of countermeasures and on hazard.

Hy-Stone (HS)
The simulationanalysis of rockfall propagation was performed by means of Hy-Stone, a 3D rockfall simulator that reproduces the block motion from the dynamics equations (Crosta et al. 2004;Frattini et al. 2012;Dattola et al., 2021) using a triangulated vector topography derived from Digital Terrain Models (DTMs).The model allows to simulate blocks with the shape of spheres, cylinders, ellipsoids, and discs.The stochastic nature of rockfall processes and parameters areis accommodated by slope morphology and roughness, and by the random sampling of most parameters from different probability density distributions (e.g.uniform, normal, exponential).The block trajectories are computed by splitting them in a succession of elementary motions: free fly, rolling, sliding and impacts/bouncing.When the impact process is concerned, Hy-Stone has many different models comprising the constant and not-constant restitution coefficients (Pfeiffer and Bowen 1989) and the evolution of the elasto-visco-plastic model initially proposedformulated by di Prisco and Vecchiotti (2006) and furthersubsequently extended toby introducing rotation and prismatic blocks (Dattola et al., 2021).Specific sub-model components explicitly account for the interactions between blocks and countermeasures or structures, between blocks and trees, and fragmentation (Frattini et al. 2012).

Tree-impact algorithmsub-model
The block-forest interaction is modeled through a stochastic tree-impact algorithm.Treessub-model.Tree height, trunk diameter, absorbable energy, and density (as number of trees inper 10 square meters square) are used by the algorithmas input to calculate at each cell a probability of impact, that depends on the tree density, tree size and block size, and, in case of impact of impact, a loss of block kinetic energy and a lateral deviation of the trajectories (Frattini et al., 2012) are consideredcomputed.
The block kinetic energy lost by impact on tree stems is greatest for central impacts, and decreases according to a Gaussian distribution away from the stem axis, while the angular deflection of the block on impact is assumed to vary according to the type of impact (central, lateral, scour) (Dorren et al 2004b).

Fragmentation algorithmsub-model
Hy-Stone can simulate the splitting up of a block in fragments moving independently from each other.The fragmentation occurs at impact when the kinetic energy of a block exceeds a thresholdlimit energy defined by the relationship elaborated by Yashima et al, (1987).)based on the Weibull distribution.The Yashima expression is: (1) where  , is the limit energy,  is the Poisson's coefficient,  is the Young's modulus,   is the radius of the parent block,  0 is a reference volume,  0 is the strength at the reference volume and   is the Weibull distribution parameter.Coefficients   and   are computed according to the following expression: and Coefficient   is the result of the following expression: Therefore, the threshold fragmentation energy depends on the geomechanical properties of the block and its volume (the larger the block, the lower the fragmentation energy).Once the fragmentation criterion is satisfied, a distribution of fragments is generated with a size determined according to a power-law distribution.: where () is the fragment size distribution,  is the fragment diameter,   is the maximum fragment diameter and  is a model parameter.The maximum fragment size is a fixed fraction of the parent block size.The number of fragments is computed according to athe mass conservation criterion(the total fragments mass must be about the mass of the parent block) and the above distribution, and the energy of each fragment is calculated assuming thatby means of the following expression: in which   is fragment mass,  , , is the translational kinetic energy available after impact is equally distributed among all fragments.Fragment projection velocity and direction are thenof the fragment and   is a model parameter. is computed according toby imposing the translational energy conservation criteria.Results show.Once the kinetic energy of fragment is known the inverse formula gives the fragment ejection velocity modulus.Fragment ejection velocity direction is computed stochastically within a cone whose aperture is a model parameter.Frattini et al. (2012) showed that block fragmentation has an effect on the runout extent and on the spatial distribution of velocities and heights of the flying rocks.The largest fragments, however, display a behaviorbehaviour that is more similar to that of the parent blocks.

Rockfall hazard assessment
To assess rockfall hazard, we propose a new revised PRHA (Probabilistic Rockfall Hazard Analysis), based on Lari et al. (2014), to build rockfall hazard curves starting from a set of block-volume scenario simulations.This methodology owes its idea on Cornell's (1968) probabilistic seismic hazard analysis (PSHA), which considers all possible earthquake scenarios to provide the exceedance probability of a certain level of ground motion at a site within a defined time frame.For each blockvolume scenario, s, the probability of exceeding a certain value of intensity (i, i.e. the reach of a specific value of kinetic energy), for each position along the slope (z) is: (1) where p(I)  (  ) is the probability density function of kinetic energy at the position z. for the scenario .Multiplying the exceedance probability by the annual frequency of occurrence (f),  ), we obtain the annual rate at which i is exceeded, (2) The annual frequency of occurrence (f)  ) of each scenario combines the onset frequency (fo)  ) and the transit frequency (ft) , ) at a certain position and for the specific scenario: The onset frequency (fo) , ) of blocks with a certain volume, V,  , can be expressed in terms of magnitude -frequency relationships (Hungr et al. 1999;Dussauge et al. 2003;Rosser et al., 2007).
where ()(  ) is the cumulative number of individual blocks with volume larger than V;  for the scenarios ; parameter  depends on both the area extent and the overall susceptibility of the cliff, whereas the power law exponent, b, mainly depends on lithology and geological structure (Hungr et al., 1999).To properly account for the frequency of individual blocks that propagate on the slope, it is necessary to combine the volume frequency relationship of rockfall events with the volume frequency relationship of blocks (Hantz et al., 2018;Hantz et al., 2020).The first relationship can be developed from surveyed historical events (e.g.Dussauge-Peisser et al. 2002;Chau et al., 2003;Guzzetti et al., 2003;Guthrie and Evans, 2004;Malamud et al., 2004) and provides annual frequencies of released rockfall volumes.However, these volumes should not be used for hazard analysis because single rockfall events disaggregate or fragment (Ruiz-Carulla et al., 2017) soon after the detachment and during propagation into a distribution of smaller individual blocks.On the other hand, the volume frequency relationship of blocks can be derived from the rock mass fracture network or directly from already stopped blocks, both in the talus (Ruiz-Carulla et al., 2017), along roads (Hungr et al. 1999), and caught by rockfall nets within a certain range of time (Matasci et al., 2015;Moos et al., 2018).However, these distributions usually lack the temporal frame that allows to correctly estimate the annual frequency.The combination of the two distributions can be achieved by calculating the total volume of the event (integrating the first distribution) and by calculating the  parameter of the second distribution by assuming the total volume to be equal to the first one (Hantz et al., 2018).
The transit relative frequency (ft) can be calculated for the rockfall simulation and corresponds to the ratio between the number of potential paths passing through a position (t)  ) and the total number of simulated paths from the rockfall trajectories, (ttot): For rockfall scenarios with different magnitude that occur in a certain position along the slope, the total annual rate at which  is exceeded,   (  >  ̅  ), derives from the sum of these scenarios, : By assuming a homogeneous, stationary Poisson process for the occurrence of the events (Crovelli, 2000), the probability of exceeding each intensity i in the next T years from this annual rate, Ppoiss,  , is: This represents the hazard curve at each position along the slope.
With respect to Lari et al. (2014), the revised PRHA method adopts a more flexible non-parametric approach for the kinetic energy probability distribution.Moreover, the new PRHA implements the approach proposed by Hantz et al (2016Hantz et al ( , 2019) ) for the calculation of the onset frequency (fo),  ), using the frequency-size distribution of the blocks observed along the talus to downscale the magnitude-frequency distribution of larger study areas.

Demonstration case studies
The application of potential rockfall scenarios was performed at the two representative sites that were recently affected by rockfall events in the Aosta Valley Region (Western Italian Alps) showing a significant role of forest and fragmentation at Saint Oyen and Roisan (Figure 1), respectively.During both the events, the rockfalls impacted roads and buildings, thus requiring a practical implementation of hazard assessment (for zonation) and the design of protection barriers (for mitigation).
Saint Oyen and Roisan are located in the Western Alps, within the Austroalpine-Pennidic collisional prism, consisting of overburden layers formed by continental crust and fragments of oceanic lithosphere, strongly reworked by the Alpine tectonometamorphic processes (Dal Piaz et al., 2016).
In the Saint Oyen case study (45°48'59.0"N7°12'21.0"E),about 17,500 m 3 of Ruitor micascists detached in March 2020, and reached a service road and the playing field in the lower part of the slope, passing through a mature fir forest.The presence of the forest significantly influenced the blocks distribution along the slope, increasing the lateral dispersion of trajectories and reducing their mobility.The case study is well documented by UAV flights conducted by the Regional Authority soon after the events, allowing for a detailed mapping of arrested blocks on the slope.We adopted this case study to investigate the role of forest, which has been fundamental for the rockfall dynamic, as observed in the field.Although minor fragmentation may have occurred during the event, we neglected it during the simulation to focus on tree-impact only.
Less than 10 km far, at Roisan (45°47'49.3"N7°18'49.0"E),about 1,050 m 3 of Arolla gneiss toppled in October 2019 and impacted after 20 m of free fall (Polino et al., 2015) against a bench.While the main body of the rockfall stopped in a relatively flat area close to the source area, two blocks reached the foot of the slope causing the interruption of a municipal road.The event is documented by a post-event UAV flight, and by a detailed field survey of the blocks.For this case study, the presence of forest was minor due to the size and age of the trees and it has been neglected in order to reveal better the role of fragmentation.

Analysis and results
For both case studies, we firstly back-calibrated the model parameters on the rockfall events in order to simulate several volume scenarios from local-scale rockfall source areas (with and without the use of specific algorithms for tree-impact and fragmentation) to quantify the differences in terms of dynamics, spreading and rockfall hazard.We simulated both all the scenarios by using spherical blocks and a 3D topography derived from the available 1 x 1 m Lidar DTM of Aosta Valley Region.The characteristics of each simulation, the number of simulated blocks, and the parameters adopted when using the two algorithms are reported in supplementary Table S1, 2 and 3.

Calibration by back-analysis
The calibration of model parameters was obtained by fitting the longitudinal and lateral extent of rockfall eventstrajectories and deposits by using the Hy-Stone model with and without tree-impact and fragmentation.In particular, we simulated the following scenarios (Figure 1): -SO_HS (Saint-Oyen tree impact implicit): the values of parameters are modifiedset to account for the forest, e.g.
increasing rolling friction and reducing the tangential restitution coefficient.This is the most classical approach adopted in the practice to "simulate" the effect of forest with an implicit approach.
-SO_HStree (Saint-Oyen tree impact explicit): the values of parameters are calibrated by adopting the Hy-Stone treeimpact algorithm that explicitly simulates the effect of forest; in this case, the motion parameters used in the simulation do not account for the forest.
-R_HS (Roisan Fragmentation Implicit): the values of parameters are modifiedset to allow the model to replicate the spreading of the event, including the most distal blocks, implicitly accounting for the possibility of fragmentation.
-R_HSfrag (Roisan Fragmentation Explicit): the values of parameters are calibrated by adopting the Hy-Stone fragmentation algorithm that explicitly simulates the distal blocks as fragments.
For Roisan, we experimented a different calibration strategy that replicates the spreading of the main deposit only (R_HSshort), neglecting most distal blocks (Figure 1E).Although this strategy is physically correct to simulate non-fragmenting blocks, it provides an overall spreading that strongly underestimate the possible reach distance of fragments and the hazard level, accordingly.
For Saint-Oyen, both the simulations (RSO_HS, RSO_HStree) provide a good match with the main deposit of the 2020 event (Figure 1 A and B), with a slightly larger spreading when using the tree-impact algorithm, consistently with the fact that the impact with trees adds a component of lateral dispersion to the trajectories.
For Roisan, we can observe a good match between the longitudinal and lateral extent of the main deposit from the 2019 event and the simulations (R_HS, R_HSfrag) but we observe an overestimation of the blocks reaching the paved road (18 blocks modelled, while just 2 blocks during the event) when the fragmentation algorithm is not used (Figure 1A1C).The comparison with simulated stopping points shows that the model without fragmentation is able to reach the maximum distance, but not in the right location, since trajectories are strongly controlled by topography.This does not happen with the fragmentation algorithm, which is able to replicate the right position of the distal blocks in the meadow (Figure 1B).1D).
In addition, the volume distribution is also considered.Table 1 Table 2and Table 2 report the values of normal and tangential restitution coefficients and of the friction coefficient for the different slope materials used in the rockfall numerical simulations, in the cases of SO_HS and supplementaryR_HS.
Supplementary Table S1 and Table S2 report the parameters of tree-impact and fragmentation algorithms.used in the cases of SO_HStree and R_HSfrag.

Effect of tree-impact and fragmentation algorithmsub-models on kinetic energy
To quantify the effect of explicitly simulating tree-impact and fragmentation in rockfall modelling, we performed simulations for five volume scenarios in which the released volumes are changed (Table 3), using the modelingmodelling parameters that were back-calibrated from the events as previously described.The volume scenarios range from 0.001  3 to 100  3 to encompass the block sizes surveyed on the field at the two sites.
For the spatial analysis, we divided the slope into a 10 x 10 meters square lattice and we calculated statistics of kinetic energy within each square.

Effect of tree-impact algorithmsub-model
Figure 2Figure 2 shows the 95 th percentile of the blocks kinetic energy in each 10 m square.considering the first and fifth scenario with and without the forest sub-model.This statistic variable has been chosen since it is frequently used for designing defensive works (UNI 11211;Maciotta et al., 2015;Lambert et al., 2021).In the case of small volume blocks, the simulation without tree-impact algorithmsub-model (Figure 2a2A) shows a central sector characterized by the highest kinetic energies (from 2,500 kJ up to over 10,000 kJ for the 95 th percentile), and a distal zone characterized by lower values.Trajectories are able to reach the base of the slope, the unpaved road, buildings, and playing field, and overpass the location of the outermost blocks of the 2020 event.When using the tree-impact algorithmsub-model (Figure 2B) the number of trajectories passing through the central sector of the slope decreases dramatically.The trajectories that reach the base of the slope are concentrated in the area affected by the 2020 event where the forest is damaged.These trajectories reach only the unpaved road, with associated 95 th percentile kinetic energy values of less than 2,500 J.For large blocks (Figure 2 C and Dfifth scenario), the kinetic energy is high enough to nullify the effect of forest, and the two scenariosanalysis without and with tree-impact algorithmsub-models become similar.(Figure 2 C and D).
From these results, it is evident that the use of the tree impact algorithmsub-model is relevant in the case of small volume blocks, for which the simulated trees are able to interrupt most of the computed trajectories, and in any case to decrease the kinetic energies.On the contrary, tree-impact simulationanalysis is almost irrelevant for large-volume blocks.
AnalyzingAnalysing the distribution of kinetic energies along the road and the blocks number at the foot of the slope without (HS) and with (HStree) the tree impact algorithmsub-model, we systematically observe lower values of energy for the HS modelsanalyses (Figure 3).InIndeed, in these models, the effect of the forest is simulated by modifyingreducing the restitution and friction coefficients and increasing friction coefficient, calibrated on the range of kinetic energies of the event.However, this coefficients modification is independent of the kinetic energiesblock mass of the simulated blocks, and, therefore, it is not possible to observe the scale effect revealed in the HStree analyses.
When the kinetic energies are lower both than the calibrated kinetic energies and the kinetic absorption energies of the trees (scenarios S1 and S2), the classical HS approach overestimates the runout (see the large number of blocks intersecting the road after crossing the forest in Figure 2 A).Instead, the HStree algorithm intercepts, slows, and stops the least energetic blocks, allowing only the most energetic to reach the lower part of the slope.As a result, few transits are obtained, but with much higher kinetic energies due to the filtering effect of the forest (Figure 3).
WhenIn contrast, when the kinetic energies grow beyond the calibration range (scenarios S4 and S5), the classical HS approach continues to apply the forest effect (through the modified parameters) even though the kinetic energies are well above the tree absorption energies, underestimating the runout (the number of blocks intercepting the road remains about the same as in the low energy scenarios) and the kinetic energies (Figure 3).In contrast,For the HStree algorithmanalyses in these scenarios S4 through S5, showsshow higher kinetic energies and a high number of transits (compared to the lower-volume scenarios) because the effect of the trees becomes negligible, as it should be (Figure 3).
In the intermediate S3 scenario, greater congruence between the two approaches HS and HStree is observed (Figure 3) because the simulated volumes are more coincident withsimilar to the calibration range (between 0.001 m 3 and 34 m 3 ).

Effect of fragmentation algorithmsub-model
Figure 4 shows the 95 th percentile of the blocks kinetic energy in each 10 m square grid with and without the adoption of the fragmentation algorithm.sub-modelconsidering the first and fifth scenarios.The behaviorbehaviour of the modelsanalyses with small or large block volumes is extremely different.In the case of small volume blocks, the adoption of fragmentation algorithm is almost negligible, because blocks are too small to undergo fragmentation.In Figure 4 A the highest 95 th percentile values of kinetic energy for the first scenario and without fragmentation are concentrated in the area located just below the modelled source and at the highest escarpment, and only four trajectories characterized by values up to 3 kJ reach and cross the paved road.In Figure 4 B we observe that the highest 95 th percentile values considering the fragmentation are concentrated in the area close to the wallcliff, but only one trajectory passes the road, characterized by 95 th percentile of kinetic energy much lower (up to 1.5 kJ).
For larger blocks (S5 scenario), the difference with and without fragmentation is much more significant because more than half of the blocks are fragmented (612 out of 2646, 23%).In Figure 4C without fragmentation the runout achieved by blocks does not exceed that of Figure 4a4A, but with much larger values associated with the 95 th percentile of kinetic energy reached all over the slope.The area located just below the modelled source and in the highest escarpment is characterized by kinetic energy values greater than 50,000 kJ at the intersection with the unpaved road.Values remain high also at the intersection with the paved road.In Figure 4D in which the fragmentation is considered there is an increase in the number of blocks crossing the roads, a consequent spread of trajectories with longer runouts (more than those actually achieved during the event) and a decrease in kinetic energy due to block fragmentation.On the unpaved road, the values associated with the 95 th percentile drop to 50,000 kJ, and where the event boulder stopped it decreases to 8,000 kJ.At the intersection with the paved road, percentile values are more frequently lower than 15,000 kJ except in an isolated section where they reach 50,000 kJ and over.
AnalyzingAnalysing the distribution of Kinetic Energieskinetic energies along the paved road at the foot of the slope without (HS) and with (HSfrag) the fragmentation algorithmsub-model, we systematically observe higher values of energy for the HS modelsanalyses (Figure 5).Although during the event very few blocks crossed the paved road and only two of them reached the meadow at the foot of the slope, the calibration of the model without fragmentation was accomplished by adjusting the parameters in order to reach the maximum runout.This causes a strong overestimation of the number of blocks crossing the paved road, a general overestimation of the landslide runout, and therefore also an overestimation of the kinetic energies at the element at risk.As already said for the Saint-Oyen case study, the runout in the HS models is almost independent from the kinetic energy of the blocks.blockmass.Therefore, the number of transits is roughly constant in all sixfive scenarios.
InInstead, in the HSfrag approach, the kinetic energy at the element at risk is systematically lower because the model is calibrated in order to allow only trajectedejected fragments (characterized by much lower volumes with respect to original blocks) to reach and cross the paved road as occurred during the event.The number of fragments reaching the road increases significantly through different volumes scenarios (from S1 to S5).This depends on the relationship between block size and fracture energy (Yashima et al., 1987); according to this relationship, the fracture energy scales with the radius of the block by an exponent that depends on the Weibull's coefficient of uniformity, and is always lower than 3, which is the scaling of the kinetic energy with radius.Hence, the larger the block, the higher is the probability of fracturing for a certain velocity.
The two approaches HS and HSfrag provide similar results in the S4 scenario (characterized by simulated volumes that are more coincident withsimilar to the calibration range, between 0.5 m 3 and 23 m 3 ) both in terms of number of blocks intersecting the road and in terms of kinetic energies: compared to all other scenarios, less than an order of magnitude separates the two 95 th percentile values of kinetic energy.

Rockfall Hazard
TheAs explained in the metholody section, the assessment of rockfall hazard requires the onset frequencies  0  0, for each magnitude scenario, the transit frequency, ft , and the distribution of kinetic energy in each position along the slope.
For the calibration of onset frequency parameters (eq.5),10), we adopted the methodology of Hantz et al (2018) to obtain the frequency distribution of blocks for different volume classes, by combiningwho related the magnitude -frequency relationship of all rockfall events within a fixed site with the size -frequency relationship of blocks along the talus.The first was for a specific event in the same site.
We obtained magnitude-frequency relationship by analyzinganalysing the available rockfall database of the Valle d'Aosta region, which includes 306 events with volume information (Figure 6).Among them, only 25 belongs to the same catchment of the case studies (Buthier catchment, Figure 6).ThisSince this subsample appears to be insufficient to characterize the magnitude -frequency curve, especially for smaller volumes, and that are not recorded, we therefore adopted the entire inventory that we fitted with a maximum likelihood approach, obtaining a good power-law fitting ( 2 = 0.99) for rockfalls larger than 10 m 3 , with a scaling exponent of 0.56.We believe that this parameter value is reliable also for the Buthier catchment, since the fitting curve has the same slope of larger rockfalls volumes (with a volume greater than 500 m 3 ) within the subsample.TheTherefore, this parameter  is adopted for the two cases studies.
For the parameter , we used the size -frequency relationship of blocks along the talus was obtained by image analysis only for Saint-Oyen due to a larger number of blocks and a better quality of the imagery.Figure 7 shows an excellent power-law fitting ( 2 = 0.99) for blocks larger than 0.2 m 3 , with a scaling exponent  equals to 1.22.Eventually, by combiningrelating the two distributionsmagnitude-frequency size-frequency and accounting for the potential unstable area of both case studies, we obtained an  value of 0.0072 and of 0.0021 for Saint-Oyen and Roisan, respectively.The resulting onset frequencies for the different volume scenarios are reported in Table 3 for both case studies.
Both transit frequency (ft) , ) and the distribution of the kinetic energy come from the rockfall simulation trajectories sampled within 10x10 m cells.In order to characterize the kinetic energy distribution, we tested the hypothesis adopted by Lari et al, (2014), who assumed the logarithm of the kinetic energy to be normally distributed, and obtained the kinetic energy probability density p(I)  (  ) by using the mean and standard deviation statistics.The Kolmogorov-Smirnov's test (Figure 8) shows that the normality is not rejected for more than 50% of the 10x10 m cells when using Hy-Stone without additional algorithms.
However, this percentage is lower when using the tree-impact and fragmentation algorithms, suggesting that a non-parametric approach should be adopted when the level of complexity increases.
By combining the various scenarios and taking into account their associated probabilities, (Equation 12), we constructed the hazard curves (by equation 8),13), which show the probability of exceeding a certain level of intensity in 50 years.
Figure 9 Example of hazard curves characterized by a non-logarithmic trend, calculated in five cells of R_HS model.
Figure 9 shows hazard curves only for some representative cells.We can assert that hazard curvesthey do not always have a logarithmic distribution, and that some curves (not here reported) do not reach the exceedance probability of 0.1 due to a very low transit frequency.ThenSubsequently, for each location along the slope and for each model runanalyses (SO_HS, SO_HStree, R_HS, and R_HSfrag), we reducedcomputed from each hazard curve to a single value in order to represent, fixing the hazard through a hazard map.Asexceedance probability in 50 years at 10% as done by Lari et al (2014), we chose the corresponding kinetic energy with a 10% chance of being exceeded in 50 years, and we extracted this value from each hazard curve.which is used to represent the hazard through a hazard map (Figure 10).
Compared to the SO_HS model, (Figure 10A), in SO_HStree (Figure 10B) the hazard decreases because kinetic energy is significantly lowered, except in correspondence of the two sectors most affected by the event (see calibration in Figure 1)Figure 1) where it remains similar (Figure 10(Figure 10 A and B).The total area involved remains about the same, although with slightly lower runout.However, if only the areas with Ek>1 kJ > 1 are considered, hazard decreases significantly along the road at the foot of the slope.
For the Roisan case study, compared to the R_HS model, in R_HSfrag the hazard decreases because the kinetic energy is significantly lowered, but note that the area involved increases (Figure 10 (Figure 10C and D).Analysing the distribution of the hazard values (Figure 11)(Figure 11) at the foot of the slope obtained by the different approaches without and with the tree impact and fragmentation algorithms, we observe an overestimation of the potential hazard in both case studies.In the Roisan case study, the overestimation is particularly high because the chance to fragment the blocks into smaller fragments greatly reduces the kinetic energy of those.Moreover, the distribution is less sparse because the only blocks with an energy value higher than the minimum energy value (1 kJ) that are able to reach the foot of the slope are few and localized in a 10-meter corridor.

Discussion
When hazard and risk need to be assessed, it is required to have a repeatable procedure and possibly a unique result.This study demonstrates that different modelling approach can influence both the final result of hazard analysis and risk mitigationthe design of countermeasures, but also points out the problems involved in advanced modelling, leading to necessary discussions on the topic.

Tree impact
The classical approach for modelling rockfall propagating along forested slopes is based on the modification of restitution and friction coefficients, calibrated on the extent of block propagation.This study shows that the adoption of this set of modified restitution coefficients provides a correct replication of the maximum lateral spreading and longitudinal runout, but inaccurate dynamicenergy of blocks.In fact, the modification of the restitution coefficients is independent on the size of the blocks and can slow down even those blocks that are large enough to be actually unaffected by the presence of the forest.This leads to an overestimation of rockfall runout and of the number of blocks reaching the elements at risk.When the protective role played by the forest is explicitly simulated (HStree), the hazard decreases due to the forest protection, but the high-percentiles of kinetic energy become higher.This occurs because the trees interceptstop the blocks with lower kinetic energy, withgenerating a filtering effect of the larger blocks, leading to the risk of considering, paradoxically, the presence of the forest as more dangerous.These considerations open an important discussion on the opportunity to design the defensive works only based on percentiles of the kinetic energy.

Fragmentation
In case of rockfalls charcaterized by dynamic fragmentation, the classical approach for calibrating the model with this events is based on a conservative adjustment of the parameters in order to reach the maximum runout of single fragments.We demonstrate that this approach leads to a strong overestimation of the number of transits (Figure 4), the overall landslide runout, the kinetic energy of blocks impacting the elements at risk (Figure 5),(Figure 5), and the hazard (Figure 10).(Figure10).On the other side, the alternative approach to replicate only the main deposit, neglecting the most distal blocks would result in an underestimation of all these quantities (supplementary Figure S5).Therefore, regardless of whether deciding to simulate only the blocks that have stopped in the main deposit (Figure 1E(Figure 1E) or to extend the trajectories to the maximum fragment extent (Figure 1C(Figure 1C), this study demonstrates that the result is fundamentally incorrect, especially for the design of defensive works.On the other side, the explicit modelling of fragmentation is still challenging from both a theoretical point of view (Ruiz-Carulla et al., 2017;Shen et al., 2017;Guccione et al., 2022) and a practical point of view, due to the difficulty to calibrate the geotechnical parameters that control fragmentation.This adds further uncertainty in the analysis of rockfall dynamics and hazard.
From theThe results of hazard computation with and without when using the fragmentation algorithm, which show a decrease of the hazard decreases.However, as in case of forest, this results from the fact that more weightimportance is given to the kinetic energy of the blocks and not the frequency.Is it correct to infer that the hazard decreases, mostly due to the decrease of kinetic energy, even if the frequency increases, the kinetic energy of the blocks decreases, but and the trajectories are much more dispersed much more?We believe ?Our belief is that itthis inference is not completely correct because at spatial level, more zones are involvedentirely accurate since the area affected by rockfalls is larger, even if fromblocks are smaller blocks.
This discussion leaves room for future new studies.

Probabilistic rockfall hazard
The PRHA approach allows to quantify rockfall hazard in terms of hazard curves, thus describing the probability of exceeding a certain level of hazard.For each magnitude scenario, the approach overcomes the need of selecting a statistic of the kinetic energy at a certain positon along the slope (Agliardi et al, 2009;Farvacque et al, 2021), and allows to consider the full energy distribution within a certain grid cell.With respect to Lari et al. (2014), the revised PRHA method presents two improvements: (i) it adoptsthe adoption a more flexible non-parametric approach for the kinetic energy probability distribution, instead of assuming a log-normal distribution, that we demonstrate in this paper to be frequently violated if tree impacts and fragmentation existsubsist (Figure 8); (ii) it implementsthe implementation of the approach proposed by Hantz et al (2016,     (1) Only for the explicit approach (HStree) To this purpose the In-situ Block Size Distribution (IBSD) at the cliff considered in the previous numerical simulations was obtained previously by the Geological survey of Valle d'Aosta by means of a terrestrial laser scanner survey.In the Figure S1 a comparison of the Rock Block Size Distribution (RBSD) obtained with the ortophotos at the toe of the slope, and the distributions obtained with the scenario R_HS, R_HSfrag is shown.The comparison reveals a good agreement since the curves are parallel each other although the Hy-Stone distributions overestimate the in situ one.

Figure 1 AFigure 2 Figure 3 Figure 4 Figure 5 Figure 6
Figure 1 A) Location of the two case studies in Aosta Valley region.The other panels show the back calibration of the rockfall events: A) and B) simulation of Saint-Oyen rockfall (A) with parameters modified to account for the forest (SO_HS) and (B) 615 by adopting the Hy-Stone tree-impact algorithm (SO_HStree); C) and D) simulation of Roisan rockfall (C) with parameters

Figure 7 Figure 8
Figure 7 The size frequency relationship of blocks along the talus obtained by image analysis for the Saint-Oyen event.

Figure 9
Figure 9 Example of hazard curves characterized by a non-logarithmic trend, calculated in five cells of R_HS model.655 ) and tangential restitution (  ) coefficients and of the friction coefficient (  ) for the different slope materials used in the rockfall numerical simulations for the Saint-Oyen case study.

Table 1
Values of normal (

Table 2
Values of normal (  ) and tangential restitution (  ) coefficients and of the friction coefficient (  ) for the different slope materials used in the rockfall numerical simulations for the Roisan case study.