UAV-Based Evaluation of Rockfall Hazard in the Cultural Heritage Area of Kipinas Monastery, Greece

: Rockfall events consist one of the most hazardous geological phenomena in mountainous landscapes, with the potential to turn catastrophic if they occur near an anthropogenic environment. Rockfall hazard and risk assessments are recognized as some of the most challenging surveys among the geoengineering society, due to the urgent need for accurate foresight of likely rockfall areas, together with their magnitude and impact. In recent decades, with the introduction of remote sensing technologies, such as Unmanned Aerial Vehicles, the construction of qualitative and quantitative analyses for rockfall events became more precise. This study primarily aims to take advantage of the UAV’s capabilities, in order to produce a detailed hazard and risk assessment via the proposition of a new semi-quantitative rating system. The area of application is located in the cultural heritage area of Kipinas Monastery in Epirus, Greece, which is characterized by the absence of pre-existing data regarding previous rockfall events. As an outcome, it was shown that the suggested methodology, with the combination of innovative remote sensing technologies with traditional engineering geological ﬁeld surveys, can lead to the extraction of all the necessary quantitative data input for the proposed rating system for any natural slope.


Introduction
Rockfalls, in general, are classified as a type of landslide that consists of a detachment of a rock block (or several individual blocks) from a vertical or even sub-vertical cliff, followed by rapid downward motion characterized by four main phases: free-falling, bouncing, rolling and sliding [1]. Although rockfalls have a much lower level of economic risk than large-scale landslides, the rapid velocities associated with rockfall events result in a similar number of fatalities (in the same order of magnitude) as ones caused by all other landslide types [2]. Therefore, rockfall events consist one of the most hazardous geological phenomena in mountainous landscapes, with the potential to turn catastrophic if they occur near an anthropogenic environment. Thus, the accurate foresight of likely to rockfall areas and also their magnitude and impact is of high significance. In general, the above terms can be expressed with the parameters of susceptibility, hazard and risk, respectively.
Firstly, susceptibility is described as the likelihood of an event to occur in a specific area based on the local terrain conditions [3]. The term susceptibility represents the natural tendency of an area to be affected by a potential future hazardous event, and results in a spatial-based estimation of where rockfalls are more likely to occur [4][5][6][7].
In rockfall terminology, rockfall hazard refers to the probability of occurrence of a hazardous event with a specific magnitude or severity within a given area during a

Engineering Geological Field Survey
The first and most important method used was an in-depth engineering geological survey in the area of interest. The necessity of detailed and intense field work is derived from the urgent need to gain a primary understanding of the engineering geological model of the site. This enables the identification of failure modes and possible volumes and of their possible impact. With the engineering geological model established, the evaluation of all the computer-based measurements from the UAV-derived data can be achieved and potential errors can be revised. Furthermore, via the field survey, many important parameters for the hazard assessment, which could not be obtained remotely, were estimated.
From an engineering geological point of view, the limestone rock mass is moderately jointed, intersected by numerous major vertical fractures, which ultimately form the local face of the cliff. As mentioned before, the limestone is karstified and voids of large dimensions are formed, undermining the rock slope.

Engineering Geological Field Survey
The first and most important method used was an in-depth engineering geological survey in the area of interest. The necessity of detailed and intense field work is derived from the urgent need to gain a primary understanding of the engineering geological model of the site. This enables the identification of failure modes and possible volumes and of their possible impact. With the engineering geological model established, the evaluation of all the computer-based measurements from the UAV-derived data can be achieved and potential errors can be revised. Furthermore, via the field survey, many important parameters for the hazard assessment, which could not be obtained remotely, were estimated.
From an engineering geological point of view, the limestone rock mass is moderately jointed, intersected by numerous major vertical fractures, which ultimately form the local face of the cliff. As mentioned before, the limestone is karstified and voids of large dimensions are formed, undermining the rock slope.
According to the findings of the engineering geological field survey, two (2) main orientations were detected for the slope's face, due to its constant dip direction alternations, alongside the investigated area. Additionally, mainly three (3) major discontinuity sets and the bedding were identified, as presented in the Schmidt stereographic projection in Figure  2. The three major discontinuity sets are characterized as very steep in general (dip angle ≥80°) and two of them are transversal to the slope plane; thus, they form the rock slope face in most places. The bedding is presented as very mild up to horizontal in some (dip angle ≤15°) sub-regions of the slope. According to the findings of the engineering geological field survey, two (2) main orientations were detected for the slope's face, due to its constant dip direction alternations, alongside the investigated area. Additionally, mainly three (3) major discontinuity sets and the bedding were identified, as presented in the Schmidt stereographic projection in Figure 2. The three major discontinuity sets are characterized as very steep in general (dip angle ≥ 80 • ) and two of them are transversal to the slope plane; thus, they form the rock slope face in most places. The bedding is presented as very mild up to horizontal in some (dip angle ≤ 15 • ) sub-regions of the slope.
It is worth noting that due to the difficult terrain of the study area and the large scale of the hazardous regions, it was decided that the capabilities of the UAV should be taken advantage of, thus the estimation of both the spacing and the distance of the discontinuities were taken from the 3D point cloud model in accordance of course with the field-based observations. The spacing of the discontinuities has relatively high values (more than 5 m only in specific locations), hence the sizes of the rock blocks vary from moderate to very large (1.5 m 3 up to larger than 5.5 m 3 ). In places, mainly due to stress relief, the size of the rock blocks is significantly smaller (0.5 m 3 up to 2.5 m 3 ), especially on the eastern part of the slope. The distance of the discontinuity planes varies significantly depending on the degree of fracturing of the rock. In general, it varies in a range between 2 m and 5 m, although locally the distance can be lower than 2 m. It is worth noting that due to the difficult terrain of the study area and the lar of the hazardous regions, it was decided that the capabilities of the UAV should b advantage of, thus the estimation of both the spacing and the distance of the disc ties were taken from the 3D point cloud model in accordance of course with th based observations. The spacing of the discontinuities has relatively high value than 5 m only in specific locations), hence the sizes of the rock blocks vary from m to very large (1.5 m 3 up to larger than 5.5 m 3 ). In places, mainly due to stress relief, of the rock blocks is significantly smaller (0.5 m 3 up to 2.5 m 3 ), especially on the part of the slope. The distance of the discontinuity planes varies significantly dep on the degree of fracturing of the rock. In general, it varies in a range between 2 m m, although locally the distance can be lower than 2 m.
The joint wall compressive strength is high (estimated up to 24 MPa), w roughness of the discontinuity planes can be characterized as rough (JRC ranges b 8 and 10). In general, there is no infilling material in the discontinuities' surfaces basic angle of friction along the discontinuities was adopted to be equal to φ = 30 Studying the engineering geological conditions that control the investigate the face of the slope is generally characterized by relaxation of the rock mass, du high inclination. In conjunction with the spacing and the distance of discontinuitie erate-to-large prone blocks are formatted. The main reason that these instabilities widespread across the slope is the "nature" of the rock mass itself and its super chanical properties, which restrain the prone-to-failure blocks, causing a tempora ity on the slope.

