SfM-MVS Photogrammetry for Rockfall Analysis and Hazard Assessment Along the Ancient Roman Via Flaminia Road at the Furlo Gorge (Italy)

: Rockfall events represent signiﬁcant hazards for areas characterized by high and steep slopes and therefore e ﬀ ective mitigation controls are essential to control their e ﬀ ect. There are a lot of examples all over the world of anthropic areas at risk because of their proximity to a rock slope. A rockfall runout analysis is a typical 3D problem, but for many years, because of the lack of speciﬁc software, powerful computers, and economic reasons, a 2D approach was normally adopted. However, in recent years the use of 3D software has become quite widespread and di ﬀ erent runout working approaches have been developed. The contribution and potential use of photogrammetry in this context is undoubtedly great. This paper describes the application of a 3D hybrid working approach, which considers the integrated use of traditional geological methods, Terrestrial Laser Scanning, and drone based Digital Photogrammetry. Such approach was undertaken in order to perform the study of rockfall runout and geological hazard in a natural slope in Italy in correspondence of an archaeological area. Results show the rockfall hazard in the study area and highlights the importance of using photogrammetry for the correct and complete geometrical reconstruction of slope, joints, and block geometries, which is essential for the analysis and design of proper remediation measures.


Introduction
Rockfall is a common natural process acting on rock slopes. Such phenomenon is usually characterized by unpredictability and high velocities, representing a significant hazard especially near communication routes, villages, workplaces (for instance either open pit mines or dams), and other popular areas (such as archaeological, historical, or religious sites). A rockfall event can be very dangerous, causing casualties even when a relatively small mass is mobilized [1,2], and can be triggered by one or a combination of different factors, such as climate, seismic activity, and vegetation growth. These factors can cause the opening of joints, pore pressure increase, freeze-thaw processes, or chemical weathering of rock [3]. Once a rockfall has occurred, the trajectory and the energy of falling blocks are controlled by a complex interaction of multiple factors including the shape and size of the block, micro-and macro-geometry of the slope, and its surface properties (e.g., roughness, material properties, land cover). Hence, it is important to study the possible trajectory and energy of a rockfall hazard varies both spatially and temporally, statistical based approaches are often preferred [11]. In this paper, a hybrid approach that combines both deterministic and statistical methods was adopted to perform 3D runout analysis of a natural rock slope in Italy ( Figure 1). The site, where some rockfall events have already occurred in the past, was analyzed because of the risk posed by the presence of the local ancient Roman Via Flaminia road, close to an important archaeological site of Roman age, and the Church named "Santa Maria delle Grazie." A specific multi-source survey, using a Total Station (TS), GNSS, TLS, and UAV, was performed to characterize the steep rock slope (practically inaccessible without mountaineering-climbing equipment and relative capacity) and determine the location and geometry of potentially unstable blocks visually identified by photointerpretation. The collected data was then utilized to perform specific slope stability studies, 3D runout analyses, and to produce a new rockfall hazard map. Specifically, the aim of this study was to make available to professionals of Enel Green Power O&M Hydro Italy with hydro-civil function of the Furlo dam, a deterministic analysis containing rock block volumes, weights, stability SFs, and runout analyses to be used as quantitative data for protection measures planning. often preferred [11]. In this paper, a hybrid approach that combines both deterministic and statistical methods was adopted to perform 3D runout analysis of a natural rock slope in Italy ( Figure 1). The site, where some rockfall events have already occurred in the past, was analyzed because of the risk posed by the presence of the local ancient Roman Via Flaminia road, close to an important archaeological site of Roman age, and the Church named "Santa Maria delle Grazie." A specific multisource survey, using a Total Station (TS), GNSS, TLS, and UAV, was performed to characterize the steep rock slope (practically inaccessible without mountaineering-climbing equipment and relative capacity) and determine the location and geometry of potentially unstable blocks visually identified by photointerpretation. The collected data was then utilized to perform specific slope stability studies, 3D runout analyses, and to produce a new rockfall hazard map. Specifically, the aim of this study was to make available to professionals of Enel Green Power O&M Hydro Italy with hydro-civil function of the Furlo dam, a deterministic analysis containing rock block volumes, weights, stability SFs, and runout analyses to be used as quantitative data for protection measures planning.

Study Area
The study area is in the Municipality of Fermignano, Province of Pesaro-Urbino. The Furlo Gorge has developed between the Pietralata Mount (889 m a.s.l.) and the Paganuccio Mount (976 m a.s.l.) due to the erosional power of the Candigliano River; during the years the river has reached a considerable depth in this area, which is currently no longer visible due to the presence of the dam, built in 1922, which reduced the impetuous river to a placid lake ( Figure 1).
In Roman times, with the aim of connecting people in such an impervious area, the Emperor Vespasian dug a tunnel in the narrowest point of the gorge that was called "petra pertusa" or "forulum" (small hole), from which the name "Furlo" was generated. Next to the Furlo, there is an earlier Etruscan passageway, 8 m long, 3.30 m wide, and 4.45 m high, and a small Church, also called "della Botte," once inhabited by a hermit. The tunnel was completely carved into the rock by a chisel, in which graded cuts are still visible.
The slopes are made up of limestone of marine origin belonging to the Mesozoic Umbria-Marche Succession of the Italian Apennine [59]. Limestones mainly include the Calcare Massiccio Formation dated to Early Jurassic (i.e., Hettangian-Early Sinemurian stages). This formation is characterized by the absence of stratification with large banks typical of carbonate platform environment but presents numerous high-angle fault systems and discontinuities, both tectonic and arising from unloading [60][61][62]. The color of the Calcare Massiccio Formation is typically white or greyish, while its thickness reaches up to 300 m [59].

