Evaluation of a 3-D rockfall module within a forest patch model

Many slopes in the Alps are prone to rockfall and forests play a vital role in protecting objects such as (rail) roads and infrastructure against rockfall. Decision support tools are required to assess rockfall processes and to quantify the rockfall protection effect of forest stands. This paper presents results of an iterative sequence of tests and improvements of a coupled rockfall and forest dynamics model with focus on the rockfall module. As evaluation data a real-size rockfall experiment in the French Alps and two 2-D rockfall trajectories from Austria and Switzerland were used. Modification of the rebound algorithm and the inclusion of an algorithm accounting for the sudden halt of falling rocks due to surface roughness greatly improved the correspondence between simulated and observed key rockfall variables like run-out distances, rebound heights and jump lengths for the real-size rockfall experiment. Moreover, the observed jump lengths and run-out distances of the 2-D trajectories were well within the stochastic range of variation yielded by the simulations. Based on evaluation results it is concluded that the rockfall model can be employed to assess the protective effect of forest vegetation. Correspondence to: W. Rammer (werner.rammer@boku.ac.at)


Introduction
Forests in mountainous regions provide important protection functions for society and particularly the protection against rockfall has attracted considerable attention recently.Reasons are that rockfall is one of the most common natural hazards in mountainous landscapes, and that the protective effect of many European mountain forests may decrease in the future due to abundant old-growth phases with a lack of regeneration (e.g., Ott et al., 1997;BFW, 2004).To study the development of the protective effect as well as the influence of forest management, reliable simulation tools are required which are able to take into account the spatial pattern of rockfall processes on slopes as well as the effect of forest vegetation composition and structure on run-out distances and kinetic energies of falling rocks (Stoffel et al., 2006;Dorren et al., 2005).Due to the long time horizons of regeneration processes in mountain forests, models must be able to cope with time scales of several decades.Of particular importance are trade-off relationships between rockfall protection and other forest services and functions within the framework of multifunctional forest management.
Rockfall is defined here as a fast gravitational movement of rock boulders, rolling, tumbling or sliding down a hillslope (Selby, 1995).Boulders act rather independently and the boulder size considered for trajectory analysis is typically below 5 m 3 .
A wide range of rockfall models exists, covering different spatial scales and levels of complexity.Historically, among the earliest approaches are 2-D rockfall models that simulate rock trajectories along a slope profile (e.g.Bozzolo and Published by Copernicus Publications on behalf of the European Geosciences Union. W. Rammer et al.: Evaluation of a 3-D rockfall model Pamini, 1986;Azzoni et al., 1995;Pfeiffer and Bowenm, 1989;Pfeiffer et al., 1993).Typically, such models have been applied to support the design of technical rockfall protection measures.With the evolution of digital mapping and spatial modeling techniques a new class of 3-D models, simulating individual trajectories based on a digital elevation model, became increasingly popular (e.g.Agliardi and Crosta, 2003;Dorren et al., 2004;Lan et al., 2006).Such models are typically used on local to regional scales.
However, although forests can offer an effective protection against rockfall (e.g.Bigot et al., 2009) only a small number of rockfall models take into account the mitigating effect of trees on forested slopes.These models consider the effect of forest vegetation either implicitely (e.g.Brauner et al., 2005;Wehrli et al., 2006;Berger and Dorren, 2007), using aggregated stand variables like stand density or average diameter at breast height, or spatially explicit (Dorren et al., 2006), simulating individual tree hits along a rocks trajectory.Wehrli et al. (2006) conclude that a spatially explicit representation of forest structure would enhance the simulation approach.
Recently Woltjer et al. (2008) developed a spatially explicit 3-D rockfall model and coupled it with the process oriented 3-D forest patch model PICUS v1.4 (Lexer and Hönninger, 2001;Seidl et al., 2005).In Woltjer et al. (2008) the new rockfall model was tested for parameter sensitivity and for its ability to assess rockfall protection effects in line with static formulations from the literature (Gsteiger, 1993;Brauner et al., 2005).Moreover, the model showed plausible sensitivity to forest management as implemented in the forest model.While the forest model has been evaluated rigorously in several previous studies (e.g., Seidl et al., 2005Seidl et al., , 2009)), an evaluation of the rockfall module against empirical data is presented in this article.
The objectives of this study are twofold.First, recent model improvements are introduced, and second, the model is evaluated against empirical rockfall data from four test cases from France, Switzerland and Austria.In an iterative manner results from comparisons of model simulations and evaluation data feed back into model enhancement.Finally, conclusions on the applicability of the model are drawn.