Development of the 3D Model
Considering the unique location of the monastery and the steep morpholog rock slope, it was decided to conduct a photogrammetry survey via a manuall UAV in September of 2020, in order to construct an accurate and detailed 3D mod investigated site.
The camera view direction for the UAV-based photogrammetric survey was mented as [64] proposed for steep slopes. In [64], it is suggested that for steep slope The joint wall compressive strength is high (estimated up to 24 MPa), while the roughness of the discontinuity planes can be characterized as rough (JRC ranges between 8 and 10). In general, there is no infilling material in the discontinuities' surfaces and the basic angle of friction along the discontinuities was adopted to be equal to ϕ = 30 • .
Studying the engineering geological conditions that control the investigated slope, the face of the slope is generally characterized by relaxation of the rock mass, due to its high inclination. In conjunction with the spacing and the distance of discontinuities, moderate-to-large prone blocks are formatted. The main reason that these instabilities are not widespread across the slope is the "nature" of the rock mass itself and its superior mechanical properties, which restrain the prone-to-failure blocks, causing a temporal stability on the slope.

Development of the 3D Model
Considering the unique location of the monastery and the steep morphology of the rock slope, it was decided to conduct a photogrammetry survey via a manually flown UAV in September of 2020, in order to construct an accurate and detailed 3D model of the investigated site.
The camera view direction for the UAV-based photogrammetric survey was implemented as [64] proposed for steep slopes. In [64], it is suggested that for steep slopes (slope angle > 40 • , usually rock slopes), an oblique or even horizontal camera view should be applied in order to achieve the optimal results during the acquisition of the imagery data.
The UAV gear, via which the imagery data were collected, was a DJI Phantom 4 Pro V2.0 equipped with a 1-inch 20-megapixel sensor and an FC6310 camera of 8.8 mm focal length and a manually adjustable aperture from F2.8 to F11. A total of 362 images were acquired, in order to cover the whole slope face. The derived images from the UAV were processed using the Pix4Dmapper software v. 4.5.6. (Pix4D, Lausanne, Switzerland), resulting in a dense 3D point cloud dataset consisting of about 104 million points with a ground resolution of 1.27 cm/pixel and a reprojection error of 0.157 pixels. Further processing of the data enabled the construction of a detailed DEM (Digital Elevation Model) with the same pixel dimensions as the point cloud. The generated point cloud was georeferenced using the UAV sensor GPS data and a set of 5 artificial GCP (Ground Control Points) according to the methodology proposed by [41,65]. Considering the steep geometry of the slope and the relevant hazard to access higher parts, the GCP were placed in accessible lower parts across the road that leads to the monastery in some distinctive areas, as shown in Figure 3. acquired, in order to cover the whole slope face. The derived images from the UAV were processed using the Pix4Dmapper software v. 4.5.6. (Pix4D, Lausanne, Switzerland), resulting in a dense 3D point cloud dataset consisting of about 104 million points with a ground resolution of 1.27 cm/pixel and a reprojection error of 0.157 pixels. Further processing of the data enabled the construction of a detailed DEM (Digital Elevation Model) with the same pixel dimensions as the point cloud. The generated point cloud was georeferenced using the UAV sensor GPS data and a set of 5 artificial GCP (Ground Control Points) according to the methodology proposed by [41,65]. Considering the steep geometry of the slope and the relevant hazard to access higher parts, the GCP were placed in accessible lower parts across the road that leads to the monastery in some distinctive areas, as shown in Figure 3. Following the extraction, the point cloud was processed using the CloudCompare TM software v.2.11.1 (Anoia). In order to proceed to the structural analysis, the whole 3D model was trimmed of vegetation, constructions and "noise" in general. This prevented their presence in the final point cloud from affecting the correct recognition of discontinuities.
Then, the point cloud was segmented into four (4) different sub-regions depending on two (2) main criteria. The sub-regions were divided in order for each to preserve a consistent dip direction at the face of the slope and their downstream area to sustain a similar impact level, in the potential event of rockfall. This segmentation of the 3D point cloud into the 4 sub-regions is shown in Figure 4. The main characteristics of each area that the initial point cloud was segmented into are listed in Table 1. Following the extraction, the point cloud was processed using the CloudCompare TM software v.2.11.1 (Anoia). In order to proceed to the structural analysis, the whole 3D model was trimmed of vegetation, constructions and "noise" in general. This prevented their presence in the final point cloud from affecting the correct recognition of discontinuities.
Then, the point cloud was segmented into four (4) different sub-regions depending on two (2) main criteria. The sub-regions were divided in order for each to preserve a consistent dip direction at the face of the slope and their downstream area to sustain a similar impact level, in the potential event of rockfall. This segmentation of the 3D point cloud into the 4 sub-regions is shown in Figure 4. The main characteristics of each area that the initial point cloud was segmented into are listed in Table 1.

Structural Analysis
For each sub-region of the point cloud, a supervised automated analysis was performed with the open-source DSE (Discontinuity Set Extractor) software (developed by [66]). This approach, proposed by [66], calculates the normal vector for each point and its "knn" closest neighbor points. The visual representation of this procedure is via their corresponding poles in a stereographic projection. Then, statistical analysis was carried out for these poles, providing the stereogram and enabling the extraction of the principal discontinuity sets.
Specifically, the PCA (Principal Component Analysis) algorithm was used to calculate these normal vectors and a tolerance (n) test was performed in order to define if each set of points could be considered as a candidate of a planar area, or as non-planar and be discarded. Afterwards, in order to estimate the principal orientation, the Kernel Density Estimation [67] was used. The density function was estimated automatically, and with a manual identification of the peaks of the most represented orientations, the number and the orientation of principal planes was determined. Consequently, a discontinuity set was assigned to each input point from point cloud, enabling that way the isolation of all the co-belonging points to a specific set. Then, a spatial cluster analysis was performed with the DBSCAN (Density-Based Spatial Clustering of Applications with Noise) algorithm. This resulted in clusters of points, belonging to the same set, forming plane continuous surfaces which, hence, could be considered as planes, so that the automatic calculation of planes' equations (spacing [66] and persistence [68]) was enabled.