Study Area
The study area is in the Municipality of Fermignano, Province of Pesaro-Urbino. The Furlo Gorge has developed between the Pietralata Mount (889 m a.s.l.) and the Paganuccio Mount (976 m a.s.l.) due to the erosional power of the Candigliano River; during the years the river has reached a considerable depth in this area, which is currently no longer visible due to the presence of the dam, built in 1922, which reduced the impetuous river to a placid lake ( Figure 1).
In Roman times, with the aim of connecting people in such an impervious area, the Emperor Vespasian dug a tunnel in the narrowest point of the gorge that was called "petra pertusa" or "forulum" (small hole), from which the name "Furlo" was generated. Next to the Furlo, there is an earlier Etruscan passageway, 8 m long, 3.30 m wide, and 4.45 m high, and a small Church, also called "della Botte," once inhabited by a hermit. The tunnel was completely carved into the rock by a chisel, in which graded cuts are still visible.
The slopes are made up of limestone of marine origin belonging to the Mesozoic Umbria-Marche Succession of the Italian Apennine [59]. Limestones mainly include the Calcare Massiccio Formation dated to Early Jurassic (i.e., Hettangian-Early Sinemurian stages). This formation is characterized by the absence of stratification with large banks typical of carbonate platform environment but presents numerous high-angle fault systems and discontinuities, both tectonic and arising from unloading [60][61][62]. The color of the Calcare Massiccio Formation is typically white or greyish, while its thickness reaches up to 300 m [59].
From a structural-geological perspective, the study area is located within a wide anticline of Apennine origin which has at the core the Calcare Massiccio Formation; the compressive stresses responsible of the Apennine chain formation acting since at late Paleogene period (i.e., late Oligocene epoch) generated a NW-SE oriented anticline with a sub-vertical axial plane, sub-horizontal axis, and an open cylindrical geometry [61].
The core shows a collapse structure identified by two direct faults NW-SE oriented that lower the central portion of the fold and determine the formation of an elongated graben along the direction of the axial plane. The fault displacements are of a few tens of meters. The anticline is an uprooted structure, translated to NE above a basal thrust which is temporally extended towards NE, from the bottom upwards, involving more recent stratigraphic intervals in the deformation [63]. The anticline shows a radius of curvature greater than 2 km and is characterized by a symmetry with orientation N 125-130 and a barely noticeable gradient towards NE.
In the Furlo Gorge collapses are extremely widespread with several sub-vertical fronts of various origins. Slope detachments occur very often due to physical-chemical weathering along stratification beds, networks of tectonic and/or newly formed joints produced by unloading processes [64].
Previous research on rockfalls at the Furlo Gorge report of past events that often affected the Via Flaminia Road, sometimes resulting in its closure while restoration occurs [63]. The publications also include engineering-geological data (e.g. rock unit weight, rock mass classification, joint friction angle) and rockfall runout tests which have been used to integrate and validate the results of this research.

