Application of Unmanned Aerial Vehicle Data and Discrete Fracture Network Models for Improved Rockfall Simulations

In this research, we present a new approach to define the distribution of block volumes during rockfall simulations. Unmanned aerial vehicles (UAVs) are utilized to generate high-accuracy 3D models of the inaccessible SW flank of the Mount Rava (Italy), to provide improved definition of data gathered from conventional geomechanical surveys and to also denote important changes in the fracture intensity. These changes are likely related to the variation of the bedding thickness and to the presence of fracture corridors in fault damage zones in some areas of the slope. The dataset obtained integrating UAV and conventional surveys is then utilized to create and validate two accurate 3D discrete fracture network models, representative of high and low fracture intensity areas, respectively. From these, the ranges of block volumes characterizing the in situ rock mass are extracted, providing important input for rockfall simulations. Initially, rockfall simulations were performed assuming a uniform block volume variation for each release cell. However, subsequent simulations used a more realistic nonuniform distribution of block volumes, based on the relative block volume frequency extracted from discrete fracture network (DFN) models. The results of the simulations were validated against recent rockfall events and show that it is possible to integrate into rockfall simulations a more realistic relative frequency distribution of block volumes using the results of DFN analyses.


Introduction
Rockfalls are a major hazard to persons and property, especially in proximity of infrastructure such as roads, railways and housing. In this research, we analyzed the area of the Scanno Lake, at the foot of the Mount Rava. The area is characterized by deposits of paleolandslide and more recent landslide/rockfall events. During such recent rockfall events, rock blocks of volume varying between ca. 0.01 to 500 m 3 detached from the SW flank of the Mount Rava and travelled up until the base of the slope as well as in proximity of infrastructures such the Frattura road, the Frattura Vecchia village and the Frattura graveyard ( Figure 1). During this study, this section of the mountain was thoroughly investigated through conventional and remote sensing surveys, with the final objective to define an innovative approach based on the integration of structural geomechanical and rockfall trajectory analyses, and to provide improved understanding of the rockfall hazard posed by the geological setting. The steps involved in the analysis of rockfalls usually include the survey and characterization of rock outcrops, kinematic assessment and/or engineering classification of rock masses and, where necessary, stability and rockfall runout simulations.
The survey and characterization of rock outcrops allow the acquisition of information about geometrical and physical characteristics of rock masses, such as rock strength, slope and discontinuity attitude, discontinuity spacing, persistence, roughness etc. [1]. Such parameters are usually collected through conventional engineering geological (geomechanical) surveys that, in the case of inaccessible or high steep slopes, are sometimes combined with more innovative remote sensing techniques. The advent of these new remote sensing technologies for the survey of geological features has led to step-change increases in the quality of data available for slope/geomechanical studies. Laser scanning (LS) and digital photogrammetry (DP) have been the most widely used remote sensing techniques for landslide studies and characterization [2][3][4][5][6]. Lato et al. [2] showed how to improve the use of LS data for the automated structural evaluation of discontinuities in rock slopes, while Francioni et al. [4] illustrated the use of DP for the characterization and stability analysis with limit equilibrium methods of a coastal area in Cornwall (UK). Bonneau and Hutchinson [5] showed the use of terrestrial laser scanning for the characterization of a cliff-talus system in the Thompson River Valley (British Columbia, Canada) and Kromer et al. [6] developed an automated fixed-location time lapse photogrammetric rock slope monitoring system. A critical overview of some of the limitations of terrestrial DP and LS when The steps involved in the analysis of rockfalls usually include the survey and characterization of rock outcrops, kinematic assessment and/or engineering classification of rock masses and, where necessary, stability and rockfall runout simulations.
The survey and characterization of rock outcrops allow the acquisition of information about geometrical and physical characteristics of rock masses, such as rock strength, slope and discontinuity attitude, discontinuity spacing, persistence, roughness etc. [1]. Such parameters are usually collected through conventional engineering geological (geomechanical) surveys that, in the case of inaccessible or high steep slopes, are sometimes combined with more innovative remote sensing techniques. The advent of these new remote sensing technologies for the survey of geological features has led to step-change increases in the quality of data available for slope/geomechanical studies. Laser scanning (LS) and digital photogrammetry (DP) have been the most widely used remote sensing techniques for landslide studies and characterization [2][3][4][5][6]. Lato et al. [2] showed how to improve the use of LS data for the automated structural evaluation of discontinuities in rock slopes, while Francioni et al. [4] illustrated the use of DP for the characterization and stability analysis with limit equilibrium methods of a coastal area in Cornwall (UK). Bonneau and Hutchinson [5] showed the use of terrestrial laser scanning for the characterization of a cliff-talus system in the Thompson River Valley (British Columbia, Canada) and Kromer et al. [6] developed an automated fixed-location time lapse photogrammetric rock slope monitoring system. A critical overview of some of the limitations of terrestrial DP and LS when dealing with high steep rock slopes was presented by Sturzenegger and Stead [7]. Many of the limitations discussed by Sturzenegger and Stead [3] have now largely been overcome by the increasing use of aerial platforms such as Unmanned Aerial Vehicles (UAV). The introduction of such platforms has dramatically improved the application of these systems [8,9], making DP and LS even drainage basin, between the Montagna Grande and Mt. Genzana ridges. It is in one of the highest average elevation areas in the entire Central Apennines. Radiocarbon dating of soil samples collected from the paleolandslide debris accumulation suggests an age of c. 12,800 years [45].
The Scanno paleolandslide has been investigated by several scholars over the last few decades. Nicoletti et al. [46] described the debris accumulation area and deposition mechanisms. Bianchi-Fasani et al. [45] and Della Seta et al. [47] focused their study on the landslide failure mechanisms, describing the Scanno landslide as a slow-moving rock avalanche where the bedding planes represent the sliding surface. More recently, Francioni et al. [48] proposed a new interpretation, suggesting the landslide scar is controlled by low-and high-angle normal faults associated with the Difesa-Monte Genzana-Vallone delle Masserie (DMG) fault zone. Their research highlighted how the high-angle faults (SW dipping, F2-F3) and two joint sets (dipping toward SW and SE) represent the backscarp surfaces and lateral release surfaces of the Scanno paleolandslide, respectively.
Geologically, the study area is characterized by Jurassic-Paleogene marine limestone, with very thick beds of calcarenites, and marly limestone with thin clayey marly layers (Mount Genzana, MG, unit). The middle-lower part of the valley is formed by pelitic-arenaceous siliciclastic rocks (Neogene foredeep deposits). Quaternary clastic continental deposits (slope breccia deposits, alluvial fan conglomerate) largely cover the bedrock [49]. On the SW flank of the Monte Rava, within the Scanno landslide scar, the bedrock is specifically covered by talus slope and cone deposits resulting from the recent degradation of the landslide crown. Figure 2 shows a simplified geological map of the study area, after Francioni et al. [48].