Structural Analysis
For each sub-region of the point cloud, a supervised automated analysis was performed with the open-source DSE (Discontinuity Set Extractor) software (developed by [66]). This approach, proposed by [66], calculates the normal vector for each point and its "knn" closest neighbor points. The visual representation of this procedure is via their corresponding poles in a stereographic projection. Then, statistical analysis was carried out for these poles, providing the stereogram and enabling the extraction of the principal discontinuity sets.
Specifically, the PCA (Principal Component Analysis) algorithm was used to calculate these normal vectors and a tolerance (n) test was performed in order to define if each set of points could be considered as a candidate of a planar area, or as non-planar and be discarded. Afterwards, in order to estimate the principal orientation, the Kernel Density Estimation [67] was used. The density function was estimated automatically, and with a manual identification of the peaks of the most represented orientations, the number and the orientation of principal planes was determined. Consequently, a discontinuity set was assigned to each input point from point cloud, enabling that way the isolation of all the co-belonging points to a specific set. Then, a spatial cluster analysis was performed with the DBSCAN (Density-Based Spatial Clustering of Applications with Noise) algorithm. This resulted in clusters of points, belonging to the same set, forming plane continuous surfaces which, hence, could be considered as planes, so that the automatic calculation of planes' equations (spacing [66] and persistence [68]) was enabled.
In this particular analysis, the used parameters for all 4 sub-regions were the same. Specifically, for the calculation of the normal vectors of each point and its corresponding poles, two parameters were defined; the number of nearest neighbors (knn = 30) and the tolerance for the coplanarity test (n = 0.2). For the calculation of the density of the poles, the number of bins for the kernel density estimation (nbins = 64) and the minimum angle between normal vectors of discontinuity sets (ϕ = 20 • ) were determined. Then, for the assignment of a discontinuity set to each point, the parameter that was imported was the minimum angle between the normal vector of a discontinuity set and the normal vector of the point (cone = 30 • ). Lastly, for the cluster analysis, an assumption was made about all cluster members of a discontinuity set having the same normal vector, and the parameter used for test if two clusters should be merged (ksigma = 1.5) was defined.
The discontinuity sets, for each sub-region, as they were semi-automatically extracted based on DSE, numbered four (4) and are presented visually in Figures 5-9 for each sub-region separately, while their main characteristics are listed in Table 2. of the point (cone = 30). Lastly, for the cluster analysis, an assumption was made cluster members of a discontinuity set having the same normal vector, and the p used for test if two clusters should be merged (ksigma = 1.5) was defined. The discontinuity sets, for each sub-region, as they were semi-automat tracted based on DSE, numbered four (4) and are presented visually in Figure  each sub-region separately, while their main characteristics are listed in Table 2.     It is important to underline that three of them comprise major joint sets (J1, J2 and J3) and one (1) the bedding (Β).

Kinematic Analysis
In general, from the field survey, it is preliminarily comprehended that the main failure mode in the whole area comprises rock block detachment followed mainly by toppling or wedge failure and secondarily by planar sliding failure. The fallen blocks reveal fresh surfaces and as a result, the lateral blocks exhibit a lack of confining stress. Consequently, new cracks emerge and the same process accelerates.
In order to determine the optimal failure model for each sub-region, a kinematic analysis was performed based on the calculated major discontinuity sets from DSE software, which were previously validated according to the field-obtained structural data [60,62]. It is important to underline that three of them comprise major joint sets (J1, J2 and J3) and one (1) the bedding (Β).

Kinematic Analysis
In general, from the field survey, it is preliminarily comprehended that the main failure mode in the whole area comprises rock block detachment followed mainly by toppling or wedge failure and secondarily by planar sliding failure. The fallen blocks reveal fresh surfaces and as a result, the lateral blocks exhibit a lack of confining stress. Consequently, new cracks emerge and the same process accelerates.
In order to determine the optimal failure model for each sub-region, a kinematic analysis was performed based on the calculated major discontinuity sets from DSE software, It is important to underline that three of them comprise major joint sets (J 1 , J 2 and J 3 ) and one (1) the bedding (B).

Kinematic Analysis
In general, from the field survey, it is preliminarily comprehended that the main failure mode in the whole area comprises rock block detachment followed mainly by toppling or wedge failure and secondarily by planar sliding failure. The fallen blocks reveal fresh surfaces and as a result, the lateral blocks exhibit a lack of confining stress. Consequently, new cracks emerge and the same process accelerates.
In order to determine the optimal failure model for each sub-region, a kinematic analysis was performed based on the calculated major discontinuity sets from DSE software, which were previously validated according to the field-obtained structural data [60,62]. This analysis was executed, for each area, via DIPS v.7 software (RocScience Inc., Toronto, ON, Canada). The geometrical characteristics of each slope and of each discontinuity set were imported, in addition to the friction angle of the bed rock, which was estimated as ϕ = 30 • , and a stereographic projection occurred for every failure type. These produced stereonet plots, identified any potential failures and highlighted the discontinuity set(s) which contributed to those. The results of this process are demonstrated in Table 3 and the stereographic projections for each sub-region are shown in Figures 10-13. Table 3. A concise review of the kinematic analysis per sub-region.

Sub-Region
Type of Failure Contributed Discontinuity Sets Appl. Sci. 2021, 11, x FOR PEER REVIEW 11 of 34 ON, Canada). The geometrical characteristics of each slope and of each discontinuity set were imported, in addition to the friction angle of the bed rock, which was estimated as φ = 30, and a stereographic projection occurred for every failure type. These produced stereonet plots, identified any potential failures and highlighted the discontinuity set(s) which contributed to those. The results of this process are demonstrated in Table 3 and the stereographic projections for each sub-region are shown in Figures 10-13. Table 3. A concise review of the kinematic analysis per sub-region.

Rockfall Susceptibility Analysis
After the generation of the optimal failure models for each sub-region, the monitoring and, subsequently, the quantification of two (2) important elements for this hazard assessment were made possible. These were the already failed rock blocks and the potentially unstable ones for each sub-region. With the spatial distribution of these blocks monitored, important measures were exported regarding their magnitude (volume, mass) and therefore, the frequency of each magnitude was calculated per sub-region. Additionally, the density of each category of rock blocks in every sub-region was estimated with the combination of each slope's size and the number of the monitored blocks.

Rockfall Susceptibility Analysis
After the generation of the optimal failure models for each sub-region, the monitoring and, subsequently, the quantification of two (2) important elements for this hazard assessment were made possible. These were the already failed rock blocks and the potentially unstable ones for each sub-region. With the spatial distribution of these blocks monitored, important measures were exported regarding their magnitude (volume, mass) and therefore, the frequency of each magnitude was calculated per sub-region. Additionally, the density of each category of rock blocks in every sub-region was estimated with the combination of each slope's size and the number of the monitored blocks.

Rockfall Susceptibility Analysis
After the generation of the optimal failure models for each sub-region, the monitoring and, subsequently, the quantification of two (2) important elements for this hazard assessment were made possible. These were the already failed rock blocks and the potentially unstable ones for each sub-region. With the spatial distribution of these blocks monitored, important measures were exported regarding their magnitude (volume, mass) and therefore, the frequency of each magnitude was calculated per sub-region. Additionally, the density of each category of rock blocks in every sub-region was estimated with the combination of each slope's size and the number of the monitored blocks.

