Identification of potential rockfall sources using UAV-derived point cloud

Recent advances in remote sensing techniques and computer algorithms allow accurate, abundant, and high-resolution geometric information retrieval for rock mass characterization from 3D point clouds. The automatic application of the extracted information for local scale rockfall susceptibility assessment, where discontinuities characteristics play a major role in rocky slope stability, requires step by step logical procedures. This paper presents a novel methodology to use the extracted discontinuity set characteristics for a local scale rockfall susceptibility assessment, tailored for Uncrewed Aerial Vehicle (UAV) data acquisition. The method consists of 4 steps: (i) 3D slope model reconstruction using UAV digital photogrammetry, (ii) automatic characterization of discontinuity sets, (iii) slope stability analysis, and (iv) susceptibility assessment using a new Rockfall Susceptibility Index. The proposed method was applied to a road cut rocky slope in a mountainous area of the Samaria National Park, in Crete Island, Greece. Visual validation indicates that the areas of higher and moderate rockfall susceptibility on the 3D model of the rocky slope are adjacent to rockfall source areas marked by the presence of fallen blocks on the foot of the slope. The proposed methodological workflow presents novelties related to the use of point clouds for the estimation of the Rock Quality Designation (RQD) index, the visualization of discontinuity set spacing, the evaluation of the persistence and the Slope Mass Rating (SMR) index, as well as the incorporation of the persistence of overhangs into the rockfall susceptibility assessment and visualization.