Study Area
Investigations were undertaken in the SW flank of the Monte Rava (Abruzzi Apennines, Italy). The SW flank of the Monte Rava is characterized by the enormous scar left by the famous Scanno paleolandslide, which dammed the Tasso River and formed Scanno Lake, one of the most famous examples of a naturally dammed lake in Italy (Figure 1). The area is located in the Sagittario River drainage basin, between the Montagna Grande and Mt. Genzana ridges. It is in one of the highest average elevation areas in the entire Central Apennines. Radiocarbon dating of soil samples collected from the paleolandslide debris accumulation suggests an age of c. 12,800 years [45].
The Scanno paleolandslide has been investigated by several scholars over the last few decades. Nicoletti et al. [46] described the debris accumulation area and deposition mechanisms. Bianchi-Fasani et al. [45] and Della Seta et al. [47] focused their study on the landslide failure mechanisms, describing the Scanno landslide as a slow-moving rock avalanche where the bedding planes represent the sliding surface. More recently, Francioni et al. [48] proposed a new interpretation, suggesting the landslide scar is controlled by low-and high-angle normal faults associated with the Difesa-Monte Genzana-Vallone delle Masserie (DMG) fault zone. Their research highlighted how the high-angle faults (SW dipping, F2-F3) and two joint sets (dipping toward SW and SE) represent the backscarp surfaces and lateral release surfaces of the Scanno paleolandslide, respectively.
Geologically, the study area is characterized by Jurassic-Paleogene marine limestone, with very thick beds of calcarenites, and marly limestone with thin clayey marly layers (Mount Genzana, MG, unit). The middle-lower part of the valley is formed by pelitic-arenaceous siliciclastic rocks (Neogene foredeep deposits). Quaternary clastic continental deposits (slope breccia deposits, alluvial fan conglomerate) largely cover the bedrock [49]. On the SW flank of the Monte Rava, within the Scanno landslide scar, the bedrock is specifically covered by talus slope and cone deposits resulting from the recent degradation of the landslide crown. Figure 2 shows a simplified geological map of the study area, after Francioni et al. [48]. The main DMG fault zone characterizes the area forming the tectonic contact between the MG unit (limestone and marly limestone in the footwall) with siliciclastic deposits in the hanging wall. The DMG fault zone consists of several faults forming horst structures of different scale and The main DMG fault zone characterizes the area forming the tectonic contact between the MG unit (limestone and marly limestone in the footwall) with siliciclastic deposits in the hanging wall. The DMG fault zone consists of several faults forming horst structures of different scale and sigmoidal geometries [48]. The bedrock shows a general NW-SE trend with irregular attitude related to the fault zones.
The occurrence of the Scanno landslide deeply modified the geometry of the SW flank of the Monte Rava, forming a concave and steep scar (with a maximum elevation a.s.l. of 1860 m at the top of Monte Rava) and a rugged accumulation area. The landslide accumulation area features a gentle slope in the upper part, with evidence of slump blocks (with back tilted faces) in the lower part along the valley bottom and upslope on the opposite valley side [48]. The bottom of the valley is represented by a flat alluvial area. The area is also characterized by fluvial-alluvial fan deposits (at the bottom of the valley) and by deposits of more recent landslide/rockfall events.

Materials and Methods
In this research, we develop a new approach involving the integration of discrete fracture network (DFN) simulations to support and improve rockfall analyses. This is based on a statistical analysis of fractures and derived DFNs to define the ranges of block volumes within a hypothetical/virtual rock mass. From the conventional geomechanical surveys, it was possible to calculate the orientation, Fisher constant K (a measure of the degree of discontinuity poles clustering), spacing and persistence of each discontinuity set. The 3D models obtained from UAV surveys were utilized to improve the data gathered from conventional surveys and to develop several sampling windows. These were then analyzed through the recently developed freeware software FracPaq [50], which allowed definition of the ranges of fracture intensity for each discontinuity set. Through these data, 3D DFN models of the slope and the range of rock block volumes characterizing the rock mass were established. These ranges were subsequently used as input parameters in 3D rockfall trajectory simulations through the software Rockyfor3D (RF3D), a GIS-based model based on the integration of both statistic-and process-based approaches [51].
Two different types of rockfall analysis were performed, rockfall analysis 1 and 2. During the rockfall analysis the released rock blocks class had a Gaussian uniform distribution, defined by fixed presets available in RF3D. Once the results gathered from this analysis were calibrated against a validation set of end points (arrest locations of rock blocks) mapped from aerial RGB orthophotographs, a novel method to define the distribution of rock block volumes was developed and tested in rockfall analysis 2. This analysis was based on the use of a more robust nonuniform block volume distribution using the relative block volume frequency extracted from DFN models.

Conventional and DP (UAV-Based) Geomechanical Surveys
Conventional geomechanical surveys were performed at the toe of the Scanno paleolandslide scar (SW flank of the Mount Rava) to determine the characteristics of the main discontinuity sets. Due to the difficulties in accessing most of the outcrops in the vicinity of the scar zone, only five small geomechanical scanlines were carried out ( Figure 3A). By flying a UAV, it was possible to reach previously inaccessible areas and obtain topographical data representative of the entire slope. Three main areas within the landslide scar area were analyzed (UAV 1, UAV 2 and UAV 3, Figure 3A). The slopes on the three areas, UAV 1-3, present different orientations that allowed the reduction of orientation bias during the measurements of discontinuity attitude. Dip and dip direction of the slopes within these areas are 65 • /249 • for UAV 1, 60 • /168 • for UAV 2 and 56 • /194 • for UAV 3. The drone used for the surveys was equipped with a camera with the following characteristics: 12.4 megapixel camera resolution; 2.8 mm focal length; 6.16 mm wide and 4.62 mm high camera sensor size. In order to achieve a high-resolution 3D model of the outcrops, photographs were taken from an average distance from the slope of ca. 35 and 50 m in locations with complex geological conditions (very high fracture intensity areas associated with thin bedded marly limestone or fracture corridors related to fault damage zones, such as area UAV 1 in Figure 3A) and in less fractured areas with very thick bedded calcarenites (UAV 2 and 3 in Figure 3A), respectively. The UAV DP survey was undertaken through multiple vertical photographic strips, flying at an ascending and descending speed of 5.0 km/h. Figure 3B shows an example flight plan for the area covered by UAV 3. To guarantee a vertical overlap of ca. 80%, photographs were acquired every 2 seconds when the distance from the slope was ca. 35 m, and every 4 seconds when the distance was 50 m. The distance between each vertical strip was set to ca. 5 m to ensure a lateral overlap between strips of ca. 80%.
Remote Sens. 2020, 12, 2053 6 of 28 seconds when the distance from the slope was ca. 35 m, and every 4 seconds when the distance was 50 m. The distance between each vertical strip was set to ca. 5 m to ensure a lateral overlap between strips of ca. 80%. The photographs gathered from UAV surveys were processed using Agisoft Metashape [52] to create the 3D point cloud models of the outcrops. The freeware software CloudCompare [53] was then utilized to manage the point cloud data and define the main discontinuity sets in the rock mass and their respective spacing and persistence. For verification and validation purposes, the remote sensing rock discontinuity dataset was compared with the data measured in conventional contact scanline surveys. Other geomechanical parameters, such as joint roughness (JRC), joint compressive strength (JCS), joint alteration and aperture, were defined during conventional scanline surveys (G1-G5 in Figure 3) and integrated with the UAV data.