Rockfall Susceptibility Analysis
After the generation of the optimal failure models for each sub-region, the monitoring and, subsequently, the quantification of two (2) important elements for this hazard assessment were made possible. These were the already failed rock blocks and the potentially unstable ones for each sub-region. With the spatial distribution of these blocks monitored, important measures were exported regarding their magnitude (volume, mass) and therefore, the frequency of each magnitude was calculated per sub-region. Additionally, the density of each category of rock blocks in every sub-region was estimated with the combination of each slope's size and the number of the monitored blocks.
The followed methodology was based initially on the approach in [69], but it was further developed to utilize the maximum of the detailed UAV-derived dataset's capabilities. Specifically, the proposed methodology started with the processing of every sub-region's point cloud, with the CloudCompare TM software v.2.11.1 (Anoia), in order to perform two separate semi-automated detections for the failed and the likely-to-fail rock blocks, respectively. These detections were based on the visual interpretation of the main discontinuity sets (DSE exports) and their contribution to each slope's failure model (Kinematic Analysis exports). Hence, the geometric characteristics, the spatial distribution, the magnitude and the density of presence of each rock block were specified. Due to the lack of available rock block volumes to be analyzed, the identification and the quantification of the already failed blocks were extracted via the analysis of the fresh surfaces across each sub-region's face slope, that attested to the existence of a prior failure. In total, 455 rock blocks were detected, 326 of them as already detached and 129 as potentially unstable. Their detailed allocation per sub-region is demonstrated in Table 4 and their visual interpretation is shown in Figure 14. The followed methodology was based initially on the approach in [69], but it was further developed to utilize the maximum of the detailed UAV-derived dataset's capabilities. Specifically, the proposed methodology started with the processing of every subregion's point cloud, with the CloudCompare TM software v.2.11.1 (Anoia), in order to perform two separate semi-automated detections for the failed and the likely-to-fail rock blocks, respectively. These detections were based on the visual interpretation of the main discontinuity sets (DSE exports) and their contribution to each slope's failure model (Kinematic Analysis exports). Hence, the geometric characteristics, the spatial distribution, the magnitude and the density of presence of each rock block were specified. Due to the lack of available rock block volumes to be analyzed, the identification and the quantification of the already failed blocks were extracted via the analysis of the fresh surfaces across each sub-region's face slope, that attested to the existence of a prior failure. In total, 455 rock blocks were detected, 326 of them as already detached and 129 as potentially unstable. Their detailed allocation per sub-region is demonstrated in Table 4 and their visual interpretation is shown in Figures 14.  A supervised automated analysis was then performed for the calculation of each block's volume via the 2.5D Volume Tool, provided by the CloudCompare TM software v.2.11.1 (Anoia). In order to estimate the credibility and highlight potential outliers in these automated calculations, the manual calculation of 20% of the total rock block vol- A supervised automated analysis was then performed for the calculation of each block's volume via the 2.5D Volume Tool, provided by the CloudCompare TM software v.2.11.1 (Anoia). In order to estimate the credibility and highlight potential outliers in these automated calculations, the manual calculation of 20% of the total rock block volumes for each sub-region via the extracted point clouds was decided. The followed methodology was inspired by relevant previous surveys, such as [70], and was further developed.
In this particular case, the volumes of both the already failed and the potentially unstable rock blocks varied. In general, their volumes fluctuated in moderate-to-large ranges (≤4 m 3 ), but it is noteworthy to highlight the existence of some individual cases in which these values were surpassed. Of course, these excesses have not occurred extensively along the slopes, but they comprise isolated cases. In particular, most of these isolated cases were identified in sub-region U3, where many of the potentially unstable block volumes were estimated to have greater values than the average.
On the contrary, the spatial distribution of the examined rock blocks (both the already failed and the potentially unstable) seems to follow certain directions across the face of each slope. This could easily be understood with the predetermined optimal failure models for each sub-region, as the formation of these rock blocks followed the directions of the discontinuity set(s) which contributed to any of the potential failures.
Specifically, for sub-region U1, the main potentially unstable block volumes that were identified, across the slope's profile, ranged from small to moderate with values between 0.5 m 3 and 1.6 m 3 . In sub-region U2, the potentially unstable block's volume distribution mainly fluctuated between 0.6 m 3 and 1.3 m 3 . In the third sub-region (U3), the volume distribution mode displayed its highest range for any of the investigated slopes, with values mainly varying from 1.2 m 3 up to 5.2 m 3 . Of course, as mentioned before, some of those volume excesses (>5.1 m 3 ) comprised isolated cases and they were not considered as a distinctive feature for this sub-region. At last, for sub-region U4, the main potentially unstable block volumes that were identified, across the slope's profile, varied from very small to small with values mainly between 0.2 m 3 and 1.1 m 3 .
Overall, the correlation of both the spatial distribution and the estimated volumes between the already failed and the potentially unstable rock blocks seemed successful by the results of this analysis. Thus, the identification and the quantification of the potentially unstable rock blocks, which had the biggest margin of error in this particular case, appeared to be well-established. In Figures 15-18 is shown the volume distribution (frequency %) for each investigated sub-region analytically and in Figure 19 the cumulative volume distribution for the whole investigated site.
umes for each sub-region via the extracted point clouds was decided. The followed methodology was inspired by relevant previous surveys, such as [70], and was further developed.
In this particular case, the volumes of both the already failed and the potentially unstable rock blocks varied. In general, their volumes fluctuated in moderate-to-large ranges (≤4 m 3 ), but it is noteworthy to highlight the existence of some individual cases in which these values were surpassed. Of course, these excesses have not occurred extensively along the slopes, but they comprise isolated cases. In particular, most of these isolated cases were identified in sub-region U3, where many of the potentially unstable block volumes were estimated to have greater values than the average.
On the contrary, the spatial distribution of the examined rock blocks (both the already failed and the potentially unstable) seems to follow certain directions across the face of each slope. This could easily be understood with the predetermined optimal failure models for each sub-region, as the formation of these rock blocks followed the directions of the discontinuity set(s) which contributed to any of the potential failures.
Specifically, for sub-region U1, the main potentially unstable block volumes that were identified, across the slope's profile, ranged from small to moderate with values between 0.5 m 3 and 1.6 m 3 . In sub-region U2, the potentially unstable block's volume distribution mainly fluctuated between 0.6 m 3 and 1.3 m 3 . In the third sub-region (U3), the volume distribution mode displayed its highest range for any of the investigated slopes, with values mainly varying from 1.2 m 3 up to 5.2 m 3 . Of course, as mentioned before, some of those volume excesses (>5.1 m 3 ) comprised isolated cases and they were not considered as a distinctive feature for this sub-region. At last, for sub-region U4, the main potentially unstable block volumes that were identified, across the slope's profile, varied from very small to small with values mainly between 0.2 m 3 and 1.1 m 3 .
Overall, the correlation of both the spatial distribution and the estimated volumes between the already failed and the potentially unstable rock blocks seemed successful by the results of this analysis. Thus, the identification and the quantification of the potentially unstable rock blocks, which had the biggest margin of error in this particular case, appeared to be well-established. In Figures 15-18 is shown the volume distribution (frequency %) for each investigated sub-region analytically and in Figure 19 the cumulative volume distribution for the whole investigated site.        After the identification of the population, the magnitude and the spatial distribution of the potentially unstable rock blocks with a 3D approach, a rockfall analysis based on trajectometry simulations was executed via RocFall software v.8.0.12 (RocScience Inc., Toronto, ON, Canada), in order to predict the potential rockfall hazards in the area. RocFall software can be characterized as a 2D statistical analysis program designed to assist with the assessment of slopes at risk for rockfalls. Energy, velocity and "bounce height" were determined by using the software for each sub-region. Distributions and comprehensive statistics of these parameters were automatically calculated along each slope's profile and projected as graphs. This software can also assist in determining remedial measures. Information about the kinetic energy and location of impact on a barrier can help determine the effectiveness (capacity, size and location) of barriers [71]. The barriers can be custom, user-defined or predefined from the software's library. As a result, the extracted information is vital for the estimation of each sub-region's hazard level, based on the proposed Rockfall Hazard and Risk Rating System. After the identification of the population, the magnitude and the spatial distribution of the potentially unstable rock blocks with a 3D approach, a rockfall analysis based on trajectometry simulations was executed via RocFall software v.8.0.12 (RocScience Inc., Toronto, ON, Canada), in order to predict the potential rockfall hazards in the area. RocFall software can be characterized as a 2D statistical analysis program designed to assist with the assessment of slopes at risk for rockfalls. Energy, velocity and "bounce height" were determined by using the software for each sub-region. Distributions and comprehensive statistics of these parameters were automatically calculated along each slope's profile and projected as graphs. This software can also assist in determining remedial measures. Information about the kinetic energy and location of impact on a barrier can help determine the effectiveness (capacity, size and location) of barriers [71]. The barriers can be custom, Figure 19. Volume distribution for both the fallen (blue) and the potentially unstable (red) blocks for the whole slope.
In general, a detached block descends very rapidly through the air by falling, bouncing and rolling [72], while its trajectory mainly depends on the geometry of the slope. In order to model the trajectory of the falling rocks and design the rockfall barriers, it was necessary to determine some parameters for all four (4) sub-regions. These were the normal and tangential coefficients of restitution (R n and R t , respectively) of the bed rock and each slope's most critical profile. The values of two coefficients of restitution (R n and R t ) for the limestone were based on previous publications, for limestone rocks found in similar geological settings with similar properties as the one in the investigated site, estimated as 0.5 and 0.5, respectively [73]. The critical slope profile for each of the four (4) sub-regions was extracted from the georeferenced 3D point-cloud via CloudCompare TM software v.2.11.1 (Anoia). Specifically, a polyline was constructed for each profile followed the morphology. Then, from the attribute table, the 3D coordinates of every point, that consisted of each polyline (X-Y-Z), were converted to 2D (X-Y). The 2D information could then be imported to RocFall software v.8.0.12 (RocScience Inc., Toronto, ON, Canada) in order to achieve the most precise visual interpretation of each slope's profile.
In this particular case, for sub-region U1 three (3) main rockfall release points were identified across the slope's profile. The result of the simulation indicates that the failure from two (2) of them (point a and point b) cannot threaten the monastery, and, consequently, the rockfall risk can be neglected. The third identified release point can be harmful for the monastery and the maximum total kinetic energy of the rockfall at the location of impact with the monastery was estimated as 69 kJ. In sub-region U2, there were identified two (2) main rockfall release points that could both, potentially, produce threatening events for the road below the slope. The maximum total kinetic energy of the rockfall at the location of impact with the road was estimated to be 429 kJ for the first release point (point a) and 747 kJ for the second (point b). In the third sub-region (U3), two (2) main detachment points for rockfalls were identified as potentially harmful (point a and point b), and one (1) cannot threaten the road below (point c). The maximum total kinetic energy of the rockfall at the location of impact with the road was estimated as 3621 kJ for the first release point (point a) and 1386 kJ for the second (point b). These were the highest values for any of the investigated slopes. At last, in sub-region U4, two (2) main rockfall release points were identified. Both of them could potentially threaten the road below the slope. The maximum total kinetic energy of the rockfall at the location of impact with the road was estimated as 290 kJ for the first release point (point a) and 102 kJ for the second (point b). The detailed results from this analysis, per sub-region, are demonstrated in Table 5, while a visual interpretation of the rock blocks trajectories is presented in Figure 20.