The forest dynamics model
The forest model used in PICUS Rock'n'Roll is the hybrid forest patch model PICUS v1.4, which incorporates elements of a classical patch model as well as a stand level production model based on physiological principles.A detailed description of the hybridization is provided in Seidl et al. (2005).
Here, just a brief overview of the core model concept is given.
Basic modelling entity is the individual tree above a threshold of 1cm diameter at breast height.Tree biomass is arranged horizontally on 10×10 m 2 patches and vertically in 5 m canopy layers.A three-dimensional light model, allowing for the explicit consideration of direct and diffuse radiation within the canopy, is used to model inter-and intra tree species competition.The incorporation of a module that estimates NPP (net primary production) based on the concept of radiation-use efficiency (Landsberg and Waring, 1997) enhanced the robustness in predicting growth along environmental gradients.It also improved the physiological foundation of the model with regard to changing environmental conditions.The model requires input of incoming radiation, temperature, precipitation and vapor pressure deficit in at least monthly resolution.

The rockfall model
This contribution focuses on the evaluation of the rockfall model Rock'n'Roll as presented in detail in Woltjer et al. (2008).This section presents a brief overview of the model concept.The rockfall model simulates spatially explicit trajectories of spherical rocks of variable size (diameter, mass) on a three-dimensional slope.A rock trajectory consists of a series of rolling and/or jumping motions.The model uses a lumped mass approach for simulating free fall and ground impacts, and a rigid body approach for the simulation of rolling motions and tree impacts.Kinetic energy of a falling boulder is dissipated by (i) rolling friction for rolling motions, (ii) ground impacts for jumping motions, and (iii) tree impacts during rolling and jumping motions.
In a simulation rock trajectories start with a predefined initial velocity either from predefined or from randomly selected positions within the simulated area.A rock stops, when its kinetic energy falls to zero.
The slope surface is characterized by coefficients of restitution for jumping mode (i.e.rebound parameters) and by rolling friction coefficients for rolling mode.For the current study the rolling friction coefficient was kept constant at a value of 0.5 which is considered as a sensible default value for various ground types (see Azzoni et al., 1995).Coefficients of restitution describe the ratio of the normal (or tangential) velocity component of the center of mass before and after an impact (e.g.Kharaz et al., 2001).Parameter values for various surface types are retrieved from the literature (e.g., Azzoni et al., 1995;Azzoni and de Freitas, 1995;Chau et al., 2002;Schweigl et al., 2003;Heidenreich, 2004).
The rockfall model includes several stochastic components: (i) the lateral deviation and the jump angle after ground impacts, and (ii) the lateral deviation after tree impacts.Moreover, the coefficients of restitution are subject to uncertainty and may vary within a specified range of variability.Consequently, simulated trajectories may show substantial variability with regard to key variables like jump lengths or rebound heights.
In order to utilize frequently available 2-D trajectories for model validation and calibration a special mode was implemented in the rockfall model.Documented 2-D trajectories from single rockfall events usually report jump lengths and slope angles along the actual trajectory.In a simulation the rockfall model uses the trajectory profile as slope topography and forces the falling boulder along the profile.Within that "2-D mode" all modeled rockfall processes except tree hits are applicable and energy conservation laws are satisfied.
Rock'n'Roll is implemented in the C++ programming language.The rockfall model is integrated with the PICUS forest model and can be invoked every simulated year to assess the effect of vegetation on rockfall activity.The software is optimized for high calculation performance.For instance, the algorithms of rock movement are not based on fixed time steps but rather rely on analytical equations of motion to foster efficient computations.
Based on findings of preliminary model evaluation exercises and conceptual limitations of the modeling approach some components have been revised and extended.In the following sections these model enhancements are introduced.

Representation of surface topography
A major enhancement compared to the model variant in Woltjer et al. (2008) is that in the current version surface topography is modeled as a triangulated irregular network (TIN).A TIN describes a surface as a set of non-overlapping triangles with variable density.Besides raster based digital elevation models, TINs are frequently used for the representation of three dimensional surfaces, and modern GIS software packages usually support the conversion from raster based digital elevation models (DEM) to TINs and vice versa.While the vast majority of all available rockfall models use raster based DEM, we implemented a TIN-based approach.Reasons for the decision to use TINs are, inter alia, that the variable density of triangles allows a memory and time efficient computation of slope regions with more uniform terrain features, and that the flexible spatial resolution of a TIN is easier to couple with the fixed 10×10 m patch size of the forest model.
The spatial distribution of different surface properties, e.g.coefficients of restitution, is incorporated in the model by using a raster based approach.Surface properties are mapped in a GIS environment using arbitrarily shaped polygons which are imported as grids into the coupled rockfall-forest model.