Sampling Windows and 3D DFN
The UAV-extracted 3D models were also used to create sampling windows in different slope areas. From these observations, we derived the fracture intensity, expressed as the length of fractures per unit area (P21) [54]. The discontinuities identified in the sampling windows ( Figure 4A) were imported into the FracPaq software ( Figure 4B), where it was possible to calculate P21. As opposed to conventional methods used for the calculation of P21, which considers the total length of fractures over the total area of the sampling windows, FracPaq can define multiple P21 values for every sampling window and calculate an average value and its standard deviation. This was The photographs gathered from UAV surveys were processed using Agisoft Metashape [52] to create the 3D point cloud models of the outcrops. The freeware software CloudCompare [53] was then utilized to manage the point cloud data and define the main discontinuity sets in the rock mass and their respective spacing and persistence. For verification and validation purposes, the remote sensing rock discontinuity dataset was compared with the data measured in conventional contact scanline surveys. Other geomechanical parameters, such as joint roughness (JRC), joint compressive strength (JCS), joint alteration and aperture, were defined during conventional scanline surveys (G1-G5 in Figure 3) and integrated with the UAV data.

Sampling Windows and 3D DFN
The UAV-extracted 3D models were also used to create sampling windows in different slope areas. From these observations, we derived the fracture intensity, expressed as the length of fractures per unit area (P 21 ) [54]. The discontinuities identified in the sampling windows ( Figure 4A) were imported into the FracPaq software ( Figure 4B), where it was possible to calculate P 21 . As opposed to conventional methods used for the calculation of P 21 , which considers the total length of fractures over the total area of the sampling windows, FracPaq can define multiple P 21 values for every sampling window and calculate an average value and its standard deviation. This was performed by calculating the total length of discontinuities within 16 test circles ( Figure 4C). The use of 16 circles for each sampling window decreased the uncertainties associated with the fracture intensity computation and observed variance within each sampling window.
Remote Sens. 2020, 12, 2053 7 of 28 performed by calculating the total length of discontinuities within 16 test circles ( Figure 4C). The use of 16 circles for each sampling window decreased the uncertainties associated with the fracture intensity computation and observed variance within each sampling window. The P21 values were then used to develop a 3D DFN model using Move software [55]. The DFN was created using the discontinuity sets recorded during the surveys and, for each set, the Fisher value (K), the fracture length and aperture were used to establish the fracture network. In the software Move, the fracture intensity is input as P32 (fracture area per unit volume). Since the P32 is not readily available and cannot be measured directly from the surface, the DFN analysis within Move must be carried out using an iterative approach. An initial hypothetical value of P32 is used to develop a DFN model ( Figure 5A). Several cross sections of this model were then exported ( Figure  5B shows an example cross section) and the fracture traces intersecting the test planes were used to calculate P21, utilizing the same procedure illustrated in Figure 4 for in situ sampling windows ( Figure 5C). The average P21 value was then compared with the value gathered from the sampling windows and the P32 calibrated when the P21 extracted from the DFN cross sections matched the values measured from the UAV-extracted sampling windows. Once validated, it was possible to extract the range of rock block volumes and the relative block volume frequency (i.e., the percentage of blocks with specific volumes within the rock mass considered) from the DFN model.  The P 21 values were then used to develop a 3D DFN model using Move software [55]. The DFN was created using the discontinuity sets recorded during the surveys and, for each set, the Fisher value (K), the fracture length and aperture were used to establish the fracture network. In the software Move, the fracture intensity is input as P 32 (fracture area per unit volume). Since the P 32 is not readily available and cannot be measured directly from the surface, the DFN analysis within Move must be carried out using an iterative approach. An initial hypothetical value of P 32 is used to develop a DFN model ( Figure 5A). Several cross sections of this model were then exported ( Figure 5B shows an example cross section) and the fracture traces intersecting the test planes were used to calculate P 21 , utilizing the same procedure illustrated in Figure 4 for in situ sampling windows ( Figure 5C). The average P 21 value was then compared with the value gathered from the sampling windows and the P 32 calibrated when the P 21 extracted from the DFN cross sections matched the values measured from the UAV-extracted sampling windows. Once validated, it was possible to extract the range of rock block volumes and the relative block volume frequency (i.e., the percentage of blocks with specific volumes within the rock mass considered) from the DFN model. Remote Sens. 2020, 12, 2053 7 of 28 performed by calculating the total length of discontinuities within 16 test circles ( Figure 4C). The use of 16 circles for each sampling window decreased the uncertainties associated with the fracture intensity computation and observed variance within each sampling window. The P21 values were then used to develop a 3D DFN model using Move software [55]. The DFN was created using the discontinuity sets recorded during the surveys and, for each set, the Fisher value (K), the fracture length and aperture were used to establish the fracture network. In the software Move, the fracture intensity is input as P32 (fracture area per unit volume). Since the P32 is not readily available and cannot be measured directly from the surface, the DFN analysis within Move must be carried out using an iterative approach. An initial hypothetical value of P32 is used to develop a DFN model ( Figure 5A). Several cross sections of this model were then exported ( Figure  5B shows an example cross section) and the fracture traces intersecting the test planes were used to calculate P21, utilizing the same procedure illustrated in Figure 4 for in situ sampling windows ( Figure 5C). The average P21 value was then compared with the value gathered from the sampling windows and the P32 calibrated when the P21 extracted from the DFN cross sections matched the values measured from the UAV-extracted sampling windows. Once validated, it was possible to extract the range of rock block volumes and the relative block volume frequency (i.e., the percentage of blocks with specific volumes within the rock mass considered) from the DFN model.