General
The assessment of rockfall hazards and risks along roads and in other human-related locations is of great importance [24]. In order to determine the value of these parameters, a number of rating systems have been developed based on either qualitative or quantitative criteria. These approaches differ entirely in terms of input data, adopted analysis procedure and of course the final outputs.
Qualitative-based methods define hazards, elements at risk and their vulnerabilities with the usage of qualitative descriptors, such as ranked attributes, weighted indices, rating systems, scoring schemes and classifications [74]. This type of evaluation can be based on either objective (statistical or mathematical) estimates, subjective (judgments or assumptions) estimates, or even a combination of both. Detailed ratings usually use several judging factors to provide an evaluation and/or score for the slope. The results of this

General
The assessment of rockfall hazards and risks along roads and in other human-related locations is of great importance [24]. In order to determine the value of these parameters, a number of rating systems have been developed based on either qualitative or quantitative criteria. These approaches differ entirely in terms of input data, adopted analysis procedure and of course the final outputs.
Qualitative-based methods define hazards, elements at risk and their vulnerabilities with the usage of qualitative descriptors, such as ranked attributes, weighted indices, rating systems, scoring schemes and classifications [74]. This type of evaluation can be based on either objective (statistical or mathematical) estimates, subjective (judgments or assumptions) estimates, or even a combination of both. Detailed ratings usually use several judging factors to provide an evaluation and/or score for the slope. The results of this method usually are expressed using relative terms, such as high, moderate and low. Qualitative analyses provide useful information for hazard and risk management, for relative comparisons of different sites and for facilitating the prioritization of mitigation measures [75].
On the contrary, quantitative-based methods use numerical values or ranges of them instead of qualitative relative terms. Generally, quantitative hazard and risk analyses approach the evaluation of hazard in the form of a numerical probability method, which estimates the frequency of failures, including the intensity of a specific spatial distribution and attempts to quantify risk with the monetary expression of the potential damage. Even though the hazard evaluation, in qualitative terms, seems prevalent and the reliability of its results occur from the availability, quantity, quality and credibility of the derived data, the same do not apply for quantitative risk analysis. Several recent efforts have attempted to establish standard procedures for quantifying risk in terms of official national recommendations or guidelines [76]. In most quantitative risk analyses, in order to quantify the probability of potential losses that are related to the occurrence of a hazardous event, such as rockfalls, the parameters that are taken into consideration are the number of demolished structures, injured people and fatalities. However, this does not constitute a standard procedure for every risk analysis because it depends on the researcher's approach. Furthermore, the significance of the exported results varies as the acceptable and tolerable risk thresholds differ between countries and regulations.
In general, qualitative analyses are more commonly used because they are easy to use and can be performed quickly [77]. For this reason, most rockfall hazard and risk assessments start with a qualitative approach to the investigated site as a whole and then adopt a detailed quantitative character, with thorough ratings to numerically differentiate the hazard and/or risk, respectively, at the predetermined specific region of the site. In other words, qualitative analysis represents an initial analysis step of the dominant hazards at a given site to evaluate the most dangerous areas further via quantitative methods. Despite this, nowadays, quantitative rockfall hazard and risk assessments are only possible when a detailed historical inventory of the rockfall events that have occurred in the study area is available [78]. However, the documentation of prior rockfall events is usually poorly made or even absent at most sites. Furthermore, complete rockfall inventories are rare in most parts of the world, due to the lack of reporting of small and medium events, which are usually not related to significant impacts [79]. Moreover, due to shortage of available data (geomechanical, statistical and historical) of previous rockfall activities, the accuracy of the results in quantitative-based assessments varies. The accuracy of a rockfall hazard and risk assessment depends on the quality of the "components" that constitute the analysis, and their processing refinement, whether the data consist of qualitative relative descriptors or absolute numerical values. Hence, a quantitative evaluation may not be more precise than a qualitative one [80].
Consequently, determining which method best suits a specific site is not an easy task. For this paper, a qualitative-based approach was implemented and a new rockfall hazard and risk rating system was proposed.