Tree hits
In the original model variant the algorithm for the simulation of tree hits was designed as a set of mechanical and geometrical formulations based on fracture energy experiments in the laboratory.However, in reality a rock impacting a tree triggers far more energy-consuming processes, like deformation of the root-soil system or the swaying of the whole tree.Recently published empirical data on that issue (Dorren and Berger, 2006) and results from a numerical tree impact model (Jonsson, 2007) allow the implicit inclusion of these processes.
Our implementation of the tree hit algorithm is very close to the published model of Dorren et al. (2005; see also Dorren and Berger, 2006).The maximum amount of energy that can be dissipated by a tree during an impact (E diss,max ) depends on tree species and is related to the diameter of the tree at breast height via an exponential relationship (Eq.1).
For impact positions off the central bole axis only a fraction of the dissipation potential is exploited.This is expressed by a sigmoid function relating the consumable fraction of the dissipation potential to the impact position.If the kinetic energy of a rock during a frontal impact is higher than the dissipation potential, the tree breaks.
The horizontal deviation of the rock after the tree hit is modeled by deviation matrices, which define deviation ranges for three different impact types (frontal, lateral and scratch).Within this range a random number determines the realized deviation angle.For situations where the energy of the falling rock is high compared to the energy dissipation potential of the tree an additional upper deviation limit avoids unrealistic deviation angles.This limit is estimated as the hypothetical maximum deviation of the rocks' trajectory given the energy dissipation potential of the focal tree.