Rockfall Simulations
In order to capture the complex and unpredictable behavior of a failing mass of rock from a slope, a number of rockfall trajectory analysis codes were developed throughout the years, and provide an established method to assess rockfall hazard [16,17,20,56,57]. Fundamental information needed to generate models that can accurately describe rockfall trajectories is the topographical representation of the slope, a high-accuracy digital twin, either in the form of a raster DEM or a vertical profile. In the case of GIS-based methods, digital topography serves the purpose of providing a virtual surface from which to compute the types of motion (freefall, bounce, roll or drag) of the falling rock blocks. Broadly speaking, rockfall modeling software treats impact theory in two different ways: the lumped mass (LM) approach versus the rigid body (RB) approach. The lumped mass approach considers the mass being concentrated in a single point, whereas the rigid body approach uses a defined geometry to model the rock block [18,19].
As mentioned above, in this research we adopted the GIS-based code Rockyfor3D. The software uses a three-dimensional rigid-body impact model (RB) that calculates trajectories of single, individually falling rocks with discrete geometry (RB). Rockyfor3D can be used for regional-, local-and slope-scale rockfall simulations [51]. The input parameters that define rock blocks are the release cell location, the rock density, the shape and volume (with the possibility to set a statistical variation range) and initial vertical velocity. The local slope surface roughness is represented by a parameter defined as height of a representative obstacle (MOH), expressed in m. This parameter does not add a local variation of height, but accounts for the uncertain nature of the bedrock cover and its geometry, as well as its mechanical properties. Typical MOH values, as suggested by Dorren [51], which are encountered by a falling rock are represented by three statistical classes, rg70%, rg20%, and rg10%. During each rebound calculation, the MOH value in a cell is randomly chosen from the three representative values according to their probabilities of occurrence [51]. Finally, the soil type is defined through a raster map identifying the type of bedrock exposed. Once the soil-type slope is defined, the normal (Rn) and tangential (Rt) coefficients of restitution are set for each position within the DEM; these values are responsible for setting the energy transfer functions (i.e., the inelastic impacts or the energy loss upon every impact).
In this research, the simulations were performed using a 5 × 5 m resolution DEM. The locations of release points were selected based on the DEM's slope values. We defined a lower slope threshold of 55 • and pixels with a value exceeding that threshold were identified as candidate release cells. Slope roughness (MOH values) and soil type (Rn) were determined by field inspection and geomorphological analysis of high-resolution aerial images. The tangential coefficient of restitution (Rt) is automatically calculated by Rockyfor3D through the composition and size of the material covering the surface and the radius of the falling block itself [51]. The rock density, shape and volume were defined by combining field and remote sensed geomechanical data. In particular, the use of DFN models allowed definition of the possible ranges of block volumes. In Rockyfor3D, the volume of the blocks to be released has to be defined in each release cell; it can be uniform in all the cells or can vary within a predefined percentage (±5%, ±10%, and ±50%). This random variation is the same for all three block dimensions.
Two different rockfall analyses, rockfall analysis 1 and 2, were carried out. During rockfall analysis 1, to verify the reliability of the input parameters, we performed initial calibration tests utilizing a uniform ranges of block volume representative of DFN 1 (0.008-0.1 m 3 ) and eight uniform ranges of block volumes representative of DFN 2 (0.6-1.5, 1.6-3, 3.1-7.5, 7.6-20, 21-55, 56-120, 121-260, 261-585 m 3 ). These ranges were decided in relation to the total range of block volume gathered from the DFN 1 and 2 analyses and the ability to vary a predefined volume of a maximum percentage of 50%. Validation of the rockfall model was obtained iteratively by comparing the rockfall simulation results with a map representing the distribution of endpoints for recent rockfall events. This map, shown in Figure 6, was developed in a GIS environment, integrating the analysis of high-resolution orthophotographs and field inspection. Over 600 blocks with volume over 1 m 3 were identified and digitalized in the map and divided using the same ranges of block volumes utilized for the calibration Remote Sens. 2020, 12, 2053 9 of 28 of the rockfall model ( Figure 6A). Examples of block identification are represented in Figure 6B-F, which highlights the shape and volume of blocks interpreted in different areas at the toe of the slope.
In particular, Figure 6B shows the biggest block identified in the study area, with a volume of ca. 1000 m 3 .
Remote Sens. 2020, 12, 2053 9 of 28 represented in Figure 6B-F, which highlights the shape and volume of blocks interpreted in different areas at the toe of the slope. In particular, Figure 6B shows the biggest block identified in the study area, with a volume of ca. 1000 m 3 .  After achieving a satisfactory result (rockfall analysis 1 showed good agreement with a map of recent rockfall events), a new method to input nonuniform ranges of block volumes, based on the block size distribution (relative block volume frequency) gathered from DFN analysis, was developed and implemented in the simulations in rockfall analysis 2. Once the relative frequency of potential block volumes was extracted from the DFN, the same distribution was used to define the block sizes in the release cells. For example, if the DFN analysis shows that 10% of blocks have a volume of 10 m 3 , and considering a total of 3000 release cells, 300 of these cells (randomly chosen) will be populated with 10 m 3 blocks. This procedure was carried out through the integrated use of GIS spatial analysis technique and spreadsheet software.
With this newly proposed approach, it was possible to perform a simulation with a customized relative frequency distribution of rock sizes and volumes, rather than using a standard Gaussian distribution (i.e., as used in rockfall analysis 1), hence utilizing a release hazard scenario based on the statistics extracted from the DFN model.
Considering that the choice of the release cells should be random, different simulations were carried out to verify possible scenarios. The determination/choice of the best scenario could be established by comparing the results of the simulation with a map of recent rockfalls (as was done in this research). When such a map is unavailable, all the scenarios then need to be considered potentially realistic.