Background
The first rockfall hazard assessment methodology was proposed for railway networks by [81] in 1976. This method was based on a qualitative approach to evaluate the hazard level and prioritize the remediation need with alphabetical descriptors from "A" (the most hazardous regions which require immediate intervention) up to "F" (the most stable regions which do not require any intervention). This first attempt at a rating system was initially applied along the Canadian Pacific Railroad [81] and was then used as the cornerstone for any rating system methodology, related to rockfall hazards and risks, ever developed afterwards.
One of the first and most well-known rockfall hazard rating methods is called the Rockfall Hazard Rating System (RHRS) and was developed by [82] in Oregon (USA), in order to provide the local Department of Transportation with a simple, repeatable and easy-to-use standard rockfall hazard system along railways [83]. This system is considered as the basis for most rating systems, which further developed it in order to satisfy the requirements of variable transportation departments. Just to mention a few, the most widely used modified RHRS methods are: the Rock Slope Rating Procedure for New York (USA) [84][85][86][87], the mRHRS for Italy [24], the MORFH RS for Missouri (USA) [88][89][90], the TRHRS for Tennessee (USA) [91][92][93][94], the CRHRS for Colorado (USA) [95][96][97] and the RHRON for Ontario (Canada) [98][99][100].
Other important rating systems include the Rock Slope Hazard Index (RSHI) system, which was developed for rock slope inspection in highways in Scotland from [101]. The authors in [10] applied a quantitative hazard and risk methodology for the highways and railways in British Columbia, which was based on the magnitude and frequency distribution of rockfall events. The authors in [102] also proposed a new rockfall hazard rating system for use along a railway, named Quantitative Risk Assessment (QRA). This system was developed to quantify the probability of potential losses that are related to a hazardous event, with absolute results that can be used to compare different locations, thus providing a benchmark for prioritizing mitigation measures [102]. To this day, QRA constitutes one of the most popular quantitative methodologies.
A common ground between all the systems that were mentioned above is that all of them are orientated to anthropogenic cut slopes for highways and railways and not to natural ones. Some researchers have proposed systems that are applicable to natural slopes, but most of them are constructed for specific locations only. The authors in [103] proposed a hazard estimation methodology, specifically for a part of the transportation network across Central Italy, with 3D rockfall trajectories and their spatial distribution in order to estimate the risk. In a related way, [104] presented another quantitative approach with 3D rockfall trajectories for risk evaluation at the Sola d'Andorra slope (Andorra Principality) before and after the implementation of risk mitigation works. Lastly, [105] proposed a qualitative rockfall risk rating system, applicable for the estimation of risk for both natural and anthropogenic slopes, which was developed for a case study at the historical site of Monemvasia, Greece.

Rockfall Hazard and Risk Rating System
In the present study, a qualitative rockfall hazard and risk rating system is proposed. This rating system was designed with an emphasis on the rating of natural rock slopes, which indicate a susceptibility to rockfall events with a significant impact on human structures and activities. The above criteria in this particular case were fulfilled by the investigated site, the cultural heritage area of Kipinas Monastery in Epirus, Greece.
In general, four (4) major parameter categories are defined, which consist of a variety of different ranking values, in order for the total risk level for each sub-region to be estimated. Three (3) of these categories are related to the hazard level and one (1) to the impact of a potential rockfall event. In this assessment, every rating parameter has a different weight, which is based on the significance of each individual contribution to either the hazard and/or the risk level for every sub-region. The four main categories are presented in Table 6. Table 6. The main categories of rating factors and their weight for the proposed rating system.

Major Category Weight in Risk Level (%)
A. Geometry Factors 30

B. Geoengineering Factors 30
C. Triggering Factors 10 D. Impact Factors 30 In particular, the first category (A) of the rating system is related to the geometry and the morphology of the slope in general. It involves parameters such as the height, the angle, the roughness and the vegetation of the slope, the main height of the rockfall source area and the availability of obstacles and vegetation in the catchment zone at the footway of the slope. Category B involves the engineering geological factors of each slope. These factors describe the condition of the rock discontinuities, the intact rock strength, the rock mass structure, the weathering level of the rock mass, the spatial distribution of the predetermined rock blocks (already fallen and prone) per slope, the volume of the dominant rock blocks (already fallen and prone) and the potential kinetic energy of a rockfall at the location of impact with the ground. The parameters in the third major category (C) are related to the potential triggering factors, and feature the rainfall conditions and their severity, the drainage conditions of the rock mass and the seismicity of the area. The last category (D) refers to the impacts of a potential rockfall event and contains factors related to the human presence and the value of the pre-existing structures. As mentioned before, the weight for each category varies, depending on the importance of the parameters involved. Specifically, in this assessment, category A, B and D are given a weight of 30% each and category C 10% to the total risk score for every investigated sub-region. The factors of each category and their individual weights are presented in Table 7.  Each examined factor can be evaluated between 10 and 100, as the conditions evolve from favorable to adverse, and then is multiplied by a respective predetermined weight. The total hazard score is calculated by summing the individual score of each factor, from the first three (A, B and C) major categories. At last, the total risk score is calculated by adding the fourth category (D) to the hazard score. Therefore, a slope with the highest hazard level will have a total weighted score of 70 and, one with the highest risk level a total weighted score of 100. The final result of this assessment is the construction of two (2) separate classifications, one hazard and one risk related, with five categories each (very low up to very high) for all four distinct sub-regions in the investigated site. Based on these classifications, both a hazard and a risk zonation of the whole investigated cliff against rockfall occurrence were made possible to be constructed. These classifications are shown in Tables 8 and 9. Lastly, it is important to highlight that the proposed rating system was developed on an empirical basis, with the involved parameters, their rating and their weight ranges based on the available data, reasonable geoengineering judgment and the known facts.