Introduction
In the last decades, the rocky slope mapping and stability has transitioned from two dimensional (Baillifard et al. 2003;Irigaray et al. 2003;Günther et al. 2012;Yilmaz et al. 2012) to three dimensional (Abellán et al. 2006(Abellán et al. , 2014(Abellán et al. , 2016Gigli et al. 2014;Bonilla-Sierra et al. 2015;Assali et al. 2016;Menegoni et al. 2019;Zhang et al. 2019b;Buyer et al. 2020;Papathanassiou et al. 2020). Three-dimensional models can provide critical information for rockfall source detection and susceptibility identification along the slope height, and with a high spatial resolution (Dunham et al. 2017;Matasci et al. 2018). Given the primary importance of source identification for the rockfall hazard assessment, and of the effect of block volume for the rockfall propagation (trajectory, energy, rebound height), methodologies based on the use of 3D model data need to be developed towards this end.
Moving from regional to local scale, rockfall susceptibility assessment methods variate from the GIS spatial analysis tools (Saroglou 2019) that exploit data layers including lithology and presence of geological faults, slope angle, and intensity of triggering mechanisms (earthquake and rainfall intensity maps) towards those based on the analysis of the local geological structure and strength of the rock mass (Andrea et al. 2010;Irigaray et al. 2003;Yilmaz et al. 2012). Local rockfall susceptibility assessment commonly takes place through the application of indicators such as the Slope Mass Rating index (SMR) of Romana (1993), which considers the angular relationship between the slope and discontinuities alongside basic rock mass characteristics such as discontinuity spacing, persistence, roughness, aperture, and infilling. These strategies can lead to the detection of potential source areas as points or lines to be used for the rockfall run-out simulation. Using such points or lines as potential rockfall initiation areas, although conservative, thus on the safe side, is not often accurate since sources areas that are eventually stable and will not generate block detachment are considered (Matasci et al. 2018). For more reliable rockfall simulation results, accurate locations of potential source areas are required, which implies better characterization of the rock mass at a local scale (Gigli et al. 2014).
The use of modern surveying techniques such as Light Detection and Ranging (LiDAR) and Uncrewed Aerial Vehicle (UAV)-based imagery and its processing with digital photogrammetry brought a new potential for rockfall sources and instability identification. These techniques allow the generation of high-resolution 3D models of the slope surface in the format of point clouds or 3D meshes, where accurate geometric information from rock masses can be extracted using manual procedures or automatic routines via computer algorithms, among them open-source software (Guo et al. 2019;Riquelme et al. 2014aRiquelme et al. , 2015Riquelme et al. , 2018Slob 2010;Sturzenegger et al. 2011;Ünlüsoy and Süzen 2020;Zhang 2019a).
Although many publications deal with the feature extraction of the rock mass characteristics from point clouds, few works, such as Matasci et al. (2018) and Dunham et al. (2017), have developed a methodology on how to apply these features for automatic rockfall source identification. Papathanassiou et al. (2020) provided a comprehensive methodology to characterize a 3D model (point cloud) of a rock mass by means of the Rock Mass Rating and Slope Mass Rating indices; nevertheless, the identification of potential rockfall source areas were made only by visual assessment without incorporating the extracted features of the 3D model in an integrated workflow.
Moreover, attempts to identify the extent of potentially detachable rock masses directly on a 3D point cloud environment are scarce (Bonneau et al. 2019;Chen et al. 2017;Farmakis et al. 2020;Menegoni et al. 2020). For instance, recent studies on kinematic stability analysis consider the angular relationship between the rock mass discontinuities and the slope in stereogram plotting (Alameda-Hernández et al. 2019;Menegoni et al. 2019;Wang et al. 2019), but the extent of the areas to be affected on the point cloud cannot be directly deducted from that. Such information is important for an effective rockfall hazard analysis, where usually the source areas are not refined across the slope, and could therefore compromise the propagation results.
Recent contributions for the rockfall susceptibility assessment on LiDAR-derived point clouds have been made by Dunham et al. (2017) and Matasci et al. (2018). Dunham et al. (2017) developed the Rockfall Activity Index (RAI), a slope morphology-based index, allowing the centimeter accurate point cloud to be classified into increasing rockfall activity levels based on a simple logic tree algorithm. Matasci et al. (2018) developed a routine to quantify the rockfall susceptibility at a local to a sub-regional scale, covering up to thousands of square meters. The final output is a classified point cloud, providing the locations of susceptible areas, considering planar, wedge, and toppling failure mechanisms and overhanging parts of the rock cliff. The routine to quantify rockfall susceptibility was later applied in an augmented reality environment by Zhang et al. (2019b). Even though the works of Dunham et al. (2017) and Matasci et al. (2018) consider the presence of overhangs for the susceptibility analysis, the extension of the overhangs was not taken into account. The extension of the overhangs is an indication of how large the unstable areas are and therefore constitute relevant information for rockfall hazard.
The objective of this work is to present a methodology that uses UAV imagery-derived point clouds to identify and visualize potential rockfall sources on rocky walls. The proposed methodology consists of four steps: (i) 3D point cloud generation, (ii) discontinuity set characterization on the point cloud, (iii) slope stability analysis by means of the SMR index, and (iv) calculation of the Rockfall Susceptibility Index on the 3D model of the rocky slope, as indicative of rockfall stability. The Rockfall Susceptibility Index incorporates as parameters the overhangs, the SMR index of Romana (1993), the discontinuity set persistence, and spacing.
Different from the previous works, the proposed methodology consists of using an index for rockfall susceptibility assessment, which is based on data that can be extracted from already available discontinuity characterization algorithms for point clouds. However as these algorithms rarely incorporate this information into a unified methodology to visually indicate potential rockfall sources on the slope model, this work focuses on that. A procedure to visualize the spacing, persistence, and the SMR index values directly on the point cloud is shown, instead of the commonly used non-visual (numerical) data. As a part of the presented methodology, we also propose the assessment of the Rock Quality Designation (RQD), according to Deere and Deere et al. (1989), exclusively by data obtained by the point cloud model of the slope surface. Finally, for the visualization of the rockfall sources, we suggest the incorporation of the extension of overhangs, considering that larger overhang areas are indicative of higher hazard.
As the aim here was to develop a methodology which is mostly based on data extracted from the 3D model of the slope, the innovative use of techniques for the calculation or visualization of the involved parameters is also described and explained in detail in the following sections.
The proposed methodology was applied to a road cut rocky slope in the wider mountainous area of the Samaria National Park, in Crete Island, Greece. The validation of the susceptibility results is based on the visual observation of past events, which indicates that the areas calculated as being of higher and moderate rockfall susceptibility on the slope are adjacent to past rockfall sources, that have led to relatively more and larger deposited blocks on the foot of the slope. The approach herein proposed is suitable for rocky slopes along transportation corridors such as road cuts and railways, open-pit mines, and rock cliffs at a local scale (up to hundreds of square meters).

Case study
To develop and test the proposed methodology, we worked on a road cut situated in the broader context of Samaria National Park, in Crete, Greece. The latter is a park located in western Crete, on the southern slope of the White Mountains (Fig. 1), a World's Biosphere Reserve as established by UNESCO in 1981. Tectonic activity, erosion, and karstification shaped the mountainous landscape of the area (Vogiatzakis and Rackham 2008), creating plenty of narrow, tall, and vertical rocky slopes, with Samaria Gorge being the steepest, tallest, and narrowest opening (Spanos et al. 2008).
The island of Crete is part of the Hellenic arc, a lithospheric structure originated by an orogenic process due to the convergence of the African and Eurasian tectonic plates (Manutsoglu et al. 2003). This lithospheric framework favors seismicity in the island, with earthquake intensity ranging from 0.12 to 0.24 of acceleration coefficient g (Saroglou 2019). In the past century, the closest epicenters to Crete had a magnitude of up to 6.3 according to the earthquake catalogue of Makropoulos et al. (2012). Prevalent in the geological setting are the high pressurelow temperature metamorphic rocks belonging to the Plattenkalk unit overlain by the Phyllite-Quartzite unit (Seidel 2003). At a regional scale, the Plattenkalk unit incorporates the Mavri formation (Lower Liassic) in the lower part of the stratigraphic column up to the Aloides formation (Eocene) in the upper part (Manutsoglu et al. 2003), which consists mainly of cherty calcite marbles, dolomite, phyllite, and a calcareous metaflysh (Seidel 2003).
The rainfall season is mainly concentrated between September and April, with a mean annual precipitation of 1244 mm estimated from 7-year data recorded by the meteorological station in the Samaria Gorge (National Observatory Athens 2020).

Field data collection
The studied rocky slope is located on a road cut approximately 1620 m above mean sea level, between the entrance of the Gorge and the Kallergis mountain refugee (Fig. 1), a place to host hikers and climbers. The rocky slope is approximately 14 m in height, 46 m in lateral extension, and has an average slope angle of 51° towards the road. It consists of dark grey platy limestone with bedding planes of decametric thickness, and with occasional centimetric intercalations of quartz or calcite.
Four principal discontinuity sets were identified (Fig. 2) and their orientation was obtained manually using a geological compass, except DS 4 (Table 4) since it was not visually identified during the field investigations as being significantly different from other discontinuity sets. The compass measurements were obtained along 2 m of a total 4 m long scanline in the upper part of the rocky slope as well as on specifically chosen surfaces since the scanline did not intersect all the discontinuity sets.
The discontinuity sets are almost perpendicular to each other, forming parallelepiped shaped rock blocks. The most representative discontinuity set is the bedding plane (DS 1) with the same dip as the slope surface (51°). Small rockfalls constitute the main process of mass wasting, where the rock blocks bounded by DS 2 laterally and by DS 3 and DS 4 in the front and back, slide on the bedding plane surface (DS 1) in a planar failure mechanism followed by a free fall. The block fragments measured at the bottom of the slope have a volume ranging from 0.064 m 3 to 0.001 m 3 , or edge dimensions from 40 to 5 cm. Vegetation is abundant above the slope and some small bushes are present on the slope surface too. Chemical weathering is not prominent and no seepage was observed during the fieldwork.
The risk of rockfalls in this area is related to vehicle transiting the road connecting the village of Omalos, next to the entrance of Samaria Gorge and the Kallergis Refugee. Although this road is not highly transited, this slope was selected as an example for the application of the method and the processing of the point cloud, to be further extended to cut and natural slopes where risk is higher, and thus susceptibility and hazard analysis are essential.

Methodology
The proposed methodological workflow consists of four steps ( Fig. 3): (i) 3D slope model reconstruction using UAV digital photogrammetry, (ii) automatic characterization of discontinuity sets, (iii) slope stability analysis, and (iv) rockfall susceptibility assessment.

The 3D slope model reconstruction using UAV digital photogrammetry
The red/green/blue (RGB) images taken with a UAV platform are used as input for the reconstruction of the 3D slope model as a point cloud employing Structure from Motion (SfM), a photogrammetry technique that uses computer algorithms to extract key points in overlapping images taken from multiple view angles to create 3D models (Westoby et al. 2012). Once the point cloud is generated, points that do not represent the surface of the slope, such as those corresponding to vegetation, soil, and rock fragments in the surroundings, should be removed to reduce the size of the file and to avoid misclassification of the automatic detection of the discontinuity surface, therefore improving the results of the classification (Menegoni et al. 2019;Riquelme et al. 2018).

Automatic characterization of discontinuity sets
The discontinuity set characterization performed onto the 3D point cloud is performed using the open-source software Discontinuity Set Extractor (DSE) developed in MAT-LAB by Riquelme et al. (2014a), and consists of the semiautomatic identification and definition of the discontinuity sets (DS). In particular, for each DS it is possible to obtain the dip, the dip direction, the persistence, and the normal spacing.
The discontinuity set dip and dip direction are obtained semi-automatically in the DSE software in three main steps following the methodology proposed by Riquelme et al. (2014a): local curvature calculation, statistical analysis of the planes, and cluster analysis. In this approach, each cluster is described mathematically as a plane equation and visually as a group of 3D points with the same local normal vector orientation. In practical terms, each cluster represents an exposed surface on the rocky slope corresponding to a The scanline in a has a total of 4-m extension from A to A′ discontinuity and it is colorized with the same color as the other clusters assigned to the same discontinuity set. The output is a.txt file of the classified 3D point cloud, containing for each point, the following attributes: 3D coordinates (X, Y, Z), Discontinuity Set id (Ds id), Cluster id (Cl id), and Cluster plane equation parameters (A, B, C, D). The main limitation of this approach is that only the 3D points belonging to exposed surfaces of the rocky slope are identified as discontinuity surfaces.
An algorithm developed by Riquelme et al. (2015) and implemented in the DSE software allows the automatic normal spacing computation from the results of the discontinuity set dip and dip direction (the classified 3D point cloud). The parameter D for each cluster represents its position in space, and thus, the normal spacing is computed as the orthogonal distance between the position of neighboring clusters. The user can choose to compute the normal spacing considering two scenarios: full persistent or non-persistent discontinuities (see Riquelme et al. 2015 for full description). The output is a.txt file containing: the Discontinuity Set id (Ds id), the Pairs of neighboring clusters id (Cl id), and the Spacing between clusters as the orthogonal distance (m). This computational approach requires that the normal vector of each cluster plane is equal to the normal vector of the discontinuity  (Riquelme et al. 2015). In practical terms, it implies that the studied slope in which this method is applied should have parallel discontinuities on a given set.
The computation of the discontinuity persistence is automatically performed also from the classified 3D point cloud using the DSE software, according to the methodology developed by Riquelme et al. (2018). If the discontinuities in the studied slope are fully persistent, the user can set a k threshold to merge coplanar clusters and form a new merged cluster for the computation of persistence. If the discontinuities are non-persistent, k can be set as 0 so that each cluster is considered as an individual discontinuity for the computation. For each cluster (if k = 0) or merged cluster (if k > 0), the persistence is computed in three directions (dip, dip direction, and length of the maximum chord) and in terms of the area of the cluster. These persistence values are not visualized in the 3D point cloud but are provided numerically in a.txt file as an output from DSE software, similarly to the spacing calculation, which includes the following data: Discontinuity Set id (Ds id), Clusters id (Cl id), persistence along the dip (m), persistence along the dip direction (m), persistence along the maximum length (m), and the persistence in terms of area (m 2 ). As the output files from the discontinuity set normal spacing and persistence computation do not have the 3D position of the cluster, the visualization of these rock mass characteristics is not plotted in the point cloud environment. On step (iv) of the methodology, we propose a workflow chart (Fig. 4) for their visualization directly on the 3D model as well as on the use of these characteristics as indicators for the rockfall susceptibility assessment.

Slope stability analysis
The rocky slope stability is analyzed by means of the SMR index for each discontinuity set, using the open-access calculator SMRTool developed by Riquelme et al. (2014b). This procedure requires first the characterization of basic geomechanical parameters of the rock mass for the Rock Mass Rating-RMR b of Bieniawski (1989). For this case study, no drill cores were available for the calculation of the RQD index, one of the geomechanical parameters for the RMR b . An alternative to calculate the RQD index directly on the 3D point cloud is proposed here, which follows a similar procedure as described by Deere and Deere (1989) for drill cores. Our approach consists of selecting sections on the surface of the slope which are perpendicular to the prevalent discontinuity sets, so that the rock blocks formed by the intersection of trace lengths in the slope surface larger than 100 mm are summed and afterward divided by the total length of the section. The obtained RQD index was afterward validated by the estimated RQD index from the correlation to the volumetric discontinuity count J V for cubic shape rock fragments (Eq. 1), after Palmstrom (2005): where the J V is calculated as the inverse sum of (n) discontinuity sets having each a mean discontinuity set spacing S n as shown in Eq. 2:

Rockfall susceptibility assessment
A novel approach for the integration of geomechanical properties of the rock mass on the 3D model to provide a visualization of the rockfall susceptibility assessment in 3D is described in Fig. 4. First, the geomechanical properties are evaluated as susceptibility indicators and visualized. Then the Rockfall Susceptibility Index is calculated, integrating these indicators. Using the proposed procedure, potential rockfall sources and their susceptibility can be visualized directly on the point cloud.
The following susceptibility indicators are used: -The presence of overhangs (I o ), indicating a lack of support for the overlying rock mass. First, the classified 3D point cloud is overlaid with the RGB point cloud to evaluate which cluster (Cl id) represents overhanging surfaces belonging to the corresponding discontinuity set. For those clusters, the value of 1 is given to all its 3D points as a new attribute on the classified 3D point cloud whereas the value of 0 is given to all the other clusters. -The minimum spacing between discontinuities (I S ), as an indicator of the rock mass quality. Discontinuities with smaller spacing translate to a poorer quality rock and thus a higher probability for rock detachment. In computing the discontinuity spacing, a given cluster can have more than one normal spacing value if it has more than one neighboring cluster (Fig. 5). To account for a conservative scenario, the value of the smallest spacing is given to all the 3D points belonging to that given cluster as a new attribute on the classified 3D point cloud. -The discontinuity persistence in terms of maximum length (I P ), indicating the extent of lack of support for the overlying rock mass. Focus on the persistence of the DS 1, DS 2, DS 3, and DS 4 is also given here, as it indicates areas of weak lateral support of the rock mass and areas being more susceptible to rockfalls. To account for a conservative scenario, the length of the maximum chord computed with the DSE software was used. The degree in which this assumption is conservative or realistic has to be verified by observation of the study site in situ or virtually on the point cloud. The value of persistence of each cluster was given to all the 3D points belonging to that cluster of a particular discontinuity set, as a new attribute in the.txt file of the classified 3D point cloud. -The SMR index (I SMR ), indicating instabilities mainly based on kinematic criteria, without directly considering the shear resistance of the rock mass. In this approach, the computed SMR class of each discontinuity set is assigned to all its clusters and corresponding 3D points, also as a new attribute in the classified 3D point cloud.
Considering the aforementioned indicators, each point comes with the following attributes, with their respective range of values, as explained: -3D coordinates (X, Y, Z) The point attributes can be visualized at the 3D space by setting them as scalar fields using a point cloud visualization software. Here, CloudCompare (2020) v.2.10.2. is used. For each attribute (I o I S I P I SMR ) of a point, a score of 0, 1, or 2 is given. The scores are defined based on the influence of the attribute on the quality of the rock mass (I S ), the stability of the slope (I SMR ), and the lack of support for overlaying rock masses (I o , I P ). For instance, lower values of normal spacing indicate lower rock mass quality according to the RMR index (Bieniawski 1989), and therefore receive a higher value in the scoring for susceptibility. Discontinuities with SMR classes considered unstable according to Romana (1993) were also given a high value of scoring. The presence of overhangs scores higher too.
The scores are empirically established as shown in Table 1. For the I S and I P , the assessment was carried out by visualizing the numerical value on the 3D space and checking the range of values of I S and I P in areas of low (scoring 0), moderate (1), or high susceptibility (2). For the I SMR , the classes which are proposed by Romana (1993) were used: class I described as no failures was considered as of low susceptibility (scoring 0), class II and III described with some failure of blocks as of moderate (1), and the remaining classes IV and V with bigger failures as of high susceptibility (2). For the I o , overhang surfaces (I o = 1) were characterized as highly susceptible (scoring 2), whereas the rock Fig. 5 Schematic representation of normal spacing calculation (S) for non-persistent discontinuity of a given discontinuity set. In this example, 4 clusters (Cl) are separated by different positions in space D. Cl 2 has a smaller normal spacing S 1-2 computed to Cl 1 compared to the spacing S 2-3 to cluster Cl 3 . A similar scenario applies for Cl 3 , where the computed spacing S 2-3 to Cl 2 is smaller than the spacing S 3-4 to cluster Cl 4 . In our approach, all the 3D points of Cl 2 receive the value of S 1-2 as a new attribute whereas all the 3D points of Cl 3 receive the value of S 2-3  (0). The use of these thresholds was further verified and validated, after application to the case study.
The scored indicators are summed up with the same weight to give the Rockfall Susceptibility Index I RF (Eq. 3) ranging from 0 to 8: The classification of I RS into low, moderate, and high susceptibility levels depends on the specific site characteristics of a given slope and thus requires each time adjustment to the local conditions. In this case study, the I RS was classified into three susceptibility classes: low for scores up to 3, moderate from 3 to 5, and high from 5 to 8, assuming an equivalent distribution of the scores into the classes. A higher index in a particular area in the slope does not guarantee that a rock detachment will occur in the future, but it is rather an indication of areas in the slope more prone than others based on the considered indicators (overhangs, spacing, persistence, and SMR class). The highlighted as highly (3) I RF = I S + I P + I SMR + I O and moderate susceptible areas on the 3D point cloud can be afterward used for the next steps of rockfall hazard assessment including detailed block volume calculation, intensity, and trajectory of fallen blocks.

Application to a rocky slope in Samaria National Park, Crete
The 3D slope model reconstruction using UAV digital photogrammetry The UAV data acquisition of RBG images was carried out using the inbuilt 12.4-megapixel resolution camera on board of a DJI Phantom 4 quadcopter. A manual flight mode was performed in which the UAV stops in the air to capture the image and then moves to the next spot to acquire the next image, instead of taking images while moving. This procedure reduces the chance of having motion-blurred images. Variable oblique views of the rocky slope were preferred since it provides a better point cloud resolution, for the given geological structure, and avoids occlusion. This was particularly important in this case study to capture the overhanging areas of the slope, as they are required for the rockfall susceptibility assessment. The distance from the UAV to the  Table 2 Input parameters used for discontinuity set extraction in the DSE software Parameter Value K nearest neighbor-k nn 20 points Tolerance for the coplanarity test-ɳ max 0.3 Number of bins for the density analysis-nbin 64 User-defined number of discontinuity sets-DS 4 Minimum angle between principal poles-γ 1 20°M aximum angle between a pole and its principal pole (cone filter)-γ 2 30°M inimum number of points per cluster-ppc 100 points Cluster distribution threshold for cluster alignment-k 0 surface of the slope was less than 15 m to obtain a sufficient level of detail of the geological structures. A total of 236 RGB images were acquired and after a visual quality assessment, 221 were selected for image processing, discarding the low contrast and/or highly blurred ones. The photogrammetric software Pix4Dmapper version 4.3.31 was used to reconstruct the 3D model of the slope as a dense point cloud, using the standard processing parameters for the 3D model template. The generated point cloud contained ~ 11 million points, with ground sampling distance (GSD) of 5.26 mm/pixel, point spacing of less than 1 cm, and an average point density of 22,278 points/m 3 . The on-board satellite positioning system of the UAV provided a scaled, orientated, and georeferenced point cloud.
The points not belonging to the surface of the slope were removed using the CloudCompare software v. 2.10.2. Most of the dense vegetation was manually cropped out in the upper part of the slope but also some small bushes originating from inside the rock mass were deleted. This explains the empty areas without 3D points on the rocky slope surface. Points disconnected from the slope surface (i.e., in the air) were filtered out using the tool of Connected Components, applying octree level 9 (grid step = 0.101775) and a value of 5 for the minimum points per component. Finally, the area of interest was cropped to remove points belonging to the road in the bottom of the rocky slope and those belonging to soil and rock fragments in the surroundings. The final 3D point cloud contained 5,228,097 points (Fig. 6) and was converted to a.txt file for the discontinuity set characterization step.

Discontinuity set dip and dip direction
The parameters used in the DSE software for the semiautomatic extraction of the discontinuity sets are summarized in Table 2. The results obtained for the dip and dip direction, the number of clusters for each set, and the total number of points assigned to each set are shown in Table 3. Out of the 5,228,097 points, a total of 3,645,455 were classified into discontinuity sets. The other 1,582,642 points remained unassigned since they do not fall into the cone filter criteria (γ 2 ≤ 30°); in other words, the poles of the normal vectors from those points have a dip and dip direction > 30° to any of the discontinuity set poles. Figure 7 summarizes visually the results obtained in this step with the stereogram plot of the extracted discontinuity set and the point cloud classified in different colors according to the discontinuity sets. DS 1 is the most representative discontinuity set in the slope with the majority of points. It corresponds to the slope surface and also to the bedding plane of the geological layers in the study site. DS 3 and DS 4 correspond to the overhangs, and DS 2 is almost a sub-vertical discontinuity bounding the lateral sides of the overhanging surfaces.
The discontinuity sets that are represented by planar surfaces are well characterized by this method and the points on these surfaces are correctly assigned to their corresponding discontinuity sets. However, points belonging to less planar surfaces (i.e., due to waviness or roughness) can be assigned to two different discontinuity sets as noted in some parts for DS 2 and DS 3. This limitation can explain why DS 2 has a slightly higher difference of 14° in the dip direction and DS 3 of 12° in dip compared to the values obtained in the field (Table 4). This discrepancy can be explained by the fact that the field measurements used for comparison were obtained locally on a chosen surface of the slope rather than the average of various measurements. For the scope of this work, the field measurements were only used to guide in the semi-automatic extraction in the DSE software, and further statistical analyses to evaluate the computed results for the discontinuity sets were not made. In-depth studies for the comparison or validation of the results obtained by the DSE software and field measurements are presented by Riquelme et al. (2014a), Jordá Bordehore et al. (2017), and Menegoni et al. (2019).

Normal discontinuity spacing
The normal spacing between discontinuities was calculated automatically using the DSE software (Table 5). For each cluster, the normal distance to the closest cluster is measured. The values of spacing, therefore, correspond to pairs of clusters. This implies that for each discontinuity set, the total number of spacing values is equal to the total number of clusters minus 1. DS 3 and DS 4 have a higher mean normal spacing, whereas DS 1 has the lowest (Fig. 8). DS 1 represents the slope surface but also the bedding plane of the geological layers, and as a result, the mean spacing of this discontinuity set (0.311 m) represents on average the spacing between the bedding planes, which are corroborated by field observations as being of decametric order (0.20 to 0.30 m). The extreme minimum spacing of DS 1 (0.05 mm) is due to the similar D parameter of clusters 5 and 11, where the computation was performed. In other words, the clusters, in this case, are almost coplanar and therefore the normal distance is practically zero. If a k threshold had been used (k > 0), these clusters would have been merged and only one merged cluster would have been formed.
The mean value of the normal spacing of each discontinuity set was considered as a geomechanical parameter for the RMR b computation ("Slope stability analysis" section).

Discontinuity persistence
The persistence of each cluster of a given discontinuity set in the dip, dip direction, maximum length, and the area of the convex hull were computed in the DSE software, considering non-persistent discontinuity conditions (k = 0). As only the persistence along the length of the maximum chord was considered for the susceptibility assessment, the statistical analysis results are presented in Table 6 and histogram in Fig. 9. As expected, DS 1 has the maximum persistence value, since it represents the bedding plane. DS 2 and DS 3 have a similar range of mean persistence, whereas DS 4 has the lowest.
The mean value of the length of the maximum chord for DS 2, DS 3, and DS 4 was considered as a geomechanical parameter for the RMR b calculation ("Slope stability analysis" section). Since DS 1 is mainly represented by 6 major clusters out of 190 (Fig. 10), their mean value of maximum length was chosen (15.35 m) instead of  Only the segments with more than 100-mm length intersected by discontinuities (red) were considered the mean value considering all the 190 clusters together (0.99 m).

Slope stability analysis
For the computation of the RMR b , the geomechanical parameters were characterized by field observations (e.g., Uniaxial Compressive strength-UCS, weathering, groundwater) and scanlines (aperture, roughness, infilling), observations on the images (infilling), and on the 3D point cloud (RQD, spacing, persistence, aperture). The estimated UCS value of 50-100 MPa was obtained by using a geologic hammer to fracture specimens (limestone) and counting the required blows to break them (ISRM 1978), yielding from 2 to 5 blows.
For the RQD index calculation, three sections were selected on the RGB 3D point cloud (Fig. 6), along a profile oriented perpendicularly to the principal discontinuity set on the surface of the slope. The calculation of the RQD index followed the same procedure as for drill cores (Fig. 11) and yielded 88% for section 1, 72% for section 2, and 95% for section 3, resulting an average of 85%. To validate this approach, the mean normal discontinuity set spacing obtained using the DSE software (Table 5) was used in Eq. 2, providing J V = 8.72 discontinuities/m 3 . This value of the J V applied to Eq. 1 provided a RQD index of 88%, close to the average value obtained by our approach (85%).
The aperture of 9 out of 10 discontinuities which intersected the scanline presented a value greater than 5 mm, as well as the majority of the apertures in other areas of the rocky slope. For the roughness characterization, the only discontinuity set that presented a uniform description is DS 1, which was characterized as smooth. The other discontinuity sets presented variable roughness either in macro-scale (stepped, undulating, and planar) or microscale (rough, smooth), according to the definition by ISRM (1978) or considering the ranges of Joint Roughness Coefficient-JRC values (from 2 to 10, 14-16, and 18-20) according to the definition by Barton and Choubey (1977). Therefore, the intermediate description of slightly rough was chosen for DS 2, 3, and 4 for calculation of the RMR b .
The infilling was absent and no pronounced weathering was observed in any discontinuity neither the presence of groundwater in the rocky slope. The description of the five geomechanical parameters for the RMR b computation for each discontinuity set is presented in Table 7 and the corresponding rating in Table 8. Although the score of each discontinuity set is slightly different (63-75), all scores fall within the same class II of the rock mass, corresponding to good quality rock.
The input parameters for the SMR index computation and the results for each discontinuity set applying the SMRTool calculator are shown in Table 9. Although potential toppling failures were identified in the SMRTool, they were not considered for this case study since the visual inspection of the rocky slope did not show important evidence of these types of failures. The large extension of the bedding plane (being fully persistent) is where the majority of discontinuities terminate, preventing the discontinuities to penetrate further on the rock mass and creating the typical parallelepiped shaped blocks with sufficient height to form a toppling failure mechanism.
DS 2 and DS 3 were assigned with an SMR corresponding to class I (completely stable), DS 4 to class II (stable), and DS 1 to class III (partially stable). DS 1 configures a more unsafe situation compared to others because its orientation is parallel to the slope orientation, indicated by F 1 = 1. In other words, the bedding plane of the geologic layer has the same dip direction and angle as the sliding plane of the slope, which favors instability. It is worth noticing that DS 3 that represents overhanging configures a safer situation compared to DS 1, but only in terms of failure mechanism. This can be explained because DS 3 is not parallel to the slope orientation but rather perpendicular to it. However, this discontinuity set as well as DS 4 are tension cracks or overhangs and contribute to instability due to lack of support in the rock mass. In these areas, the detached blocks slide on the slope surface, indicating the potential occurrence of instabilities. Therefore, these discontinuity sets will be also considered as indicators (overhangs) in the susceptibility analysis.

Indicators
The results of each indicator used for the rockfall susceptibility assessment (I S -spacing, I P -persistence, I SMR -Slope  . 12 (a) RGB point cloud overlaid with overhanging surfaces (red-1) compared to non-overhanging surfaces (white-0). Closer view of the upper overhang (b) and the two lower overhang (c) areas in pink (1) Mass Rating Index, and I o -overhangs) are presented separately for visualization on the 3D point cloud, followed by their scoring and integration into the Rockfall Susceptibility Index according to Table 1 and Eq. (3). All the discontinuities extracted by the DSE software were used for this analysis except those belonging to DS 1. The reason is that these discontinuities are the most dominant in the slope and represent the sliding plane (Fig. 7), as the platy limestone detaches in platy layers. Thus, given that majority of the slope consists of sliding planes, the entire slope would result as unstable, not being useful to distinguish different susceptibility levels. The presence of overhangs however is a differential factor determining the falling of blocks. It means that the lack of support plays a major role in rock detachment for this particular case study. Wherever the slope surface profile is not parallel with the sliding planes it would be useful to include all the discontinuities. For I o , three areas with overhanging surfaces were identified on the left side of the rocky slope (Fig. 12): one in the upper part, almost above the height of 14 m (Fig. 12b), and two close to the bottom below the height of 3 m (Fig. 12c). The surfaces of the overhangs are well delineated in most of the cases, except for less flat surfaces as shown in few parts in the bottom. The discontinuity sets which contain overhanging discontinuities are the DS 3 with 11 discontinuities and the DS 4 with 9 discontinuities.
To account for a conservative scenario for the I S , the minimum normal spacing was considered whenever a discontinuity presented more than one value. Figure 13 shows a realistic example for discontinuity 26 of DS 2. The DSE software computed the normal spacing to discontinuity 118 on its right side (S 26-118 = 0.48 m) and to discontinuity 18 on the left side (S 26-18 = 2.26 m). The use of the minimum spacing for the discontinuity 26 (S 26-118 ) leads to a poorer rock mass quality (Fig. 13a) whereas the use of maximum spacing (S 26-18 ) leads to a higher rock mass quality (Fig. 13b).
After applying the minimum normal spacing for each discontinuity for DS 2, DS 3, and DS 4, the I S in a 3D space is obtained (Fig. 14a). The majority of the discontinuity sets with minimum normal spacing are located on the left part of the rock slope (bottom and top). Few discontinuities have a higher value of spacing greater than 2 m (yellowish in Fig. 14a) and are located in the central and right parts of the rock slope. To account for a conservative scenario for the I P , the maximum length of persistence computed by the DSE software was assigned to each discontinuity of DS 2, DS 3, and DS 4 (Fig. 14b). In comparison with the 3D spacing distribution, the discontinuities with higher persistence are in the central area of the slope. The rocky slope has SMR classes I (DS 2, DS 3) and II (DS 4) located mainly on the left part (Fig. 14c).

Rockfall susceptibility index
The indicators for rockfall susceptibility were scored based on the thresholds of Table 1 and summed up to provide the rockfall susceptibility index I RF Eq. (3). It is worth noticing that the SMR classes VI (unstable) and V (completely unstable) were not present in the rocky slope (Table 9) and thus the higher value of I RF is 7 instead of 8. Nevertheless, the proposed scoring of the I SMR allows a broader applicability of this methodology in rocky slopes where all the SMR classes are present.
The 3D visualization of the I RF indicates areas of higher susceptibility in the left part of the rocky slope (Fig. 15a). The areas of overhangs (bottom and upper part) show moderate to high susceptibility depending on the persistence and spacing of discontinuities (Fig. 15b, c). The persistence of the exposed surfaces is a useful indicator since it contributes Fig. 15 (a) Areas with moderate (3-5) and high (5-7) rockfall susceptibility index are mainly in the left part of the slope. Overhanging areas in the upper (b) and bottom (c), differentiated in terms of susceptibility due to spacing and persistence properties ◂ Fig. 16 Moderate susceptibility areas in dark yellow as lateral exposed surfaces of the rocky slope to the differentiation of susceptibility levels within the overhanging areas. In other words, the greater the extension of lack of support is (persistence of overhangs), the more prone to rockfall it will be compared to smaller areas, an indication of larger unstable volumes.
In the lateral view of the 3D rocky slope, there are several areas with moderate susceptibility represented by the discontinuities with higher persistence (Fig. 16), showing previous breakages of rocks that were once part of the continuous slope surface. This inference is based on the observation of the right part of the slope, where fewer lateral exposed surfaces are present, and as a consequence, the slope surface is more continuous.

Discussion
The validation of the rockfall susceptibility assessment was carried out qualitatively by visual inspection from field observation and images (Fig. 17). The areas of high and moderate susceptibility, according to the I RF score, correspond to the areas where bigger and more fallen blocks are observed, i.e., below the overhanging surfaces (Fig. 17b) and the bottom left of the foot of the slope (Fig. 17c). In contrast, the right and continuous part of the slope has fewer and smaller blocks compared to the left part. Thus, the approach proposed in this work for identifying rockfall source areas yields realistic results. Further quantitative and more accurate validation can be performed at a later stage, using change detection techniques to assess future rockfalls and compare their location with the detected sources.
The proposed methodology can be adapted for its application to other study areas in terms of the use and assessment of indicators, scoring, and integration into the proposed I RF . For instance, even though the studied rocky slope did not have SMR class VI and V, it was included in the scoring of indicators (lower susceptibility). Although the studied rocky slope only presented planar failures, the other types (toppling and wedge) are partly incorporated into SMR and thus could be applied to slopes where these failures occur.
Adaptations required for the transferability of this methodology lie in the scoring of the indicators based on the thresholds for the I S and I P in areas of low (0), moderate (1), or high susceptibility (2), which are indeed site-specific. However, Fig. 17 (a) Orthoimage of the rock slope (white dashed line) generated with the UAV photos and the locations where bigger and more fallen blocks were found (orange boxes). (b) Closer view of the fallen blocks below overhanging surfaces. (c) Closer view of the bigger blocks located on the foot of the slope the higher weighting for longer persistence, small normal spacing discontinuities, presence of overhangs, and SMR classes VI and V should remain. One relevant note when using this approach elsewhere is that the identified areas of higher susceptibility do not guarantee that a rock detachment will occur, but rather that is it more likely compared to other areas of the same rocky slope based on referred indicators only. The study of the triggering mechanisms is not part of this work.
Concerning the intrinsic characteristics of the chosen slope, it is worth noticing that the most dominant discontinuity set (DS 1) has the majority of the indicators contributing to a higher rockfall susceptibility: (i) it is parallel to the slope which configures the worst-case scenario for planar failure computed by the SMR index, (ii) has the highest persistence (Table 6), and (iii) the minimum spacing (Table 5). For these reasons, almost the entire slope would be classified as highly susceptible whereas, in reality, only specific parts are, i.e., overhangs and lateral exposed surfaces (Fig. 16). Validation is required for rocky slopes with different discontinuity set characteristics to evaluate its performance.
This methodology has the advantage that the original 3D point cloud is used from beginning to end as new information of indicators is added as a scalar field, without the need for mesh generation and thus interpolation of results. The spacing and persistence indicators (I S and I P ) consider the local information of each cluster (discontinuity). A limitation of this work is that discontinuities represented by trace lengths are not used for the computation of the persistence and the spacing, and thus not incorporated as indicators for the susceptibility analysis.

Conclusions
This work presents a novel methodology for Rockfall Susceptibility Assessment for 3D slope models in the form of point clouds. It consists of four steps: (i) 3D slope model reconstruction using UAV digital photogrammetry, (ii) automatic characterization of discontinuity sets, (iii) slope stability analysis, and (iv) susceptibility assessment using a new Rockfall Susceptibility Index. The proposed approach can be used to refine the identification of potential rockfall source areas and to improve the input for hazard assessment, including rockfall run-out simulations. The methodology was applied to a rocky slope in the mountainous area of the Samaria National Park, in Crete Island, Greece. The identified areas of higher and moderate rockfall susceptibility in the 3D point cloud correspond to the areas where the bigger and highest number of fallen blocks were found on the foot of the slope, indicating that the methodology is efficient for detecting automatically potential rockfall sources.
The methodological workflow proposed in this work uses the original 3D point cloud from beginning to end as new information of indicators for susceptibility is added. This brings the advantage of no mesh generation and interpolation, which could compromise the quality of the 3D model. Moreover, the proposed methodological workflow can be automatized to be an efficient tool for rockfall source detection. Most of the input data needed to generate the rockfall susceptibility index are derived directly from the 3D point cloud, including the alternative approach to estimate the RQD index herein proposed and validated by the correlation to volumetric joint count J V of Palmstrom (2005). The only additional source of data is fieldwork data or image analysis to complement the discontinuity characterization (i.e., aperture, roughness, infilling, weathering, and groundwater).
This work also provides an attempt to visualize normal spacing, persistence, and the SMR index information in a point cloud environment. Additional innovation relates to the incorporation of the extension of overhangs in the rockfall susceptibility assessment by considering the persistence of exposed surface as maximum length. For the overhanging areas, this procedure allows the differentiation of parts with a larger lack of support them others, and as a consequence an indication of a greater volume of rock that is prone to rockfall. Further work should focus on the validation of this methodology in rocky slopes with different characteristics in terms of discontinuity set configuration (i.e., dominant discontinuity set not parallel to slope surface), or in the case that different SMR classes are present (I to VI), and where toppling and wedge failure mechanism occurs.