Geomechanical Data and Fracture Analysis
Three main photogrammetric models were derived from the UAV surveys (UAV 1-3). Table 1 shows the details of the UAV surveys and derived point clouds. The area imaged using UAV 1 is the location with highest ground resolution (1.32 cm/pix) due to the geological characteristics and high fracture intensity. To achieve such resolution, the survey was conducted at 36 m from the slope using 1365 photographs. Images of areas UAV 2 and 3 were acquired from ca. 46 and 48 m from the slope, resulting in resolutions of 1.7-1.9 cm/pix. Figure 7 shows representative 3D models of areas UAV 1 ( Figure 7A The integration of geomechanical and DP data highlighted two main discontinuity sets, J1 and J2, NW and NE striking, respectively. Bedding planes, S0, show a variable orientation, with a generalized NW-SE orientation, and ranging from NE dipping to SW dipping in proximity to the normal faults outcropping in proximity of area UAV 1. In this area, the rock strata appear dragged towards the SW by the main SW dipping faults. Figure 7D shows the stereonet, highlighting the main discontinuity sets and the variability in the bedding plane orientations. The different dip and dip direction of slope faces on the three UAV areas reduced orientation bias during discontinuity analysis. In particular, the area UAV 1 (dip and dip direction of slope 65°/249°) was particularly suited to the study of set J2 and S0, while the areas UAV 2 and 3 (dip and dip-direction of slope 60°/168° and 56°/194°, respectively) for the study of J1 and S0. In general, it was found that the orientation of the discontinuity sets was constant in the different slope areas (with the exception of S0 in proximity of fault zones, as mentioned above).
Due to the geological characteristics of the study area, persistence and intensity of fracture sets vary within the landslide scar, especially in the proximity of thin bedded layers and/or fault zones. Due to this, it was decided to create two different DFN models, DFN 1 and DFN 2, representative of high fracture intensity domains (with high P21 values) and low fracture intensity domains (with low P21 values), respectively. The fracture intensity variation was evaluated through the analysis of several sampling windows created within the three investigated areas (UAV 1-3). The discontinuities identified in these sampling windows (Figures 8A and 9A) were analyzed through FracPaq (Figures 8B and 9B) in order to calculate the fracture intensity value, P21, of each discontinuity set (P21 was measured using the procedure illustrated in Section 3 and Figure 4). Fracture intensity values and length (persistence) of discontinuity sets given in Tables 2-4 were used to develop and constrain the 3D DFN models. In Table 2, the fracture intensity values are shown in terms of P21 extracted from sampling windows and DFN. Tables 3 and 4 show the fracture intensity (in terms of P32) and the discontinuity lengths, extracted from UAV sampling windows and included in the DFN model using mean and standard deviation or power law distribution (the option chosen was the one that best fit with the sampling window data). The results of this procedure, in both high and low fracture intensity areas, is illustrated in Figures 8 and 9, respectively. Figure 8 shows the discontinuities identified in high fracture intensity areas ( Figure 8A), their analysis in FracPaq ( Figure 8B) and the section extracted from the final DFN 1 and 2 ( Figure 8C). Figure 9 illustrates the The integration of geomechanical and DP data highlighted two main discontinuity sets, J1 and J2, NW and NE striking, respectively. Bedding planes, S0, show a variable orientation, with a generalized NW-SE orientation, and ranging from NE dipping to SW dipping in proximity to the normal faults outcropping in proximity of area UAV 1. In this area, the rock strata appear dragged towards the SW by the main SW dipping faults. Figure 7D shows the stereonet, highlighting the main discontinuity sets and the variability in the bedding plane orientations. The different dip and dip direction of slope faces on the three UAV areas reduced orientation bias during discontinuity analysis. In particular, the area UAV 1 (dip and dip direction of slope 65 • /249 • ) was particularly suited to the study of set J2 and S0, while the areas UAV 2 and 3 (dip and dip-direction of slope 60 • /168 • and 56 • /194 • , respectively) for the study of J1 and S0. In general, it was found that the orientation of the discontinuity sets was constant in the different slope areas (with the exception of S0 in proximity of fault zones, as mentioned above).
Due to the geological characteristics of the study area, persistence and intensity of fracture sets vary within the landslide scar, especially in the proximity of thin bedded layers and/or fault zones. Due to this, it was decided to create two different DFN models, DFN 1 and DFN 2, representative of high fracture intensity domains (with high P 21 values) and low fracture intensity domains (with low P 21 values), respectively. The fracture intensity variation was evaluated through the analysis of several sampling windows created within the three investigated areas (UAV 1-3). The discontinuities identified in these sampling windows (Figures 8A and 9A) were analyzed through FracPaq (Figures 8B and 9B) in order to calculate the fracture intensity value, P 21 , of each discontinuity set (P 21 was measured using the procedure illustrated in Section 3 and Figure 4). Fracture intensity values and length (persistence) of discontinuity sets given in Tables 2-4 were used to develop and constrain the 3D DFN models. In Table 2, the fracture intensity values are shown in terms of P 21 extracted from sampling windows and DFN. Tables 3 and 4 show the fracture intensity (in terms of P 32 ) and the discontinuity lengths, extracted from UAV sampling windows and included in the DFN model using mean and standard deviation or power law distribution (the option chosen was the one that best fit with the sampling window data). The results of this procedure, in both high and low fracture intensity areas, is illustrated in Figures 8 and 9, respectively. Figure 8 shows the discontinuities identified in high fracture intensity areas ( Figure 8A), their analysis in FracPaq ( Figure 8B) and the section extracted from the final DFN 1 and 2 ( Figure 8C). Figure 9 illustrates the discontinuities identified in low fracture intensity areas ( Figure 9A), their analysis in FracPaq ( Figure 9B) and the section extracted from the final DFN models representing low fracture intensity areas ( Figure 9C). It is possible to observe, following the validation procedure, that the intensity values are similar and the DFN model is representative of the fracture intensities measured on the slope.    Figure 10 shows the results of DFN 2. Figure 10A illustrates a graph of relative block volume frequency (%). Figure 10B shows the DFN model for the three discontinuity sets J1, J2 and S0, while Figure 10C shows the range of block volumes represented through the discretization of a 400 × 400 m rock mass (dimensions similar to the studied slope, colors represent different block sizes). It is possible to observe from Figure 10A that, in DFN 2, the minimum block volume is 0.6 m 3 , while the maximum can reach 1000 m 3 . However, the relative frequency of block volume can vary within this range, with the majority of the block volumes (ca. 70%) below 25 m 3 . Circa 20% of block volumes are between 25 and 100 m 3 , ~5% are between 100 and 200 m 3 and ~2% between 200 and 300 m 3 . Only 0.6% and 0.3% of blocks are made up of volumes up to 500 m 3 and 1000 m 3 , respectively.   Figure 10A illustrates a graph of relative block volume frequency (%). Figure 10B shows the DFN model for the three discontinuity sets J1, J2 and S0, while Figure 10C shows the range of block volumes represented through the discretization of a 400 × 400 m rock mass (dimensions similar to the studied slope, colors represent different block sizes). It is possible to observe from Figure 10A  The DFN 1 presents very different results ( Figure 11A-C) with more than 99% of block volumes below 0.1 m 3 ( Figure 11A), a minimum volume of 0.008 m 3 and a maximum volume slightly above 1 m 3 . This DFN can be considered representative of highly fractured areas associated with the DMG fault zone.  The DFN 1 presents very different results ( Figure 11A-C) with more than 99% of block volumes below 0.1 m 3 ( Figure 11A), a minimum volume of 0.008 m 3 and a maximum volume slightly above 1 m 3 . This DFN can be considered representative of highly fractured areas associated with the DMG fault zone.

Rockyfor3D Simulation Results
Rockfall simulations provide maps of the spatial distribution of end points of rock block trajectories. The results of these simulations are usually presented in the form of probability density functions showing the reach probability, the volume and number of blocks deposited, kinetic energy etc.
Slope roughness and soil type were initially obtained by studying the available geomorphological map, high-resolution orthophotographs and by field inspections. The MOH values were subsequently adjusted during the rockfall model calibration in relation to model response and comparison between rockfall model results, and location of rock blocks associated with recent previous rockfall events (map of Figure 6). Figure 12 highlights the soil type map obtained after the calibration process. Table 5 reports the MOH values, represented by three statistical classes, rg70, rg20 and rg10, assigned to each soil type.

Rockyfor3D Simulation Results
Rockfall simulations provide maps of the spatial distribution of end points of rock block trajectories. The results of these simulations are usually presented in the form of probability density functions showing the reach probability, the volume and number of blocks deposited, kinetic energy etc.
Slope roughness and soil type were initially obtained by studying the available geomorphological map, high-resolution orthophotographs and by field inspections. The MOH values were subsequently adjusted during the rockfall model calibration in relation to model response and comparison between rockfall model results, and location of rock blocks associated with recent previous rockfall events (map of Figure 6). Figure 12 highlights the soil type map obtained after the calibration process. Table 5 reports the MOH values, represented by three statistical classes, rg70, rg20 and rg10, assigned to each soil type.   Prior to running the simulations, the locations of release points (cells) were selected by setting a slope threshold, which was set to 55 • (all the pixels with a slope value above 55 • were initially selected as release cells). After the identification of the potential release cell locations, all the candidate pixel positions were compared to aerial pictures and erroneous release points, such as vegetated areas, were discounted. Figure 13A shows the release cells identified using this approach. Rock density was set to 2500 kg/m 3   Prior to running the simulations, the locations of release points (cells) were selected by setting a slope threshold, which was set to 55° (all the pixels with a slope value above 55° were initially selected as release cells). After the identification of the potential release cell locations, all the candidate pixel positions were compared to aerial pictures and erroneous release points, such as vegetated areas, were discounted. Figure 13A shows the release cells identified using this approach. Rock density was set to 2500 kg/m 3 Figure 13B shows the results from the first simulation of rockfall analysis 1, carried out using block volumes varying from 0.008 to 0.125 m 3 (representative of DFN 1). It is possible to see that the areas of block deposition (light green) correspond well with the debris fan and slope deposit areas, which are usually made of blocks with volumes lower than 0.1 m 3 .  Figure 13B shows the results from the first simulation of rockfall analysis 1, carried out using block volumes varying from 0.008 to 0.125 m 3 (representative of DFN 1). It is possible to see that the areas of block deposition (light green) correspond well with the debris fan and slope deposit areas, which are usually made of blocks with volumes lower than 0.1 m 3 . Figures 14A-D and 15A-D show the results of the other eight simulations of rockfall analysis 1, undertaken using the range of volumes derived from the low fracture intensity DFN 2 (from 0.6 to 500 m 3 ). The resultant locations of deposited blocks are highlighted in light green. These can be compared with the locations of rock blocks representing recent rockfall events, shown by red dots and grouped using the same range of block volume used during each simulation.
The simulations carried out with block volumes ranging from 0.6 to 1.5 m 3 ( Figure 14A) and 1.6 to 3 m 3 ( Figure 14B) show that the blocks tend to deposit in the area between the landslide scar and the Frattura Road, which links the old and historical Frattura Vecchia village to the Frattura village. When increasing the block volumes up to 7.5 m 3 ( Figure 14C) and 20 m 3 ( Figure 14D), some of the simulations suggest there is a possibility that the blocks can reach and cross the road, sometimes reaching as far as Frattura graveyard. The majority of deposited blocks, however, are still contained within the area between the slope and the Frattura Road. This changes when the volume of blocks is further increased. This can be seen in the simulations presented in Figure 15A-D, where the blocks with volumes higher than 20 m 3 often overcome the Frattura road. With volumes greater than 55 m 3 , most of the simulated rockfalls extend beyond the road and reach the Frattura graveyard.   The simulations carried out with block volumes ranging from 0.6 to 1.5 m 3 ( Figure 14A) and 1.6 to 3 m 3 ( Figure 14B) show that the blocks tend to deposit in the area between the landslide scar and the Frattura Road, which links the old and historical Frattura Vecchia village to the Frattura village. When increasing the block volumes up to 7.5 m 3 ( Figure 14C) and 20 m 3 ( Figure 14D), some of the simulations suggest there is a possibility that the blocks can reach and cross the road, sometimes reaching as far as Frattura graveyard. The majority of deposited blocks, however, are still contained within the area between the slope and the Frattura Road. This changes when the volume of blocks is further increased. This can be seen in the simulations presented in Figure 15A-D, where the blocks with volumes higher than 20 m 3 often overcome the Frattura road. With volumes greater than 55 m 3 , most of the simulated rockfalls extend beyond the road and reach the Frattura graveyard.
The results of these simulations, in terms of locations of deposited blocks, show a very good correspondence with the maps of recent rockfall events, demonstrating the validity of the rockfall models (Figures 13-15). However, it must be stressed that these simulations are based on the use of uniform ranges of block volumes, which are not representative of the in situ rock block distribution, where the relative frequency of block volume may vary. From the DFN analysis, it is possible to see that the range of block volume fluctuates from 0.008 to 1.1 m 3 in DFN 1 and from 0.6 to 500 m 3 in The results of these simulations, in terms of locations of deposited blocks, show a very good correspondence with the maps of recent rockfall events, demonstrating the validity of the rockfall models (Figures 13-15). However, it must be stressed that these simulations are based on the use of uniform ranges of block volumes, which are not representative of the in situ rock block distribution, where the relative frequency of block volume may vary. From the DFN analysis, it is possible to see that the range of block volume fluctuates from 0.008 to 1.1 m 3 in DFN 1 and from 0.6 to 500 m 3 in DFN 2 (the blocks reaching 1000 m 3 were not considered due to the very low relative frequency). The majority of block volumes varies between 0.008 and 0.16 m 3 in DFN 1 and between 0.6 and 25 m 3 in DFN 2.
To further strengthen the rockfall model, these results were compared with the relative frequency of block volume extracted from the map of recent rockfall events. Furthermore, a new map representing potential rock blocks characterizing the area UAV 2 was created and the relative frequency of these blocks calculated. The results of these two analyses, which considered block volumes ranging from 0.6 to ca. 1100 m 3 (as for DFN 2), are represented in Figure 16A,B. Figure 16A illustrates the relative frequency of block volumes extracted from the map of recent rockfall events and Figure 16B the relative frequency of block volumes calculated from the map of potential blocks characterizing the area UAV 2.

DFN 2.
To further strengthen the rockfall model, these results were compared with the relative frequency of block volume extracted from the map of recent rockfall events. Furthermore, a new map representing potential rock blocks characterizing the area UAV 2 was created and the relative frequency of these blocks calculated. The results of these two analyses, which considered block volumes ranging from 0.6 to ca. 1100 m 3 (as for DFN 2), are represented in Figure 16A,B. Figure 16A illustrates the relative frequency of block volumes extracted from the map of recent rockfall events and Figure 16B the relative frequency of block volumes calculated from the map of potential blocks characterizing the area UAV 2. It is possible to see that the frequency of block volumes of both analyses is very similar to the one derived from the DFN 2, with the majority of the blocks between 0.6 to 25 m 3 , less frequent blocks between 25 and 100 m 3 (ca. 10%-15%) and sporadic blocks over 300 m 3 . Although it was confirmed that most of the blocks have volumes below 25 m 3 , the analyses also confirm the presence of (less frequent) medium to large blocks highlighted by the DFN analysis ( Figure 10) and the map It is possible to see that the frequency of block volumes of both analyses is very similar to the one derived from the DFN 2, with the majority of the blocks between 0.6 to 25 m 3 , less frequent blocks between 25 and 100 m 3 (ca. 10-15%) and sporadic blocks over 300 m 3 . Although it was confirmed that most of the blocks have volumes below 25 m 3 , the analyses also confirm the presence of (less frequent) medium to large blocks highlighted by the DFN analysis ( Figure 10) and the map of recent rockfall events ( Figure 6). A further confirmation of the presence of such blocks is confirmed by the analysis of rockfall scars. An example of this is shown in Figure 17 where it is possible to observe the 3D model of the area UAV 2 ( Figure 17A), four examples of rockfall scars ( Figure 17B-E) and a suspended block ( Figure 17F). The calculated volume of failed blocks associated with such scars is ca. 40 m 3 ( Figure 17B), ca. 200 m 3 ( Figure 17C), ca. 60 m 3 ( Figure 17D) and ca. 230 m 3 ( Figure 17E). Figure 17F highlights the volume of an overhanging block with a volume of ca. 500 m 3 .
confirmed by the analysis of rockfall scars. An example of this is shown in Figure 17 where it is possible to observe the 3D model of the area UAV 2 ( Figure 17A), four examples of rockfall scars ( Figure 17B-E) and a suspended block ( Figure 17F). The calculated volume of failed blocks associated with such scars is ca. 40 m 3 ( Figure 17B), ca. 200 m 3 ( Figure 17C), ca. 60 m 3 ( Figure 17D) and ca. 230 m 3 ( Figure 17E). Figure 17F highlights the volume of an overhanging block with a volume of ca. 500 m 3 .  Figure 17C), ca. 60 m 3 ( Figure 17D) and ca. 230 m 3 ( Figure 17E). Figure 17F highlights the volume of overhanging block with a volume of ca. 500 m 3 .
In light of this evidence, it is clear that the use of a uniform range of volume, as commonly seen in rockfall simulations, does not represent a realistic rock mass condition. Therefore, as part of rockfall analysis 2, we developed a new approach that introduces the use of relative block volume frequency gathered from DFN models. In relation to this relative frequency and the total amount of release cells selected for rockfall simulations, we calculated the number of cells to be assigned to a specific volume. In terms of spatial location, the selection of block volumes among the release cells was randomly performed and therefore several simulations were undertaken to verify multiple possible scenarios. In each simulation, the relative frequency of block volumes was fixed, but their spatial distribution within the release cells was randomly changed. Figure 18A shows an example of the release cell maps created using the proposed relative frequency volumes; the colors represent the different volumes and it is possible to see how these are not uniform but are related to the frequency derived from DFN models. Figure 18B highlights the results of one of the simulations with relative frequency volumes. Using this approach, every map highlights the rockfall trajectories associated with the entire range of volumes gathered from the DFN model. In this case, the simulation shows that the smallest blocks arrest their travel at the toe of the slope and the larger ones can potentially reach the road and other infrastructure. The location of deposited blocks and their volumes agree with that observed during field mapping of recent rockfall events (Figures 6 and 18C), which presents the rock blocks interpreted from high-resolution orthophotographs. This can be clearly seen  Figure 17C), ca. 60 m 3 ( Figure 17D) and ca. 230 m 3 ( Figure 17E). Figure 17F highlights the volume of overhanging block with a volume of ca. 500 m 3 .
In light of this evidence, it is clear that the use of a uniform range of volume, as commonly seen in rockfall simulations, does not represent a realistic rock mass condition. Therefore, as part of rockfall analysis 2, we developed a new approach that introduces the use of relative block volume frequency gathered from DFN models. In relation to this relative frequency and the total amount of release cells selected for rockfall simulations, we calculated the number of cells to be assigned to a specific volume. In terms of spatial location, the selection of block volumes among the release cells was randomly performed and therefore several simulations were undertaken to verify multiple possible scenarios. In each simulation, the relative frequency of block volumes was fixed, but their spatial distribution within the release cells was randomly changed. Figure 18A shows an example of the release cell maps created using the proposed relative frequency volumes; the colors represent the different volumes and it is possible to see how these are not uniform but are related to the frequency derived from DFN models. Figure 18B highlights the results of one of the simulations with relative frequency volumes. Using this approach, every map highlights the rockfall trajectories associated with the entire range of volumes gathered from the DFN model. In this case, the simulation shows that the smallest blocks arrest their travel at the toe of the slope and the larger ones can potentially reach the road and other infrastructure. The location of deposited blocks and their volumes agree with that observed during field mapping of recent rockfall events (Figures 6 and 18C), which presents the rock blocks interpreted from high-resolution orthophotographs. This can be clearly seen by comparing Figure 18B,C, and highlights the good agreement between the spatial distribution of deposited blocks and their volumes. The differences in the accumulation of blocks in the two maps, i.e., the significant number of blocks with volumes between 120 and 260 m 3 in the map of Figure 18B, is related to different amount of blocks present as shown in the two maps. The map presented in Figure 18B shows the simulation of  Figure 18C shows the location of ca. 600 rock blocks interpreted from the high-resolution orthophotographs and field inspection.
Remote Sens. 2020, 12,2053 22 of 28 by comparing Figure 18B,C, and highlights the good agreement between the spatial distribution of deposited blocks and their volumes. The differences in the accumulation of blocks in the two maps, i.e., the significant number of blocks with volumes between 120 and 260 m 3 in the map of Figure 18B, is related to different amount of blocks present as shown in the two maps. The map presented in Figure 18B shows the simulation of 3100 rockfalls and their end point locations, while Figure 18C shows the location of ca. 600 rock blocks interpreted from the high-resolution orthophotographs and field inspection.