Results
The proposed Rockfall Hazard and Risk Rating System was applied at the four (4) predetermined sub-regions (U1 to U4, Figure 3) along the rock cliff of Kipinas Monastery. Since these locations are in the same area, some of the factors (mostly spatial related) have the same rating for all four, whereas some differ for each slope section. In particular, 12 out of 23 factors have the same rating. After extensive application of this rating system to the investigated site, the total scores were categorized, according to Table 8 for hazard and  Table 9 for risk estimation. The final result of the application is a hazard and risk zonation of the cliff against rockfall occurrence, presented in Table 10 and in Figure 21.  Sub-region U1 is categorized in the moderate hazard zone, mainly due to the low height of the potential rockfall release zone and also due to the small sizes of the prone blocks. However, due to the high human presence and the high structure of the involved structure (in this sub-region is located the monastery), the total risk is high. The hazard levels in the area U2 are high, mainly as a result of the moderate-to-large volume of the prone blocks and their notable total kinetic energy at the location of impact with the road. However, the lack of significant structures in this location (only a road is located in U2) and the medium human presence (location of a resting area along the main road) result in a moderate total risk level. Sub-region U3 is the only one with both a high hazard level and high risk level. This occurs primarily from the notable large volume of the prone rock blocks which could be released from a significant height (47 m) which is greater than 2/3 of the total slope height. Additionally, the residencies of the monks are located in U3; the impacts of a potential rockfall event in this area would be very important. In the last investigated location (U4), both the hazard level and the risk level are moderate. The size of the prone unstable blocks is not large, in contrast with the dense spatial distribution shown, which occurs from the small dimensions of the area. Additionally, in this area there is a catchment zone, even if it is very restricted. The value of the structures is low, because only the road is located there, and the human presence is below moderate levels.

Discussion
In this study, a new Rockfall Hazard and Risk Rating System has been proposed and the implementation of a methodology has been suggested in order for all the critical parameters of the system to be estimated. Through this, a semi-quantitative rockfall hazard and risk assessment was performed.
As presented above, many rating systems for rockfall assessment have been pro- Sub-region U1 is categorized in the moderate hazard zone, mainly due to the low height of the potential rockfall release zone and also due to the small sizes of the prone blocks. However, due to the high human presence and the high structure of the involved structure (in this sub-region is located the monastery), the total risk is high. The hazard levels in the area U2 are high, mainly as a result of the moderate-to-large volume of the prone blocks and their notable total kinetic energy at the location of impact with the road. However, the lack of significant structures in this location (only a road is located in U2) and the medium human presence (location of a resting area along the main road) result in a moderate total risk level. Sub-region U3 is the only one with both a high hazard level and high risk level. This occurs primarily from the notable large volume of the prone rock blocks which could be released from a significant height (47 m) which is greater than 2/3 of the total slope height. Additionally, the residencies of the monks are located in U3; the impacts of a potential rockfall event in this area would be very important. In the last investigated location (U4), both the hazard level and the risk level are moderate. The size of the prone unstable blocks is not large, in contrast with the dense spatial distribution shown, which occurs from the small dimensions of the area. Additionally, in this area there is a catchment zone, even if it is very restricted. The value of the structures is low, because only the road is located there, and the human presence is below moderate levels.

Discussion
In this study, a new Rockfall Hazard and Risk Rating System has been proposed and the implementation of a methodology has been suggested in order for all the critical parameters of the system to be estimated. Through this, a semi-quantitative rockfall hazard and risk assessment was performed.
As presented above, many rating systems for rockfall assessment have been proposed during the past decades, but nowadays, with the rise of remote sensing technologies in the field of geohazards, it is believed that there is still margin for improvement. Therefore, the construction of a new rating system based in quantitative parameters, exported from the UAV-derived data, was considered as an appropriate evolution in rockfall hazard and risk assessments.
The qualitative-based approach (i.e., rating system) was selected over a quantitative one, in order to cover a wider range of applications in areas susceptible to rockfall events with negligible to nonexistent data from former rockfall-related events. The main target of the proposed methodology is not to replace the detailed quantitative analyses of rockfall events, but to be considered as the initial step towards rockfall hazard and risk assessments with the implementation of it to the investigated site as a whole in order to classify it to sub-regions regarding their hazard and/or risk level, respectively. Through its application, useful information for hazard and risk management will be provided, for relative comparisons of different sites and also for facilitating the prioritization of mitigation measures. Then, of course, if it is deemed necessary, a detailed quantitative character, with thorough ratings to numerically differentiate the hazard and/or risk, respectively, could be adopted at the already predetermined susceptible regions of the site [75].
The main point that sets the current approach apart from its predecessors is the type of the value of the implemented factors. Generally, rating systems constitute a qualitative type of assessment, but this does not mean that the value of the factors that comprise them have to be qualitative-based. The main disadvantage of qualitative-based values is that they are subjective. This means that in the same area, with the same available data, the results could vary significantly between different researchers. On the contrary, the value of the factors of this proposed rating system are quantitative-based. The UAV-derived data, in accordance with the appropriate computer-based processing, enabled the calculation of important parameters in order for the properties of the complicated phenomena of rockfalls to be understood and expressed in absolute values. Factors regarding the Total Kinetic Energy, the rock block distribution across the slopes (fallen and potentially unstable) and the mode rock block volumes (fallen and potentially unstable) that are absent in most of the prior well-established rating systems such as RHRS [82] and its modifications [24,[84][85][86][87][88][89][90][91][92][93][94][95][96][97][98][99][100], RHI [11] and the rockfall risk rating systems proposed by [105], comprise the cornerstone of this newly proposed Rockfall Hazard and Risk Rating System.
The methodology that was followed in this particular study and the parameters that were decided to form the rating system may seem to be related to a quantitative-based analysis at first glance. Specifically, the detection of potentially unstable rock blocks, the estimation of their potential volumes and the probabilistic analysis of them based on trajectometry simulations are very common approaches that are followed in quantitative analysis, in sites with extensive rockfall-related multi-temporal inventories [106][107][108]. In such cases, a multi-temporal change detection procedure between the different datasets is performed, in order for topographical changes between different epochs to be monitored and characterized for modeling of the rockfall's dynamics. Consequently, the aforementioned values can be estimated. In this study, the identification of the population, the magnitude and the spatial distribution of the potentially unstable rock blocks was executed with a 3D approach, based on the accuracy of the point cloud and the knowledge of the optimal failure model.
It should be underlined that due to the absence of data from former rockfall-related events, limitations and uncertainties of this approach have been recognized. Firstly, the limitations could occur from the specifications of the 3D point cloud dataset (i.e., points density, ground resolution, reprojection error, etc.) and of course the inexperience of the user. When it comes to the uncertainties of the proposed methodology, these are subsequent to the lack of pre-existing data which could be used to validate the results. As was shown in the Kipinas Monastery case study, no limitations regarding the quality of the constructed point cloud (consisting of about 104 million points with a ground resolution of 1.27 cm/pixel and a reprojection error of 0.157 pixels) have been identified, according to the standards set by previous related analyses [46][47][48][49][50]109]. In order for the aforementioned uncertainties to be mitigated, it was decided to identify and analyze an additional parameter, the already failed rock blocks. The monitoring of the previously fallen rock blocks and the estimation of their volumes enabled the construction of a substitute for the previously non-existent inventory of prior rockfall events in the area of interest. With the establishment of this preliminary rockfall-related dataset, the distribution and volumes of the potentially unstable rock blocks were validated and, in some cases, calibrated. As a result of the reduction in the analysis incertitude, a more accurate simulation of reality is achieved.
To conclude, it is important to highlight that this proposed rating system reflects the result of an initial attempt, and further development and validation by back-analyses is encouraged for the optimum adjustment of the weight for various contained factors. Specifically, its potential application in different areas and rock slopes, where most of the parameters may have a greater range, could help to evaluate the sensitivity of each factor in the final estimation of the hazard and risk level.