Real-size rockfall experiments Vaujany
From 2001 to 2003 real-size rockfall experiments have been realized on a forested and an adjoining unforested test slope in the Forêt Communale de Vaujany in France and have been described in detail in several publications (Dorren et al., 2005(Dorren et al., , 2006;;Bourrier et al., 2009).Here we briefly introduce the experimental setup.The test sites are located on a post-glacially developed talus cone, which mainly consists of rock avalanche, snow avalanche and rockfall deposits and is situated on altitudes ranging from 1200 m to 1400 m (a.s.l.) with a mean slope gradient of 38 • .A digital elevation model with a resolution of 2.5×2.5 m was available for both sites and surface characteristics (diameters of rocks covering the slope, etc.) were grouped into surface types and mapped throughout the test slopes (see Table 1).
The unforested Site 1 (see Fig. 1 for a profile) has the morphology of a channel for the first 240 m (projected) between the starting point and a crossing forest road.From there Section B extends for 100 m (projected) to another forest road, below which the slope ends in the valley bottom (Section C).The first section after release is mainly covered by debris with diameters of 8-10 cm (type (a) in Table 1), Section B is dominated by bigger blocks between 40 and 80 cm (type (b) in Table 1) while Section C is covered by fine debris and soil material (type (c) in Table 1).
The forested Site 2 represents a typical rockfall slope in the European Alps (see Fig. 1).The main tree species with regard to basal area are Silver fir (Abies alba, 72%), beech (Fagus sylvatica -14%), Norway spruce (Picea abies (L.) Karst.-6%), and Sycamore (Acer pseudoplatanus L. -6%).The stand is characterized by a mean tree diameter at breast height of 31 cm and a stand density of 290 trees per hectare.The coordinates of each tree on Site 2 were mapped.In the upper half of the slope at Site 2 the spatial variability (i.e.patchiness) of the main surface types (e) and (f) (compare Table 1) was considerably higher compared to Site 1.
A total of 218 rocks with a mean diameter of 0.95 m (mean volume=0.49m 3 , S.D.=0.3 m 3 ) were released individually by an excavator from a forest road.Three-dimensional rock trajectories were captured by combined field measurements and video analysis.Rock velocities and rebound heights were calculated by a frame-per-frame video analysis.
For the model evaluation mean and maximum of rock velocities and rebound heights were available at both sites.The data set further included the distribution of run out distances (ROD) and the jump lengths near the camera positions at both sites (see Fig. 1).For Site 1 the data set was extended by rebound heights and velocities of individual jumps near the cameras, for Site 2 also the distribution of the heights of tree hits on the bole was available.

Single 2-D rockfall trajectories in Steg and Bad Ischl
In addition to the Vaujany data set we used two rockfall trajectories of real rockfall events for model evaluation (Berger et al., 2004).In both cases a 2-D-profile was derived by surveying impact craters of single rockfall events.Additionally, the starting points of the two rockfall trajectories were identified and the size and mass of the rocks were measured.The first trajectory was recorded in Steg, Switzerland, after a rockfall event in March 2003 on a slope with an average inclination of 31 • .The rock had a diameter of 1.84 m and a mass of 9100 kg and stopped after 37 ground impacts and 441 m (planimetric) run out distance.The slope is mainly covered by pasture and is crossed by a road and a narrow rock outcrop.The surveyed data includes qualitative characterization of surface properties along the whole trajectory (see Fig. 2 and Table 2).In Bad Ischl, Austria, in summer 2004 a rock with a diameter of 0.96 m and a mass of 1250 kg hit and damaged a rockfall net protecting a parking lot.Here, the slope with an average inclination of 40   2.

Test site Vaujany
The coefficients of restitution for each surface type from the Vaujany test sites were extracted from the literature (e.g., Azzoni et al., 1995;Azzoni and de Freitas, 1995;Chau et al., 2002;Schweigl et al., 2003;Heidenreich, 2004) (Table 1).To mimic the release of the rocks from the shovel of the excavator in the simulations the rocks were started from a height of 1.5 m above ground with a speed of 3 ms −1 and an angle of −30 • (relative to the horizon).
To account for differences in the shape of the rocks we varied the tangential coefficient of restitution r t randomly by adding a value from the range [−0.05 to 0.05] assuming a uniform probability distribution.An independent stochastic variation of both coefficients may lead to unrealistic jump angles.Due to lacking information on a possible cross correlation between r t and r n we decided to vary only r t due to its higher relevance for the calculation steps during rock rebounds.
To achieve statistically robust results two different simulation series were run: (i) a total of 10 000 rocks with a diameter of 0.95 m were released from the original starting point, and (ii) to gain insight into the variability of extreme values a series of 100 simulations, each consisting of 100 rocks with a diameter of 0.95 m, were also performed.
For the forested Site 2 we used the measured trees on their mapped positions as well as tree species, diameter and tree height.For the simulation series the dynamic components of the forest model were deactivated.Contrary to the in situ experiment, where the forest structure changed potentially with every released rock, simulated rocks always faced the same trees.For each simulated rockfall trajectory mean, minimum and maximum velocity, rebound heights, jump lengths, tree contacts as well as run out distance were recorded and compared to experimental data.

2-D trajectories in Steg and Bad Ischl
For the simulation of the 2-dimensional trajectories, the rockfall model was operated in "2-D"-mode, where the trajectories of the simulated rocks are forced along the slope profile.Also the tree hit module was disengaged because no tree hits occurred along the trajectories.Quantitative coefficients of restitution were assigned to each mapped surface cover type along the trajectories using information from in situ field tests on representative ground materials from the literature (Azzoni et al., 1995;Azzoni and de Freitas, 1995;Schweigl et al., 2003;Heidenreich, 2004).Due to missing information about the exact starting conditions all trajectories were initiated in rolling mode with a velocity of 1 ms −1 .
To account for the stochasticity of the modelled processes 10 000 rocks were simulated for each 2-D site and analyzed for the run out distance of each trajectory as well as the energy at the position of the rockfall net in the Bad Ischl case.Additionally, the profiles were divided into segments of 10 m length along the slope and the mean, minimum and maximum length of jumps that ended in that respective section were recorded.

Initial simulation runs in Vaujany
Simulations were run at both sites in Vaujany according to the setup described in Sect.2.3.1 with the initial model version.A first inspection of simulated and observed rockfall parameters yielded the following results: for both sites in the upper part of the slope (labelled as "upper part" in Fig. 1) the jump lengths were generally too low whereas in the lower parts of the slope jump lengths were clearly overestimated at both sites with regard to both mean values and variability (Fig. 3; upper panel).The deviation between observed and simulated rebound height was smaller (Fig. 3; lower panel), but with a similar tendency to underestimate on the upper slope and overestimate in the lower parts of the slope.That simulated rebound heights were partly overestimated was also confirmed by the simulated heights of tree hits at Site 2 (Fig. 3).Furthermore, the simulated maximum values for rebound heights and rock velocities were found too high (not shown).
In addition, the distribution of the run out distances (ROD) differed substantially from observations at both sites (Fig. 4).Based on these findings the rebound algorithm, one of the central modules of the rockfall model (i.e. the calculation of energy dissipation during a ground contact and the resulting jump angles after ground contact), was re-analysed.The current model version (Woltjer et al., 2008) treats the re-bound process as an inelastic point contact based on impulsemomentum law.It is well known that this is a somewhat unrealistic and oversimplifying assumption, nevertheless this approach is frequently used in rockfall models (e.g.Stevens, 1998;Guzzetti et al., 2002;Dorren et al., 2004).
When applying repeatedly the current jump-rebound module on a simplified slope with a constant slope angle and constant coefficients of restitution the rock velocity increases or decreases exponentially (Eq.2).The actual values of the velocity increment or decrement per jump-impact cycle f gain mainly depends on slope angle and the values of the specific coefficients of restitution.
f gain = relative gain/loss in velocity per jump [−], v i = starting velocity of the ith jump [ms −1 ].The velocity reached at the i-th jump is proportional to f gain raised to the power of i (Eq.3).Maximum rebound height and jump length are proportional to the squared velocity and can therefore be estimated by Eqs. ( 4) and (5).Jump angles are constant and depend on slope angle and the ratio of the coefficients of restitution.Due to the exponential behavior, trajectories may start with a long series of very short jumps and a slow velocity gain but eventually reach unrealistic velocities and rebound heights if the slope is only long enough.However, this is in contrast to empirical evidence about real rockfall trajectories.For instance, on steep sites rocks typically gain energy very quickly and reach maximum velocity after a short distance (e.g.Dorren et al., 2005;Jahn, 1988).Furthermore, reported rebound heights usually do not exceed 2 m to 4 m while simulated rebound heights are much higher (e.g.Dorren and Berger, 2006;Stoffel et al., 2006; see also Fig. 3).Re-inspection of the model formulations focused on the amount of energy that is lost during rebound.Conceptually, the loss depends on the characteristics of the contact process itself and the amount of energy gain during the flight phase, which in turn is strongly influenced by the jump angle after rebound.
The rebound process is characterized as an elasto-plastic contact in many studies (e.g., Pfeiffer and Bowen, 1989), with elastic boulder response and plastic deformation of softer ground material.Johnson (1985) proposed to extend the widely applied point-like contact concept to plastic reaction by reducing the coefficients of restitution with increasing velocity proportional to v (−0.25) to account for deformation work, elastic wave propagation and fracturing during impact (Eq.6).
This relationship is supported by investigations of elastoplastic contact reaction (Wu et al., 2003(Wu et al., , 2005;;Hayakawa and Kuninaka, 2003;Heidenreich, 2004).Similar relationships are also found in studies based on field data from rock slopes (Wu, 1985;Pfeiffer and Bowen, 1989;Pfeiffer et al., 1993).The model extension by Johnson (1985) results in increased energy dissipation at high velocities, which is in agreement with field observations (e.g.Dorren et al., 2005;Jahn, 1988).However, rock velocity does not affect the jump angle after impact.Pfeiffer and Bowen (1989) and Pfeiffer et al. (1993) developed separate scaling factors for r n and r t from field data causing both coefficients of restitution to decrease with increasing impact velocity.This concept results in both a pronounced energy dissipation and slightly lower jump angles at high impact velocities.This is in accordance with Zinggeler (1989) and Heidenreich (2004) who found for elasto-plastic contacts no strict geometric relationship between impact and jump angles due to various influential factors such as friction or material strength.

Rebound algorithm
Based on the analysis of the original rebound model as presented above the contact model was extended following the proposal by Johnson (1985) as it refers to velocity as the most influential parameter of the jump-rebound cycle (Eq.6).Generally, this approach is supported by many studies over various materials and settings (Wu, 1985;Wu et al., 2003Wu et al., , 2005;;Hayakawa and Kuninaka, 2003;Heidenreich, 2004).
While still applying the literature values for impact velocities below 10 ms −1 , both the normal and the tangential coefficients of restitution were reduced at higher velocities according to Eqs. ( 7) and ( 8).This threshold was chosen because values for coefficients of restitution are usually based on rebound field experiments with rock velocities of about 10 ms −1 or below (Broili, 1973;Azzoni and de Freitas, 1995;Chau et al., 2002).
r literature = reported coefficient of restitution in literature (tangential and normal); see Tables 1 and 2 [−], r = resulting coefficient of restitution (tangential and normal) [−], k = linear factor, v=impact velocity [ms −1 ].Furthermore, the calculation of jump angles after rebound was modified.The analysis of the jump-rebound cycle already indicated that the geometric approach to calculate the jump angle may not be sufficient to reliably represent the complex processes during actual rebounds.The real-size rockfall experiments from Vaujany provided data on rock velocity after rebound and rebound height for a series of individual jumps on the upper slope of the unforested Site 1.Using this information for each of these jumps the jump angle after rebound was back-calculated.A power function (Eq.9) was fitted to back-calculated rebound angle and rock velocity data (compare Fig. 5).
The updated rebound model uses this relationship to calculate the jump angles after each rebound.This approach aims at providing realistic jump angles while preserving the ability to consider the dampening characteristics of different surface types via the coefficients of restitution.After implementing these modifications a second set of simulations were run.Simulated velocities and rebound heights, as well as the distribution of jump lengths were much closer to observed values.However, the simulated distribution of run-out distances deviated still substantially from observations (Fig. 4).Particularly the high frequency of stops in the middle section of the unforested Site 1 in Vaujany (Section B in Fig. 1) was still not satisfactorily captured by the improved model.This section, where almost one third of the rocks stopped despite high rock velocities, is characterized by the presence of debris large enough to act as obstacles.We concluded that rocks are stopped randomly by obstacles on the surface, a process which is not accounted for by the standard rebound algorithm.

Stopping algorithm
The hypothesis in improving model behavior was that the probability for rock stops at high velocities depends on the relation of (i) the size and distribution of obstacles on the surface of a given surface type and (ii) the size of the falling rock.To establish such relationships we developed a separate software module external to the PICUS model and ran a large number of simulations of boulder impacts on virtual surfaces of block deposits.The virtual surface of block deposits is generated by imitating a natural deposition process where added blocks tend to accumulate at already present deposits.The characteristic sizes of the deposited blocks are either pre-defined or drawn randomly from observed rock size distributions.This process roughly resembles a rockfallrelated talus generation.After reaching a defined degree of surface coverage the process is stopped (see Fig. 6 for an example).
Starting with randomly distributed angles between 10 • and 20 • relative to the surface, i.e. a realistic range of impact angles for real rockfall events, rocks are thrown from random positions above the surface.A rock is considered "stopped" if an obstacle at the surface is centrally hit.This is the case when the angle between the hit point and the direction vector of the falling rock is below 25 • (Fig. 7).If no obstacle is hit a simplified rebound is simulated assuming equality between incoming and outgoing angle.Finally, the "stopping probability" for a certain combination of surface type and falling rock size is defined as the ratio of stopped to total number of events (Eq.10).
p stop = probability that a rock is stopped during ground contact, N stopped = number of rocks that were stopped, N total = total number of rocks thrown onto surface.
During trajectory calculations within the Rock'n'Roll rockfall model Eq. ( 10) is used to calculate p stop for each ground impact.If the kinetic energy of a falling boulder exceeds zero, a uniform random number R rs [0-1] decides whether a falling rock comes to a halt (a boulder is stopped if p stop > R rs ).To that end the macro-roughness of the generated virtual surface types was visually compared to surface types observed in the field by mapping personnel.First tests yielded promising results and should be corroborated by an extended quantitative comparison of simulated and observed surface roughness in the future.

Vaujany
After implementing the new model algorithms the complete set of simulations was repeated.For 3-D simulations -especially when the forest is considered explicitly -run-out distances are of particular interest.Regarding the simulated distribution of run-out distances the improved model was able to reproduce the observed pattern well (Fig. 4).Especially for the unforested Site 1 the effect of the rockstop submodule is notable for the region between 225 m and 300 m (Section B in Fig. 1).In this section the slope inclination is about 38 • , comparable to the section above, but the surface is covered with bigger blocks resulting in a higher probability of stopping impacts.
Table 3 shows that mean and maximum velocities as well as the average number of tree hits as simulated by PICUS Rock'n'Roll match the observations well.While the model only slightly underestimates the average rebound heights, the simulated maximum rebound heights deviate still substantially from observed values, especially on the forested Site 2. The average of the maximum rebound heights from 100 repetitions of a simulation series with 100 rocks overestimates the maximum rebound height by factor three.However, in the majority of simulations, the maximum rebound heights are substantially lower (see average value and value of the 90th percentile) and extremely high jumps were simulated only on few specific spots along the slope.The general plausibility of simulated rebound heights is also indicated by the good match of the height at which trees were hit on Site 2 (Fig. 3, bottom left corner).
Figure 8 shows the distribution of the number of tree hits per rock on the forested Site 2. The error bars indicate the extreme values for 100 simulation series with 100 rocks each.The simulated number of rocks that hit only one tree or no tree at all is lower than the observed value, while the number of rocks that hit two to four trees is overestimated by the model.A possible explanation for this deviation is that during the Vaujany field experiment the forest especially near the starting point was thinned out by previously released rocks, thus reducing the likelihood of early tree hits.
The match of simulated mean and standard deviation of jump lengths for the upper and lower part of the slopes improved considerably with the revised model (Fig. 3).Compared to results of the original model version the simulated jumps of the updated model in the upper part of the slope are now longer while the jump lengths in the downslope sections are shorter and exhibit less variation (Fig. 3).

2-D trajectories in Steg and Bad Ischl
The 2-D trajectories are single events which are compared to the range of outcomes from a stochastic model.The simulated run out distance in Steg with 453 m±9 m (standard  deviation) is in agreement with the observed value of 441 m.In Bad Ischl more than 99.9% of the simulated rocks reach the rockfall net with an average kinetic energy of 349 kJ±132 kJ (90th percentile=570 kJ).Keeping in mind the stochastic nature of observed jump lengths resulting from a singular rockfall event, the observed jump lengths fit well within the distribution of simulated jump lengths (see Fig. 9).

Discussion
In this paper a comprehensive model validation experiment with a 3-D rockfall and forest model was presented.Focus of the validation was on the rockfall model.The study features an iterative cycle of testing and improving the model (e.g.Vanclay and Skovsgaard, 1997).A first set of simulation runs at the sites in Vaujany indicated major deficiencies of the initial model version with regard to simulated rock ve-locities, rebound heights and run-out distances.Attempts to improve the match of simulated and observed data through calibration of the coefficients of restitution were not successful.Based on the re-analysis of the rebound algorithm two major changes in the rebound algorithm were implemented: (i) an increased damping during ground impacts with increasing rock velocity based on Johnson (1985), and (ii) a new approach to calculate jump angles after rebound based on empirical jump data.The effect of increased damping for ground contacts at higher rock velocities due to a larger effect of plastic deformation is not only suggested by literature but has been successfully implemented in existing rockfall models (e.g.Pfeiffer and Bowen, 1989).While the complex impact process is still not fully understood, the magnitude of plastic deformation is assumed to depend on ground material characteristics (mainly strength) and impact characteristics such as impact angle and contact pressure.However, in the current study the relationship of velocity and the coefficients of restitution is the same for all surface types.Additionally, an empirical relationship to estimate jump angles after rebound was derived.This modification resulted in a more realistic model behavior but also in a weakened physical foundation of the rebound algorithm.We used specific data from one of the Vaujany experimental sites itself, thus rendering the Vaujany data no longer an entirely independent data set for model evaluation.However, detailed empirical data on rockfall processes is sparse, and excluding the Vaujany data from model development deemed unduly restrictive.We are aware that the empirical relationship of velocity and jump angle based on limited data from one specific site is not a generally valid relationship.However, the modified model variant combining Johnson's approach with the empirical relationship worked well for both sites in Vaujany as well as for the two trajectories in Steg and Bad Ischl.Considering the quite diverse conditions represented by these four sites we conclude that the chosen approach has some potential which should be further explored.More empirical data on individual jumps from different surface types are required to test the generality of the approach and to provide insight into the sensitivity of simulated trajectories to variation in the rebound algorithm.In general, the necessity to calibrate rockfall models for each individual site and application may severely hamper model applicability.For instance, for assessments at larger scales such calibration procedures are very likely not feasible due to missing calibration data, thus calling for generalizable approaches.
A major model improvement was the implementation of a probabilistic stopping algorithm for falling rocks.Based on a conceptual model of the interaction of a falling rocks with boulder obstacles on the ground virtual data were generated with an external software module to derive a general meta model of the process.The advantage of the presented approach is that the virtual surface types used to establish stopping probabilities for falling boulders of specific size can be compared to surface types used for site mapping activities.Thus, a direct link can be established between input parameters used in generating the virtual surface types and observable data (e.g., block cover percentage, size distribution of deposited blocks).
Generally, the updated model showed good agreement with observed rockfall parameters, particularly with regard to rock velocities, jump length and run-out distances.Worthwhile to note is the persistent overestimation of maximum rebound heights for the Vaujany simulations despite the relatively small range of additionally applied variability in the tangential coefficient of restitution r t .However, especially on the forested Site 2 large rebound heights occurred on only few spots along the slope.A possible reason is that Rock'n'Roll calculates jump angles relative to the slope angle at the impact point.Large rebound heights can thus arise from an overestimated variability of the inclinations along the slope as a result from artifacts of the digital elevation model (Agliardi and Crosta, 2003) or the transformation to a TIN.A possible further model extension would be the consideration of soft impact with penetration of the boulder into the ground on loose soils and scree deposits.This feature could be implemented by smoothing the point inclinations along an estimated contact zone.
Both single event 2-D trajectories were contained within the range of simulated stochastic outcomes.Compared to comprehensive 3-D data sets as available in Vaujany, 2-D profiles from post-hoc analysis of individual rockfall trajectories provide less opportunity to explore and validate model behavior.However, they proved to be a readily available and efficient complementary element in our model evaluation.
In the presented validation experiments the forest dynamics model was deactivated as the rockfall processes featured in the evaluation data took place within one time step of the forest model.As already outlined in Woltjer et al. (2008) a fully coupled forest and rockfall model with feedback of rock hits on the vitality and subsequently stability and mortality risk of trees is conceptually possible within this framework.Based on data of the Austrian National Forest Inventory Vospernik (2002), for instance, showed that Norway spruce trees with rockfall wounds had an increased probability of death.However, rock fall frequency as a prerequisite for such an approach cannot be readily inferred from standard data sources, limiting a fully coupled simulation approach.

Conclusions
The current model version produces reliable results over a wide array of conditions for rockfall processes as well as for forest dynamics (e.g.Seidl et al., 2008Seidl et al., , 2005)).It can thus be applied to analyze protective effect of forest vegetation over time periods of several decades under natural forest development or different forest management regimes.The rockfall model itself is designed for high computing performance and is able to calculate approx.1000 trajectories per second (using standard hardware from 2007) which allows to efficiently extend the spatial scale of application.The coupled rockfall and forest model can simulate an area of about 40 hectares which is sufficient for most protection forest management projects.The evaluation results presented in this contribution generally support the approach and increase the credibility of the rockfall module.Hence, we are confident that the presented model can be used for operational assessments of the protective function of forests against rockfall.

Fig. 1 .
Fig. 1.Profiles of the test sites in Vaujany starting from the release points of the rockfall experiments along the steepest path.Left: unforested Site 1, right: forested Site 2. The camera equipped test sections extend from the top of the slope down to the upper forest road.The observed run-out zones extend further (maxima indicated by triangles).

WFig. 3 .
Fig. 3.Comparison of observed and simulated (initial, improved) rockfall parameters from both experimental sites in Vaujany.Jump lengths and rebound heights are taken from the upper and the lower part of the slope near the camera positions, height of tree hits from the entire slope at Site 2. Shown are mean values and standard deviation of 100 simulated rocks with 100 replicates each.

Fig. 4 .
Fig. 4. Percentage of rocks that reach a certain distance at the unforested Site 1 (left) and forested Site 2 (right) in Vaujany.The filled area denotes the observed values, the bold line shows the simulation results (N =10 000 rocks) of the final model version.The dash-dot lines indicates the initial simulation runs, the dashed lines results with improved rebound algorithm but without the rockstop submodule.
starting velocity of the i-th jump, f gain = relative velocity gain/loss per jump [−], H i = rebound height for i-th jump [m], D i = jump length (planimetric) for i-th jump [m].

Fig. 5 .
Fig. 5. Relationship of rebound angle and velocity just after rebound for 63 recorded rebounds on the unforested Site 1 in Vaujany.

Fig. 8 .
Fig.8.Observed and simulated frequency of rocks in categories for "number of tree hits per rock" for the forested Site 2 in Vaujany.The simulated tree hits are represented by the mean value of 10 000 simulated rocks, the error bars indicate the minimum/maximum range of 100 simulated rocks with 100 replicates each.

Fig. 9 .
Fig. 9. Simulated and observed jump lengths for the profiles in Steg (left) and Bad Ischl (right).Dots denote observed jump lengths.The bold line indicates mean jump length and the dashed line the minimum and maximum jump length of 10 000 simulations.

Table 1 .
Coefficients of restitution (r n = normal, r t = tangential coefficient) used in the simulation series for the Vaujany sites.

Table 2 .
Coefficients of restitution (r n = normal, r t = tangential) used in the 2-D simulation series Steg and Bad Ischl.

Table 3 .
Summary of observed and simulated (with improved rebound algorithm) rockfall parameters for both sites in Vaujany.The columns denoted with N = 100 × 100 present the mean and standard deviation of 100 simulation series with 100 rocks each.