Discussion
This paper describes the integrated use of UAV data and DFN modeling to improve the results of rockfall simulation. Mount Rava is used as the case study and is the site of the large Scanno paleolandslide. The failure left a large landslide scar that significantly changed the morphology of the mountain, forming steep rock slopes on its SW flank, which are often affected by local rockfall events. The slope is ca. 400 m high and 600 m wide and very difficult to access and survey. The use of UAV made it possible to reach inaccessible areas within the landslide scar and to significantly improve the data for geomechanical characterization. Imagery of the most geologically complex area (area UAV 1) was acquired from a distance of ca. 36 m, providing a model resolution of ca. 13 cm/pix. To obtain such a high-resolution model, however, the required processing time was very long and the model management complex. Therefore, the other two areas (UAV 2 and 3) were acquired from a distance of 46 m, obtaining a model resolution of ca. 17-19 cm/pix, which is a high resolution and suitable for the two areas, allowing recognition of all the geological features. In agreement with Della Seta et al. [47] and Francioni et al. [48], three main discontinuity sets are recognized. The 3D models and sampling windows extracted from UAV data make it possible to also denote important changes in the fracture intensity in the different areas of the landslide scar. The changes in the fracture intensity is likely related to the thickness of the bedding and the fracture corridors in fault damage zones in some areas of the slope under study. This agrees with previous studies carried out by Francioni et al. [48] and Miccadei et al. [49], which demonstrated the presence of a complex fault system.
Consideration of the variation in the fracture intensity is critical during the creation of DFN models and calculation of block volumes. Hence, we developed two different DFN models, DFN 1 and 2, representing high and low intensity fracture areas, respectively.
In this study the DFN model was developed by integrating conventional geomechanical and UAV-extracted data. The DFN were validated against the generated fracture intensity values, which were gathered from the analysis of sampling windows and the recently developed fracture analysis software FracPaq. This software allows the calculation of P 21 using several circles within the sampling windows. This approach makes the calculation of fracture intensity more rigorous, allowing for the calculation in every sampling window of an average P 21 and the analysis of their relative standard deviation and minimum and maximum values. The 2D section extracted from the 3D DFN was analyzed using the same approach/software, facilitating comparison and making the validation more rigorous. The software FracPaq was developed mainly for structural geological problems, as documented by Giuffrida et al. [58] and Watkins et al. [59]. In this paper, we adopted the software to improve the analysis of fracture intensity, P 21 . The use of this approach/software is relatively new in engineering geology and it may represent a tool to further improve the development of DFN models for the analysis of slopes.
The range of block volumes obtained from DFN 1 and DFN 2 were validated against the map of recent rockfall events and the maps of potential rock blocks in the area UAV 2 ( Figure 16A,B) and subsequently used as input in 3D rockfall simulation using the software Rockyfor3D. This software was adopted in several rockfall analyses in recent years. For example, Vanneschi et al. [34] combined the use of remote sensing techniques and Rockyfor3D to perform a study of rockfall runout and geological hazard in a natural slope in Italy. Corona et al. [60] investigated the uncertainties related to the choice of parameters accounting for energy dissipation and surface roughness. Moos et al. [61,62] carried out two different analyses to analyze the effect of specific tree species and vegetation on rockfall risk. In these examples, as in many others documented in literature, the authors performed several simulations changing the range of rock volumes released. This allows verification of the response of the rockfall model under different scenarios. In Rockyfor3D, the block dimensions defined in each source cell can be varied randomly within a predefined percent (based on a defined uniform volume variation between ±0% and ±50%) before each simulation. However, this variation has to be uniform and can reach a maximum of 50% of the chosen volume. Therefore, the use of this approach for wide ranges of block volumes is not possible and several simulations have to be carried out while gradually increasing the range of volumes. Furthermore, in cases where the distribution of block volumes is not uniform, the results may not be representative of the real rock mass condition.
In the present research, a new approach is introduced which is based on the use of relative frequency of block volume extracted from DFN models. The difference between this approach and the methods presented in literature is that once the rockfall model is calibrated and validated, every rockfall simulation will include a realistic DFN-based distribution of block volumes, without the need to perform multiple simulations changing the range of block volume. This allows presentation of the results in the form of a single final simulation, where the location of deposited blocks and their respective volume can be shown. Furthermore, the resulting map of kinetic energy will represent the entire range of potential energy dissipated upon impacts, avoiding the need to create different kinetic energy maps in relation to the change of block volumes. An example of a kinetic energy map extracted from the simulation results shown in Figure 18B is provided in Figure 19. This approach allows the analysis to include a more realistic distribution of rock volumes, where the potential block volumes are randomly distributed among the release cells. However, it is still necessary to perform multiple simulations varying the distribution of block volumes among the release cells to verify all possible scenarios. The choice of the best scenario can be established by comparing the results of the simulation with a map of recent rockfalls; where such as map is unavailable, then all the scenarios have to be considered.
Remote Sens. 2020, 12, 2053 24 of 28 out while gradually increasing the range of volumes. Furthermore, in cases where the distribution of block volumes is not uniform, the results may not be representative of the real rock mass condition. In the present research, a new approach is introduced which is based on the use of relative frequency of block volume extracted from DFN models. The difference between this approach and the methods presented in literature is that once the rockfall model is calibrated and validated, every rockfall simulation will include a realistic DFN-based distribution of block volumes, without the need to perform multiple simulations changing the range of block volume. This allows presentation of the results in the form of a single final simulation, where the location of deposited blocks and their respective volume can be shown. Furthermore, the resulting map of kinetic energy will represent the entire range of potential energy dissipated upon impacts, avoiding the need to create different kinetic energy maps in relation to the change of block volumes. An example of a kinetic energy map extracted from the simulation results shown in Figure 18B is provided in Figure 19. This approach allows the analysis to include a more realistic distribution of rock volumes, where the potential block volumes are randomly distributed among the release cells. However, it is still necessary to perform multiple simulations varying the distribution of block volumes among the release cells to verify all possible scenarios. The choice of the best scenario can be established by comparing the results of the simulation with a map of recent rockfalls; where such as map is unavailable, then all the scenarios have to be considered.