Conclusions
Rockfall hazard and risk assessments consist some of the most elaborate and challenging surveys in geoengineering. The complexity of the event's nature, combined with the importance of what is comprised in a potential failure, automatically pose their assessment as a top priority. With the introduction of remote sensing technologies, such as UAVs, some of the previous formidable deterring factors have been eased. The construction of qualitative and quantitative analyses for rockfall events became more precise and their results produce a more accurate simulation of the reality.
In general, UAVs constitute effective and efficient low-cost tools for data collection in many aspects of engineering geology [33]. Due to their ability to provide detailed information about the surface ground conditions, especially on difficult-to-access sites, many prior limitations in hazard and risk assessments can be now surpassed. With the use of digital photogrammetric techniques in the UAV-derived data, it is now possible to produce ultra-high spatial resolution datasets of hazardous regions, combined with the elimination of human risk. The further appropriate computer-based analysis of these datasets could lead to the construction of a detailed 3D visual interpretation (3D point cloud) of the investigated region that allows real-time applications and measurements of the rockfall phenomenon, which would not be attainable in that scale otherwise.
This study primarily aims to take advantage of the aforementioned capabilities that Unmanned Aerial Vehicles have to offer, in order to produce a hazard and risk assessment, with a qualitative approach, via the proposition of a new rating system. The qualitative approach was selected over a quantitative one, in order to cover a wider range of applications in areas susceptible to rockfall events with negligible to non-existent data from former rockfall-related events.
The followed methodology led to the extraction of all the necessary quantitative data that can be used for rockfall hazard and risk assessment of any natural slope, via the implementation of a new rating system. The area of application of the proposed methodology is the cultural heritage location of Kipinas Monastery in Epirus, Greece, an area characterized by the absence of pre-existing data regarding previous rockfall events, their spatial distribution along the cliff, their frequency (or the return period) of occurrence and their impact in the area of interest. Hence, this method can be suggested as a powerful approach to conduct a detailed hazard and risk assessment, in areas where negligible to no data from former rockfall-related events exist.
Initially, traditional field survey methods have been applied at the lower accessible part of the rock face, aiming to define the geometrical characteristics of the joints. Simultaneously, a photogrammetry survey via a manually flown UAV was performed. A total of 362 images were acquired from the UAV in order to cover the whole slope face. Then, the utilization of Structure from Motion (SfM) methodology enabled the production of ultra-high-resolution orthomosaic by intersecting the matched features between the overlapping offset images and interpreting them as detailed 3D point clouds. This 3D point cloud model, which consists of about 104 million points with a ground resolution of 1.27 cm/pixel, was georeferenced using the UAV sensor GPS data and a set of five artificial GCPs, and then segmented into four (4) sub-regions in order to proceed to the structural analysis. The structural analysis was performed, in a supervised, automated manner, with the open-source DSE software [66]. The result of this analysis was the extraction of the principal discontinuity sets and their main characteristics (dip, dip direction, persistence, spacing). With the geometrical characteristics of the four slopes and from the principal discontinuity sets known, a kinematic analysis was performed in order to estimate the optimal failure model for each sub-region. With the optimal failure model determined for each sub-region, a semi-automated procedure was performed in the 3D point cloud in order to monitor and quantify the magnitude (volume, mass) of both the already detached and fallen rock blocks and the prone ones for each sub-region. The final step before the implementation of the rating system was the conduction of a rockfall analysis based on trajectometry simulations. This probabilistic analysis was performed via RocFall software v.8 (RocScience Inc., Toronto, ON, Canada) for four (4) different slope profiles, one for each sub-region, in order to determine the total kinetic energy of a potential rockfall event at the moment of impact, the velocity of the rock blocks and their "bounce height".
With this methodology, all the essential data were obtained in order to provide values to all of the factors of the proposed rockfall hazard and risk rating system. In general, the proposed rating system defines four (4) major parameters categories which consist of a variety of different ranking values, in order to estimate the total risk level for each sub-region. Three (3) of these categories are related to the hazard level (geometry factors, geoengineering factors and triggering factors) and one (1) to the impact of a potential rockfall event. Every rating factor has a different weight, which is based on the significance of each individual contribution to either the hazard and/or risk level for every sub-region. Each examined factor can be evaluated between 10 and 100, as the conditions evolve from favorable to adverse, and then should multiplied by its respective predetermined weight. The total hazard score is calculated by summing the individual score of each factor, from the first three (A, B and C) major categories. At last, the total risk score is calculated by adding the fourth category (D) to the hazard score.
Lastly, in a subsequent step, in order to estimate the outcome of this rockfall hazard and risk assessment and also gain important feedback for the proposed rating system, research should continue to be performed in the cultural heritage area of Kipinas Monastery. In particular, regular monitoring of the slope, ideally per six months or annually, with a UAV is suggested in order to compare the obtained data. Through this multi-temporal change detection procedure between the different datasets, topographical changes between different epochs could be monitored and characterized for robust modeling of rockfall dynamics, as suggested by [70,106,107]. With these results, mitigation measures could be quantified and suggested for the unique needs of each sub-region. Additionally, via these outcomes, an immediate evaluation and calibration, if needed, could be performed of the proposed rating system in order be as precise and accurate as possible. This constant evaluation is necessary because the main target of this proposed rockfall hazard and risk rating system is to be considered as the transitional step between qualitative and quantitative assessments.