Topographic, TLS, and UAV Surveys
In order to recreate a detailed morphology of the slope, useful for the stability analysis and the rockfall runout modelling, an integrated TLS-UAV survey was carried out, allowing acquisition of a detailed 3D point cloud of the area. The integrated approach was fundamental given the complex morphology of the slope. This resulted in the limited use of TLS, restricting its deployment to the lower lateral parts of the rock mass, because of (i) the presence of water, which impedes data acquisition of the frontal parts of the slope, and (ii) the slope steepness, which impedes data acquisition in the higher parts of the slope. Therefore, the use of a UAV photogrammetric survey was essential to complete the data acquisition of the whole slope. On the other hand, the presence of vegetation and other obstacles precludes the possibility of a safe flight on the lateral lower parts of the rock mass. Therefore, the integration of UAV survey and TLS was essential to provide a complete model of the higher and the bottom zones of the slope and the interior of the tunnel. UAV photogrammetry was used extensively along the slope to map the inaccessible areas and provide overlapping data which increased the density and aided in the registration of the datasets.
TLS from 11 different positions was carried out by using a Trimble™ X8 device. After every scan, the scanner was alternated with a Nikon™ D7100 digital camera, properly mounted on a bracket (i.e., Leica™ Nodal Ninja 3II), to acquire high-definition (HD) images. Photos were shot using a fisheye lens (8 mm focal length); in this way, having a wide field of view (i.e., 180 • ), only six photos (every 60 • ) were necessary to guarantee an overlap of 66% and cover the whole panorama (the dip angle was set to −10 • ). In the laboratory, photos were processed using specific software for the creation of panoramic images 360 • view (PTGui™-New House Internet Services BV). PTGui has the capability to model fisheye projections and to find reference points in the images. Having the focal center of the camera lens in correspondence with the optical center of the scanner, the HD photos were aligned (i.e., texture mapping) to the 3D point clouds of every scan position. Following this, a 3D photorealistic point cloud suitable for photointerpretation was created from the TLS survey. The different scans were relatively aligned using Iterative Closest Point (ICP) algorithm [65], while the absolute alignment was carried out by measuring the coordinates of nine pre-placed mobile optical targets located in strategic positions, visible from different TLS scans. All the georeferencing and texturing processes were carried out using Trimble™ RealWorks 10.4. The adopted reference system was the Italian Gauss-Boaga, datum Monte Mario Italy 2. The coordinates of optical targets were measured by the combined use of TS (Leica™ Nova MS50) and two GNSS receivers (Leica™ GS15). The two GNSS devices firstly worked independently, in static modality, for a period longer than 3 h; then, the two measured points were processed using Leica™ Geo Office software and differential methods by combining simultaneous records from six other permanent GNSS stations of the Leica SmartNet ItalPos national network (i.e., R. San Marino/RSMARINO, Pesaro/PES2, Pennabili/PEN2, Gubbio/GUB2, Città di Castello/CIT1, Arcevia/ARCEVIA). The orthometric height of the measured points was calculated by using ConverGo software [66]. This procedure allowed a sub-centimetric accuracy for the two GNSS points, later necessary for georeferencing (i.e., roto-translation) the optical targets measured during the TS survey. In order to reach different parts of the slope, the TS was moved several times from its initial position by using the intersection method; this operation was carried out every time with an accuracy of about 1 mm.
The flights were performed using a DJI Phantom 4 PRO drone, equipped with a 20-megapixel digital camera with 1-inch sensor. The UAV surveys were carried out with direction of photo acquisition in nadiral modality (perpendicular to the lake) and in frontal modality (perpendicular and slightly oblique to the rock faces). This was done because SfM technology is based on sophisticated algorithms of image matching that use pseudo-random redundant images acquired from multiple viewpoints to reconstruct the three-dimensional geometry of an object or surface; hence, multiple images obtained from different angles help the image alignment procedure and limit non-linear deformations. Four flights were executed to cover all sectors of the slope, with UAV's distance from outcrops varied from 10 m up to 80 m, giving a nominal overlap and sidelap of at least 80% and 60%, respectively. An average estimated distance between pixel centers measured on the ground (i.e., Ground Sample Distance-GSD) of 1.8 cm was calculated.
The complex morphology of the area, with a deep and narrow gorge, did not allow for an automatic flight with waypoints: it was flown manually with an automatic triggering of photos every 2 s. Moreover, the use of a UAV equipped with a Real Time Kinematic (RTK) GNSS, which could have allowed acquisition of accurate coordinates of the camera location at every shot, was also not possible: the limited satellite visibility and the high Positional Dilution Of Precision (PDOP) within the gorge could have determined very low spatial accuracy in data. Therefore, the exterior orientation of all images (i.e., 1232 photos) was done by using 80 Ground Control Points (GCPs), which distribution is shown in Figure 2, measured again by the combining TS measurements, GNSS surveys in static modality, and differential data post-processing with a final millimeter-level positioning accuracy.
Eight artificial targets were pre-placed in the accessible zones of the study area, while the rest were obtained using natural features visible in the photos throughout the slope and by the TS during fieldwork. Their location was decided considering a balance between an optimum spatial distribution, both in space and elevation from the road, and easy identification of points on the images. The software Agisoft™ MetaShape (version 1.5.3; http://www.agisoft.com, last access: 24 June 2019) was used to process the images obtained from the UAV surveys. This software can solve the camera interior and exterior orientation parameters and generate georeferenced spatial data such as 3D point clouds, DEMs, DTMs and orthophotos. Agisoft™ MetaShape based on (i) SfM technology for the exterior orientation of images and the creation of the sparse 3D point cloud and (ii) MVS photogrammetry for the dense cloud reconstruction [67,68].
The first processing step consisted of image alignment, through which the interior and relative orientation parameters were solved. In order to improve the whole alignment process and to obtain low re-projection error, millions of tie points were automatically extracted without setting a point limit. Following image alignment, the second processing step involved georeferencing of the 3D model in such a way as to solve the exterior orientation parameters by using the GCPs coordinates measured during the GNSS-TS topographic surveys. The exterior orientation of images was necessary to measure the orientation, with respect to the North, and inclination of slopes and discontinuities which are needed for the stability analysis.
RGB data for texturing. The different point clouds were imported and no further alignment was executed. This is due to the fact the engineering-geological analyses were carried out on blocks and discontinuities only visible on the UAV-point cloud. The TLS-point cloud, instead, was fundamental for the complete characterization of the slope geometry, which is essential for the 3D rockfall runout analysis. As already mentioned, this was mainly due to the presence of trees on the lower part of the slope, which covered a portion of the rock mass in proximity to the local road. This caused distortions and holes on the UAV-point cloud, with possible consequences on the block trajectories. Therefore, the UAV-point cloud was replaced by the TLS-point cloud where needed (lower lateral parts of the slope). In this view, slight offsets between the two point clouds were considered negligible because, as explained later in Section 4.3.2, the spatial resolution of the input raster maps for the runout analysis was set equal to 50 cm/pixel, more than the possible error of mis-alignment of the two different point clouds.  Subsequently, the "optimize" tool of Agisoft™ MetaShape was utilized to adjust the estimated camera positions and to remove possible non-linear deformations, minimizing the errors due to re-projection and misalignment of the photos. Moreover, the optimization improved the model by deleting all the tie points with a re-projection error greater than one pixel.
In a subsequent step, the dense 3D point cloud was generated with medium quality and aggressive depth filtering settings. Automatic classification of cloud was done trying to remove vegetation, as much as possible, from the area of interest.
Finally, the different point clouds obtained by TLS and UAV image processing, both georeferenced in the Italian National Gauss-Boaga reference system, were unified using the open source CloudCompare™ software to produce a final complete 3D model of the rock mass including RGB data for texturing. The different point clouds were imported and no further alignment was executed. This is due to the fact the engineering-geological analyses were carried out on blocks and discontinuities only visible on the UAV-point cloud. The TLS-point cloud, instead, was fundamental for the complete characterization of the slope geometry, which is essential for the 3D rockfall runout analysis. As already mentioned, this was mainly due to the presence of trees on the lower part of the slope, which covered a portion of the rock mass in proximity to the local road. This caused distortions and holes on the UAV-point cloud, with possible consequences on the block trajectories. Therefore, the UAV-point cloud was replaced by the TLS-point cloud where needed (lower lateral parts of the slope). In this view, slight offsets between the two point clouds were considered negligible because, as explained later in Section 4.3.2, the spatial resolution of the input raster maps for the runout analysis was set equal to 50 cm/pixel, more than the possible error of mis-alignment of the two different point clouds.

Engineering-Geological Characterization of the Rock Mass and Stability Analysis
A direct engineering-geological survey was possible only on the lower parts of the slope. Unfortunately, the accessible areas allowed for the analysis of only few discontinuities, considered not fully representative of the conditions of the higher parts. For this reason, visual inspection in the field provided only a qualitative assessment of the rock mass conditions and the measurement of few joint characteristics (e.g., joint aperture, infill, weathering, compressive strength). Given the impossibility to collect all the necessary discontinuity properties, a proper characterization of the rock mass using traditional methods such as, for example, the Rock Mass Rating (RMR, [69,70]) was not possible. In the place of it, the Hoek and Brown Geological Strength Index (GSI, [71]) was determined since it bases on qualitative geological observations.
The structural analysis of the higher parts of the slope was carried out by using the Compass plugin [72] of the freeware CloudCompare™ software. Compass is a plugin which allows measurement of discontinuity orientations directly on the point cloud. In this way, the orientation of about 70 discontinuities in strategic areas of the slope was acquired to create a stereonet used to identify the different discontinuity sets. Following this, a kinematic stability analysis, based on the Markland test [73] was performed. Since the results of the kinematic analysis are influenced by geotechnical properties and topographic factors (i.e., dip and dip direction of rock outcrops), the slope was divided in four sectors ( Figure 3) with different orientations.

Engineering-Geological Characterization of the Rock Mass and Stability Analysis
A direct engineering-geological survey was possible only on the lower parts of the slope. Unfortunately, the accessible areas allowed for the analysis of only few discontinuities, considered not fully representative of the conditions of the higher parts. For this reason, visual inspection in the field provided only a qualitative assessment of the rock mass conditions and the measurement of few joint characteristics (e.g., joint aperture, infill, weathering, compressive strength). Given the impossibility to collect all the necessary discontinuity properties, a proper characterization of the rock mass using traditional methods such as, for example, the Rock Mass Rating (RMR, [69,70]) was not possible. In the place of it, the Hoek and Brown Geological Strength Index (GSI, [71]) was determined since it bases on qualitative geological observations.
The structural analysis of the higher parts of the slope was carried out by using the Compass plugin [72] of the freeware CloudCompare™ software. Compass is a plugin which allows measurement of discontinuity orientations directly on the point cloud. In this way, the orientation of about 70 discontinuities in strategic areas of the slope was acquired to create a stereonet used to identify the different discontinuity sets. Following this, a kinematic stability analysis, based on the Markland test [73] was performed. Since the results of the kinematic analysis are influenced by geotechnical properties and topographic factors (i.e., dip and dip direction of rock outcrops), the slope was divided in four sectors ( Figure 3) with different orientations.
Such analysis allowed identification of the main types of failure (i.e., toppling, planar sliding, wedge sliding). In addition, CloudCompare was used to estimate volumes, geometries, and positions of all potentially unstable blocks that could trigger rockfall events in the future. This provided a basis to identify source blocks to be used in the subsequent dynamic stability analysis and rockfall runout simulation. The geometric information can be extracted from the point cloud manually, or utilizing software to automate the process [74,75]. For example, in [76] the authors used Polyworks (InnovMetric Software, 2017) to define the sets of the discontinuities based on 3D point clouds and then used the spacing of each set to roughly calculate the average size of blocks. In their work, the Such analysis allowed identification of the main types of failure (i.e., toppling, planar sliding, wedge sliding).
In addition, CloudCompare was used to estimate volumes, geometries, and positions of all potentially unstable blocks that could trigger rockfall events in the future. This provided a basis to identify source blocks to be used in the subsequent dynamic stability analysis and rockfall runout simulation. The geometric information can be extracted from the point cloud manually, or utilizing software to automate the process [74,75]. For example, in [76] the authors used Polyworks (InnovMetric Software, 2017) to define the sets of the discontinuities based on 3D point clouds and then used the spacing of each set to roughly calculate the average size of blocks. In their work, the discontinuity information was extracted manually from the point clouds and the block size was calculated approximately. Differently, in [77] a fully automatic method for extracting rock block information from 3D point clouds was presented.
In this paper blocks were identified on the final point cloud by a careful photointerpretation that was also revised and validated by experts. The following procedure was used: (i) the point cloud was segmented in several areas containing a potentially unstable block each; (ii) every joint located below a potentially sliding block was sampled and co-planar points selected; (iii) the cloud representing each block was cut and the normals to the points, oriented toward the exterior, computed; (iv) the clouds representing both the block and the underlying joint were merged and used to build a mesh through the Poisson Surface Reconstruction Plugin [78]; (v) finally, by means of the editing tools of CloudCompare, the volume of every mesh was measured.
The dynamic stability analysis of rock blocks was carried out by Rocscience™ codes (i.e., Rocplane, Swedge, Roctopple for planar sliding, edge sliding and block toppling, respectively). SFs were calculated both in static and dynamic conditions considering the presence of water, and the local seismic acceleration. Using GeoStru™ PS software (GeoStru, 2019) the maximum horizontal acceleration "ag" of the site was derived (i.e., 0.25 g in this study). The procedure was carried out following the new Italian "Norme Tecniche per le Costruzioni" [79].

Rockfall Runout Analysis
The mathematical model that can be adopted to analyze the rockfall runout must be able to describe the behavior of a falling block in terms of trajectory direction and length, bounce height, speed, energy, and type of movement. Two principal methods can be used to solve this problem: (i) lumped mass methods and (ii) rigorous methods [80]. The lumped mass approach considers the single block as a point with specific mass and velocity. Therefore, it is a simplified model where the velocity of the block is reduced after every impact with the surface by using two coefficients Kn (normal coefficient of restitution) and Kt (tangential coefficient of restitution). Differently, a more rigorous approach considers the shape and dimension of the falling block, which must be known a priori. It is a complete method that considers the moment changes, depending on shape of block and surface, angle of trajectory, roughness of the surface, and friction between block and surface [80].
The method used in this paper follows a 3D rigorous approach. The simulation software is called Rockyfor3D [81] and combines physically based deterministic algorithms with stochastic approaches. In this sense, Rockyfor3D is also defined as "probabilistic process-based rockfall trajectory model." The trajectory is simulated as a 3D vector data by sequences of parabolic free fall through the air and rebounds on the surface. Sliding of the rock is not considered, as well as block fragmentation, and rolling of the mass is represented by a sequence of short-distance rebounds [81]. Moreover, the model considers impacts against trees. The simulations can be run by using 14 ASCII rasters characterized by the same cell size and extent, which describes the characteristics of slope surface and falling blocks. The input raster maps, created in a GIS environment using the data from the geomatic and geological surveys, are described in Table 1.
The runout simulations with the current slope configuration were performed using the input blocks identified on the UAV images and measured on the 3D point cloud. However, it has to be mentioned that, before proceeding with the rockfall runout simulation, a model calibration step was done; this process, also defined as "back analysis" in respect of a past rockfall event, which is essential to calibrate the entire 3D model together with the simulation input parameters and raster. In this case, given the steep geometry of the slope, which is characterized essentially by steep bedrock with limited vegetation overhanging water dump, it is difficult to simulate a past rockfall event, since the arrival point is generally missing (i.e., below water level). In terms of rockfall runout analysis, well-represented geometrical info were obtained: not only slopes, but also the potentially unstable blocks were well known in terms of shape and size. Physical parameters were calibrated (i.e., coefficients of normal and tangential restitution) in such a way that input data can be considered adequate for a reliable rockfall runout study. Table 1. Input raster maps necessary for using Rockyfor3D.

DTM
It represents the DTM of the slope.
Rock density It represents the rock density of every simulated source cell (i.e., falling block).
Block height, width and length These three raster maps describe the shape of a block, in terms of height, width, and length.
Block shape It defines the shape of a falling block (i.e., rectangular, ellipsoidal, spherical, and disc shaped).

Roughness 70, 20, 10
These three raster maps define the roughness of the slope surface, represented by all the obstacles lying on the slope. The value of every cell, expressed in meters, corresponds to the height of a representative obstacle that can be encountered by a block with a probability of 70% (rg70), 20% (rg20) and 10% (rg10). The values range from 0 to 100 and are used to calculate the tangential coefficient of restitution (Rt) causing an energy loss.

Soil type
This raster map defines the elasticity of the ground surface basing on the soil type; there are eight different soil types (e.g., bedrock, asphalt road, fine soil material) and every type corresponds to a certain value of the normal coefficient of restitution (Rn).

D Modelling
The final integrated point cloud (TLS + UAV) is composed of 22,308,471 points, corresponding to an average point spacing of 1.91 cm. The 3D model was built in mesh format constituting of 18,753,142 facets. Even if the final obtained 3D model is considered suitable for the goals of the study, it must be mentioned that the two point clouds come from different processing workflows (TLS and SfM-MVS photogrammetry), using different control points used for georeferencing and, therefore, have diverse spatial accuracy. The point cloud from TLS has a Registration Error of 5.45 mm (i.e., mis-alignment among the 11 scans) and a Georeferencing Error of 10 mm in respect to the artificial targets. The point cloud from SfM-MVS photogrammetry has a Root Mean Square Error (RMSE) of 80 mm in respect of GCPs. It is assumed that this higher error is due to the lack of targets properly fixed along the slope and the uncertainty in defining them on the UAV images.

Engineering-Geological Characterization
As already mentioned in Section 3.2, the traditional engineering-geological survey was limited due to inaccessibility of the study areas. Therefore, a combined approach was used with the interpretation of the point cloud from remote sensing data. The rock mass conditions, estimated through the GSI by field inspections and interpretation of photographs from UAV, show a difference between the bottom and the higher parts of the slopes. In this regard, GSI is considered to be varying from 70-80 (lower parts) to 50-60 (mid-higher parts). In fact, it was possible to observe how the rock slope is influenced (especially in the mid and higher parts) by several fractures related to different joint systems. A total of 67 discontinuities were measured, in terms of dip and dip direction, directly on the final point cloud through the CloudCompare™ Compass plugin. Then, the presence of joint sets was determined by means of a density analysis of attitudes in stereographic projection: the lower hemisphere of the Schmidt (equal area) net was used (Figure 4). Table 2 summarizes the orientation values of each joint set.
Results from the engineering-geological analysis, revealed a certain variability of the discontinuity pattern, with a difference between the basal and the mid-higher parts of the rock slope. In addition, the discontinuities dip angles increase, sometimes reversing their direction (an example is shown in Figure 5). Figure 5 shows the identification of the main discontinuities on the mid parts of the slope, the area most subjected to unstable block formation.   Results from the engineering-geological analysis, revealed a certain variability of the discontinuity pattern, with a difference between the basal and the mid-higher parts of the rock slope. In addition, the discontinuities dip angles increase, sometimes reversing their direction (an example is shown in Figure 5). Figure 5 shows the identification of the main discontinuities on the mid parts of the slope, the area most subjected to unstable block formation.    Results from the engineering-geological analysis, revealed a certain variability of the discontinuity pattern, with a difference between the basal and the mid-higher parts of the rock slope. In addition, the discontinuities dip angles increase, sometimes reversing their direction (an example is shown in Figure 5). Figure 5 shows the identification of the main discontinuities on the mid parts of the slope, the area most subjected to unstable block formation.

Kinematic Stability Analysis
After identification of the main discontinuity sets, a kinematic stability analysis was carried out to identify the possible failure conditions on the slope. This type of analysis is based on discontinuity orientations, slope direction and discontinuity friction angle. As previously explained, the discontinuity orientation and slope direction were derived from the point cloud, while the friction angle of the failure surfaces (25 • ) was chosen based on previous studies on the same area [63], the authors expertise, and conservative geotechnical assumptions. Table 3 summarizes the results of the kinematic analysis for the identified four main slope directions.
Results shown in Table 3 indicate that along V1 the formation of planar sliding (on J1) and wedge sliding (on the intersection between J1-J4) are possible (example in Figure 6a). With the same failure modes present in V2, with the addition of possible direct toppling (example in Figure 6b). Differently, on V3 only planar sliding on J1 is possible. Along V4 wedge sliding due to the intersection between J3 and J4 is possible. After identification of the main discontinuity sets, a kinematic stability analysis was carried out to identify the possible failure conditions on the slope. This type of analysis is based on discontinuity orientations, slope direction and discontinuity friction angle. As previously explained, the discontinuity orientation and slope direction were derived from the point cloud, while the friction angle of the failure surfaces (25°) was chosen based on previous studies on the same area [63], the authors expertise, and conservative geotechnical assumptions. Table 3 summarizes the results of the kinematic analysis for the identified four main slope directions.

Unstable Block Identification and Dynamic Stability Analysis
The kinematic analysis provided an understanding of possible block formations useful for the geological hazard assessment of the area. Attention was placed on the dynamic stability conditions and the rockfall runout analysis of seven potentially unstable blocks/rock masses (A, B, C, D, F, E, K in Figure 7) that were identified on the textured point cloud. Their geometrical characteristics are summarized in Table 4 together with the most critical SF as calculated in dynamic conditions considering the presence of water and seismic forces. Results shown in Table 3 indicate that along V1 the formation of planar sliding (on J1) and wedge sliding (on the intersection between J1-J4) are possible (example in Figure 6a). With the same failure modes present in V2, with the addition of possible direct toppling (example in Figure 6b). Differently, on V3 only planar sliding on J1 is possible. Along V4 wedge sliding due to the intersection between J3 and J4 is possible.

Unstable Block Identification and Dynamic Stability Analysis
The kinematic analysis provided an understanding of possible block formations useful for the geological hazard assessment of the area. Attention was placed on the dynamic stability conditions and the rockfall runout analysis of seven potentially unstable blocks/rock masses (A, B, C, D, F, E, K in Figure 7) that were identified on the textured point cloud. Their geometrical characteristics are summarized in Table 4 together with the most critical SF as calculated in dynamic conditions considering the presence of water and seismic forces.

Rockfall Runout and Hazard Assessment
As already mentioned in the Methods section, before proceeding with the rockfall runout simulation, a back analysis in respect of a past rockfall event was executed to properly calibrate the 3D model, together with the simulation input parameters and raster. This analysis was done considering reliable trajectories in which end-points were located along the Via Flaminia road, in the archaeological area, and within the lake. After this, the simulations were performed by considering each identified block, assuming real block dimensions measured directly on the point cloud. The final simulations were configured as follows: (i) number of trajectories for each block equal to 50; (ii) block volume weight 2500 kg/m 3 ; (iii) random variation of the block volume equal to ±10%. The spatial resolution of the input raster maps was set equal to 50 cm/pixel. Experiences reported in [82] showed that the preferred resolution lies between 2 × 2 m and 10 × 10 m. In addition, a 1 × 1 m resolution does not necessarily improve the quality and it increases the amount of data substantially. In respect to this, the cell size adopted in this work is already considered too detailed.
Some representative outputs of the rockfall runout simulation are shown in Figure 8. The figure shows the results in terms of:

•
Kinetic energy at the 95% confident level (translational + rotational, in kJ) of all blocks that passed over a single cell.

•
Passing height values at the 95% confident level (measured from the block barycenter in normal direction to the slope surface, in meters).
In order to evaluate the rockfall geological hazard for the area, the values of probability, kinetic energy, and passing heights of each block run in the simulations were re-classified as shown in Table 5. These values were obtained based on expertise and knowledge of characteristics of the net barriers already present in the area. The values of energy, for example, were chosen considering the potential to stop blocks within a range between 1100 kJ (as minimum) and 5000 kJ (as maximum), representative of the capacity of a standard barrier available on the market.
By the combination of the re-classified raster maps, three different classes of hazard were established (i.e., low, medium, high). A preliminary hazard class was assigned to each cell based on the combination of Height and Energy (Figure 9). Subsequently, this preliminary class was combined with the Probability value in order to define the final rockfall hazard class. For example, a block having high energy but low rebound height and low probability of reaching a certain point of the slope would be assigned to the medium hazard class. Differently, a block with medium rebound height and energy and a high reach probability, would be assigned to the high hazard class. The final hazard map obtained from the described approach is shown in Figure 10. By the combination of the re-classified raster maps, three different classes of hazard were established (i.e., low, medium, high). A preliminary hazard class was assigned to each cell based on the combination of Height and Energy ( Figure 9). Subsequently, this preliminary class was combined with the Probability value in order to define the final rockfall hazard class. For example, a block having high energy but low rebound height and low probability of reaching a certain point of the slope would be assigned to the medium hazard class. Differently, a block with medium rebound height and energy and a high reach probability, would be assigned to the high hazard class. The final hazard map obtained from the described approach is shown in Figure 10.   By the combination of the re-classified raster maps, three different classes of hazard were established (i.e., low, medium, high). A preliminary hazard class was assigned to each cell based on the combination of Height and Energy ( Figure 9). Subsequently, this preliminary class was combined with the Probability value in order to define the final rockfall hazard class. For example, a block having high energy but low rebound height and low probability of reaching a certain point of the slope would be assigned to the medium hazard class. Differently, a block with medium rebound height and energy and a high reach probability, would be assigned to the high hazard class. The final hazard map obtained from the described approach is shown in Figure 10.

Simulation with Virtual Net Counters
With the aim of proposing rockfall control defenses for the ancient Roman Flaminia road and the church, an additional simulation with virtual net counters was run. A greater number of simulations were performed (i.e., 1000 trajectories/source point) with the aim of increasing the reliability of the analysis [83]. Four net counters were included in the simulations as shown in Figure 11. The nets were virtually placed at the bottom of the rock slopes, collecting detailed data relative to (i) number of overpassing blocks, (ii) kinetic energy (95th percentile of the simulated energy), (iii) passing heights (95th percentile of the simulated passing heights), and (iv) reach probability for each net. passing heights (95th percentile of the simulated passing heights), and (iv) reach probability for each net.
The analysis indicated that it is theoretically possible to define the most suitable position where to localize rockfall barriers or other remediation measures. Results of this simulation are summarized in Table 6.    The analysis indicated that it is theoretically possible to define the most suitable position where to localize rockfall barriers or other remediation measures. Results of this simulation are summarized in Table 6.

Discussion
The study area is a steep slope characterized by the presence of four different discontinuity sets, that isolate the portion of rock mass forming blocks prone to collapse. Considering that a large part of the slope is inaccessible, the use of geomatic techniques was essential in this study for geological mapping purposes. Indeed, it was possible to remotely collect data about blocks and discontinuities to be used for a detailed geometric characterization of the slope in a spatially uniform way throughout the area. However, it should be stated that without a classical engineering-geological survey, fundamental properties of discontinuities and rock would be missing (e.g., joint aperture, infill, weathering, roughness). This is an aspect where future research could concentrate in developing new technologies and approaches, but the contribution of traditional surveys is today still essential for a correct characterization of a rock mass.
The integration of data from UAV photogrammetry and TLS was applied in this study with the aim of producing a detailed and complete 3D model of the slope, used for collecting accurate geometrical information that is needed for the stability analysis and rockfall runout simulation. To achieve this, the TLS was performed since it is a rapid technique that ensures at the same time high spatial resolution and accuracy of measured data. UAV photogrammetry was utilized because it allows the interpretation of the higher parts of the slope which are inaccessible and not visible from the Via Flaminia road. Moreover, the use of UAV was essential for the study purposes in consideration of the narrow shape of the gorge, the total inaccessibility of surrounding slopes and the lack of an open space from which a complete survey from higher position was achieved. In addition, photogrammetry provided a complete inspection and photointerpretation in a cost-effective manner, in just few hours of fieldwork. The final point cloud, representative of the whole slope including the ancient Roman Via Flaminia tunnel and the "Santa Maria delle Grazie" Church with a centimetric level of accuracy, allowed acquisition of deterministic discontinuity information (dip and dip direction). This helped in identifying the discontinuity sets (Figure 4), from which the kinematic stability analyses were derived. In particular, the slope is characterized by a minor presence of discontinuities in the lower parts, while most discontinuities are present in the mid-higher parts, where they show also higher dip slope, with few examples of local overturning ( Figure 5).
Interpretation of results from the kinematic analysis (Table 3) revealed also that the most probable failure mode is by planar sliding on J1 and wedge sliding along the intersection between J1 and J4 sets.
The final point cloud was also used to identify unstable blocks and measure their geometry. The dynamic stability analysis revealed that rock blocks A, D, E, and K are at a limit of equilibrium (SF 1.0), while the area of block B, that overheads the Church, needs cleaning and stripping of small stone elements that are present. It must be noted that the present study refers to current conditions which are unpredictably changeable (e.g., water, ice, chemical-physical alteration, interactions between different blocks). At present, the blocks and wedges identified in the area of interest clearly show the presence of rock bridges that favor the stability.
The runout analysis of the seven blocks identified thanks to the geomatic data, indicate the possibility that falling blocks could reach the local road, the archaeological site and the adjacent Church. In detail, results show low probability for the local road to be reached by a block (<4.5% for block E and <8% for block B taking as reference the virtual net counter number 4). On the contrary, the adjacent Church shows a higher theoretical probability (i.e., 85% for block E and 55% for block B taking as reference the virtual net counter number 3). A significant risk is associated to the reference virtual net counters 1 and 2 with possible rockfall of blocks A, C, D, F, and K (probability generally higher than 80% or 90%). However, these values do not consider the return period parameter and the probability that a single block could fall. Nevertheless, this can be considered as a conservative approach, given that if those parameters had been considered in the calculation, the probability would have certainly decreased. Moreover, it must be mentioned that the used code for rockfall runout simulation (i.e., Rockyfor3D) does not consider block fragmentation that could influence falling element trajectory, velocity, energy, bounce height, and end point. This is in agreement with the conservative approach of the study, and it could be included in future analyses.
The information collected during this study will be used by local authorities for developing projects of remediation works. Regarding this, the virtual net counters provide knowledge of the energy of every block as well as the passing height and velocity, which are fundamental information when planning rockfall protection measures. Results of the rockfall simulation showed how some of the falling blocks can produce high energies, high rebounds, or both. For example, the analysis showed how blocks E and D may reach the local road or the adjacent Church with a high energy, superior to the capability of a standard net barrier available on the market. On the other hand, blocks A and C showed lower energy, but high bounce heights that could pass an eventual barrier. Consequently, a possible immediate measure could be the removal or the consolidation by means of a harness with ropes, or with panels of cables laterally anchored on the stable rock, of blocks A, C, D, and E. Blocks smaller than 3 m 3 could be destroyed or removed.
More generally, the following actions could be carried out: • inspection, surface cleaning, and controlled barring for blocks smaller than 1 m 3 ; • destruction and removal of the smallest blocks; • wire mesh with reinforced ropes; • wire mesh panels; • rock bolts, harness with ropes, or wire rope panels; • sealing/consolidation with suspensions and consolidating mixtures; • reinforcing wall construction; Results from this study must be used to focus further analyses on the segments with higher probability to be reached by a block, for a proper design of protecting net barriers. This is very important in view of the fact that the hazard map obtained from the analysis ( Figure 10) shows a medium to high hazard for the entire area. This is critical as the area is a popular place given the presence of the archaeological site, the Church, and a local ancient Roman road.

Conclusions
This case study has shown an application of a combined working approach that integrates the use of traditional engineering-geological techniques and geomatic technologies with the aim of carrying out a 3D rockfall runout analysis of a natural slope in Italy. The study confirms how integration of UAV photogrammetry and TLS data can be useful to overcome the common problems found in traditional approaches, such as the lack of data of adequate scale factor, the site accessibility, the unsafe conditions, and the complete identification of possible rockfall source points. TLS or UAV photogrammetry, alone, would not have produced a complete model of this area: only their integration, carried out respecting the basic guidelines of data processing and georeferencing, could guarantee the full and reliable modelling of the slopes.
Several discontinuities and rock blocks were identified above the archaeological area and near the Church of "Santa Maria delle Grazie" along the ancient Roman Via Flaminia road. The complete 3D model of the slope provided a means to characterize the rock mass discontinuity patterns and define the possible falling trajectories of each block. At present, each wedge is characterized by theoretical limit equilibrium conditions (with possible presence of unknown rock bridges), so that the adoption of some remediation measures (e.g., barring, chemical desegregation, removal, net barriers installation) is pointed out. In the opinion of the authors, the main advantages of the adopted approach for rockfall runout analysis, even in consideration of the study aims, can be summarized as it follows: • Provides the possibility to obtain deterministic geometrical data about rock blocks (i.e., position on the slope, shape and volume) and discontinuities (i.e., dip and dip direction) in inaccessible areas with short working times and in a cost-effective manner thanks to SfM-MVS photogrammetry and UAV surveys.

•
Provides the possibility to assess the stability conditions of single blocks, to compute their SFs, and to identify possible rockfall sources, thanks to the acquisition of detailed information about rock wedges volume and weight, slope orientation, and rock mass conditions. • Provides the possibility to perform geological hazard studies that identify areas with high probability to be reached by rockfall events.

•
Provides the possibility to simulate the presence of net counters that allow collection of information regarding energy, passing heights, and reach probability of blocks that can be used as quantitative data for protection measures planning.
However, it must be said that some limitations are also present in such a study. One main issue with SfM-MVS photogrammetry is quantifying the uncertainty and accuracy without the use of a number of GCPs. Especially in zones where the narrow shape of the gorge makes impossible a proper use of UAV equipped with an RTK-GNSS, the number of GCPs should be higher than usual (such as this case study). This, combined with the height of slopes, makes the work more difficult despite an extensive and extreme use of TS. Laser scanning from a UAV represents an important tool and an opportunity in case studies as the one presented in this paper since it can acquire 3D point clouds with constant spacing and high density. Nevertheless, the use of images, and particularly the adoption of SfM-MVS photogrammetry, gives the advantage of making possible the photointerpretation and the deep study of features.
The presented case study highlights the contribution of geomatics in the interpretation and analysis of engineering-geological investigations and demonstrates a way to overcome practical issues of traditional survey techniques related to the need of accurate and detailed data and the difficult access to rock slopes.