Conclusions
This paper presents the combined use of UAV and DFN models to provide improved results of rockfall simulation. It was demonstrated that the use of UAV imaging is very important in the case of wide and inaccessible slopes, and especially in the case of very complex geological and structural settings. In this study, we highlighted, due to the presence of fault zones, how fracture intensity varies within the Scanno landslide scar. The calculation of fracture intensity was carried out by using the recently developed software FracPaq. Although this software was developed for structural geological purposes, it was shown that the approach based on the calculation of several P21 values in Figure 19. Map of the kinetic energy extracted from the simulation shown in Figure 18B.

Conclusions
This paper presents the combined use of UAV and DFN models to provide improved results of rockfall simulation. It was demonstrated that the use of UAV imaging is very important in the case of wide and inaccessible slopes, and especially in the case of very complex geological and structural settings. In this study, we highlighted, due to the presence of fault zones, how fracture intensity varies within the Scanno landslide scar. The calculation of fracture intensity was carried out by using the recently developed software FracPaq. Although this software was developed for structural geological purposes, it was shown that the approach based on the calculation of several P 21 values in every sampling window can be useful for improved slope stability studies and generation of representative DFN models. In the case of DFN 1, the resultant block volumes were very small, with majority of blocks below 0.1 m 3 . When analyzing the DFN 2, the range of block volume changed dramatically, with 70% of block volumes ranging between 0.5 and 25 m 3 . Approximately 20% of the rockfall block volumes were between 25 and 100 m 3 , ca. 5% were between 100 and 200 m 3 and ca. 2% between 200 and 300 m 3 . Some blocks generated had volumes up to 500 m 3 (0.6%) and 1000 m 3 (0.3%). The presence of blocks of such dimensions was confirmed by map analysis of recent rockfall events and the map of potential rock blocks in area UAV 2 ( Figure 16). Considering this, and the wide range of block volumes gathered from DFN analysis, we performed two suites of rockfall analysis simulations. During rockfall analysis 1, different rockfall simulations were undertaken, where there was a gradual increase in block dimensions using a conventional uniform block volume distribution. Using this procedure, we calibrated and validated the rockfall model against the map of the end point locations of recent rockfall events. The good match between this map and the results of rockfall simulations in terms of the size of the deposited blocks demonstrated the validity of the proposed models.
To further improve the quality of the rockfall simulation results, a new approach that includes the relative frequency of block volumes derived from the DFN for every single simulation was developed in the rockfall analysis 2. The results obtained from these analyses showed that it is possible to integrate into rockfall simulations a more realistic relative frequency distribution of block volumes using the results of DFN analyses. The results of the models were also compared with a map of recent rockfall events. The output gathered from this method allowed the visualization of the location of deposited blocks (classified according to volumes ranges) and their potential kinetic energy on a single map.