Rock Mass Characterization of Karstified Marbles and Evaluation of Rockfall Potential Based on Traditional and SfM-Based Methods; Case Study of Nestos, Greece

: Rockfall consists one of the most harmful geological phenomena for the man-made environment. In order to evaluate the rockfall hazard, a variety of engineering geological studies should be realized, starting from conducting a detailed ﬁeld survey and ending with simulating the trajectory of likely to fail blocks in order to evaluate the kinetic energy and the runout distance. The last decade, new technologies, i.e., remotely piloted aircraft systems (RPAS) and light detection and ranging (LiDAR) are frequently used in order to obtain and analyze the characteristics of the rock mass based on a semi-automatic or manual approach. Aiming to evaluate the rockfall hazard in the area of Nestos, Greece, we applied both traditional and structure from motion (SfM)-oriented approaches and compared the results. As an outcome, it was shown that the semi-automated approaches can accurately detect the discontinuities and deﬁne their orientation, and thus can be used in inaccessible areas. Considering the rockfall risk, it was shown that the railway line in the study area is threaten by a rockfall and consequently the construction of a rockfall netting mesh or a rock shed is recommended.


Introduction
The methodology that is usually followed for the evaluation of rockfall potential within an area is discriminated in several stages, including a detailed engineering geology field survey and analysis of the data obtained during this survey, identification of the failure type mechanism, e.g., wedge, planar slide, toppling, and the assessment of the rockfall trajectory and the evaluation of critical parameters like the kinetic energy and the bounce height of the detached rock (e.g., [1]).
The first phase of this methodology focusses on the in situ characterization of the rock mass which is based on the information provided by the field survey.Specific parameters, such as the structure of the rock mass and the quality of the discontinuities, the joint spacing, the filling material, etc., are quantified and described to assess the rock mass quality.The most applied rock mass classification systems are the rock mass rating (RMR) [2], the Q-system [3] and the geological strength index (GSI) [4,5].The first two approaches provide a semi-quantitative value of their rock mass quality through an index, while the classification based on the GSI system results to a number which, when combined with the intact rock properties, can be used for estimating the reduction in rock mass strength for different geological conditions [6].Having classified the rock mass and defined critical parameters of its structure, i.e., identification of the discontinuity sets and their orientation, a rockfall hazard analysis can follow.
A rockfall is a fragment of rock detached by sliding, toppling, or falling that falls along a vertical or sub-vertical cliff, proceeds down slope by bouncing and flying along ballistic trajectories or by rolling on talus or debris slopes [7].Although the fact that rockfall usually impacts only small areas, the damage to the infrastructure or persons directly affected may be high with serious consequences [8].Rockfalls range from small cobles to large boulders hundreds of cubic meters in size and travel at speeds ranging from few to tens of meters per second [9].Minor rockfalls affect most of the rock slopes, whereas large size ones such as cliff falls and rock avalanches, affect only great rock slopes with geological conditions favorable to instability [10].
The main factors that can trigger the detachment of a rock from a slope are rainfall, earthquakes, and human activities while the fall of a rock is determined by factors like the slope morphology [11].In addition, the joint roughness number, the joint alteration number, and the joint water reduction are important variables.The generation of rockfall is a rapid phenomenon and represents a continuous hazard in mountain areas worldwide.Rockfall presents an ongoing challenge to the safe operation of transportation infrastructure, creating hazardous conditions which can result in damage to roads and railways [12][13][14][15], as well as loss of life [16].As it has been demonstrated by [17], relatively small volumes of rocks may cause significant damage and traffic disruption, particularly in railroads.
In order to evaluate the rockfall hazard in an area and/or document the effects of a case study, remote sensing techniques were developed during the last decades such as light detection and ranging (LiDAR) (e.g., [18,19]) and photogrammetric surveys.The latter one, a photogrammetric survey, is frequently conducted using remotely piloted aircraft systems (RPAS), as well known as unmanned aerial vehicles (UAV), such as multicopters equipped with webcams, digital cameras and other sensors [20][21][22][23][24], for diverse applications such as 3D terrain analysis, slope stability, mass movement hazard, and risk management [25].RPASs provide an easy to use and low-cost remote sensing platform for rockfall studies, considered that high-resolution imagery can be acquired in short time in inaccessible areas.Using the structure from motion (SfM) image processing technique, a 3D point cloud can be created by intersecting the matched features between the overlapping, offset images [26].
The developed 3D point cloud can be proceeded for conducting a semi-automated or manually structural analysis, and for extracting information regarding the discontinuity sets, i.e., orientation and joint spacing [27][28][29][30][31][32].Having obtained these data, the evaluation of rockfall potential can be realized by performing a kinematic analysis (e.g., [33][34][35][36]).In addition, the combination of this approach with traditional engineering geological field surveys results to a rock mass classification and to slope mass rating assessment [37,38].
This study focuses on the area of Nestos, Greece, at a site where the national railway, connecting Greece and Turkey, passes in a short distance (<10 m) from a vertical cliff of marble.The marble unit is the dominant formation in the research area, consisted by massive white marbles with dolomite layers.These marbles appear to be karstic on the surface in many places.The rock slope studied for the purposes of this research is a white to light grey colored poorly bedded massive marbles of 35 m high (Figure 1).The goal of this study is twofold: (i) to perform a rock mass classification in order to evaluate the rockfall potential by assigning slope mass rating (SMR) index values, and (ii) to simulate the trajectories of likely to fail blocks in order to examine their run-out distance and to assess the rockfall hazard on the railway line.In order to achieve this goal, a RPAS survey conducted, and a 3D point cloud was generated (Figure 1).The lower part of the rock slope, i.e., site 1 was investigated both by traditional field survey and point-cloud approaches while the higher part (site 2) was analyzed based on the generated 3D model.The procedure that was followed in this study is shown in the following flowchart (Figure 2).

Engineering Geological Field Survey
The first step for the characterization of a rock mass is to perform a detailed engineering geological field survey.More specifically, for the purposes of this study, we used a geological compass in order to measure the orientation of the discontinuities, in terms of dip and dip direction, while a Schmidt hammer was used for estimating the uniaxial compressive strength (UCS) of the intact rock and the joint compressive strength (JCS).Furthermore, the roughness of the joints was assessed based on a profilometer while the aperture and joint spacing were measured using a

Engineering Geological Field Survey
The first step for the characterization of a rock mass is to perform a detailed engineering geological field survey.More specifically, for the purposes of this study, we used a geological compass in order to measure the orientation of the discontinuities, in terms of dip and dip direction, while a Schmidt hammer was used for estimating the uniaxial compressive strength (UCS) of the intact rock and the joint compressive strength (JCS).Furthermore, the roughness of the joints was assessed based on a profilometer while the aperture and joint spacing were measured using a surveyor's tape.Finally, the filling material and the weathering conditions at the area were also assessed.
Data acquisition of discontinuities is usually based either on scanline or window sampling.In our case, we "applied" the former one due to the steep morphology of the area.More specifically, only the lower part of the rock face was accessible (site 1, Figure 1), and the acquired information is related to a field survey conducted along a scanline of 4 m length (Figure 3).Even though the conducted structural analysis has limitations regarding the representation of the whole rock face, it is highlighted that an engineering geological field survey is mandatory for the rock mass characterization since crucial parameters, e.g., joint roughness and weathering, can be accurately assessed based only on this traditional approach.
Geosciences 2020, 10, x FOR PEER REVIEW 5 of 19 surveyor's tape.Finally, the filling material and the weathering conditions at the area were also assessed.Data acquisition of discontinuities is usually based either on scanline or window sampling.In our case, we "applied" the former one due to the steep morphology of the area.More specifically, only the lower part of the rock face was accessible (site 1, Figure 1), and the acquired information is related to a field survey conducted along a scanline of 4 m length (Figure 3).Even though the conducted structural analysis has limitations regarding the representation of the whole rock face, it is highlighted that an engineering geological field survey is mandatory for the rock mass characterization since crucial parameters, e.g., joint roughness and weathering, can be accurately assessed based only on this traditional approach.Having defined the orientation of joints on the field, the dip and dip direction of the dominant discontinuities were defined based on the software Dips v.7 (Rocscience, Toronto, Canada) and resulted to three main discontinuity sets ( ,  , and  ).Furthermore, 15 Schmidt hammer tests were carried out, yielding a mean value of 44 rebounds.
Moreover, having collected all the previously information on the field, it was feasible to proceed to the rock mass characterization according to the RMR system [2], where RMR = RMRbasic + R6 and RMRbasic is equal to R1 + R2 + R3 + R4 + R5.More specifically, the following five parameters have been evaluated: R1-uniaxial compressive strength of rock material; R2-rock quality designation (RQD); R3-spacing of discontinuities; R4-discontinuity conditions; and R5-groundwater conditions.The sixth parameter concerns the effect of discontinuity orientation and is not considered for the estimation of RMRbasic.The values of RMRbasic range from 0 to 100, and the rock mass is classified into five classes: 0-20 (Very poor), 21-40 (Poor), 41-60 (Fair), 61-80 (Good), and 81-100 (Very good) [2].
Taking into account the data provided by the conducted field survey, the uniaxial compressive strength of intact rock is equal to 101 MPa while the RQD is equal to 79% (R2 = 17) as it was indirectly estimated following the Equations ( 1) and (2), suggested by [40], relating the total number of joints per m 3 (Jv) with the RQD.In our case, the Jv is equal to where S1, S2, and S3 is the joint spacing for the discontinuity sets  1 ,  2 , and  3 , respectively, and RQD is equal to Having defined the orientation of joints on the field, the dip and dip direction of the dominant discontinuities were defined based on the software Dips v.7 (Rocscience, Toronto, Canada) and resulted to three main discontinuity sets (J 1 , J 2 , and J 3 ).Furthermore, 15 Schmidt hammer tests were carried out, yielding a mean value of 44 rebounds.
Moreover, having collected all the previously information on the field, it was feasible to proceed to the rock mass characterization according to the RMR system [2], where RMR = RMR basic + R6 and RMR basic is equal to R1 + R2 + R3 + R4 + R5.More specifically, the following five parameters have been evaluated: R1-uniaxial compressive strength of rock material; R2-rock quality designation (RQD); R3-spacing of discontinuities; R4-discontinuity conditions; and R5-groundwater conditions.The sixth parameter concerns the effect of discontinuity orientation and is not considered for the estimation of RMR basic .The values of RMR basic range from 0 to 100, and the rock mass is classified into five classes: 0-20 (Very poor), 21-40 (Poor), 41-60 (Fair), 61-80 (Good), and 81-100 (Very good) [2].
Taking into account the data provided by the conducted field survey, the uniaxial compressive strength of intact rock is equal to 101 MPa while the RQD is equal to 79% (R2 = 17) as it was indirectly estimated following the Equations ( 1) and (2), suggested by [40], relating the total number of joints per m 3 (J v ) with the RQD.In our case, the Jv is equal to where S1, S2, and S3 is the joint spacing for the discontinuity sets J 1 , J 2 , and J 3 , respectively, and RQD is equal to The mean values of the spacing of the three discontinuity sets are measured as J 1 = 0.3 m, J 2 = 0.18 m, and J 3 = 0.3 m, and accordingly the values of 10, 8, and 10 were assigned for the parameter R3, respectively.The surface of discontinuities is characterized as rough with unweathered walls and a separation less than 0.1 mm (R4 = 28), while the conditions regarding the presence of groundwater in joints was reported as dry; R5 = 15.Thus, an estimated range of rock mass rating (RMR) of 80-82 with a mean of 81 was assigned, and consequently the rock mass is classified as very good.In the following Table 1 are listed the assigned values of the five parameters required for the estimation of RMR basic per discontinuity set.

Development of the 3D Model
Considering the steep morphology of the rock slope, it was decided to conduct a RPAS survey in September 2019 in order to reconstruct an accurate and detailed 3D model of the studied area.RPAS imagery was collected using a DJI Phantom 4 Pro V2.0 equipped with a 1-inch 20-megapixel sensor and a FC6310S camera of 8.8 mm focal length and a manually adjustable aperture from F2. 8 to F11.A total of 74 images were acquired, covering the rock face with slightly oblique views.The images obtained from the RPAS survey were processed using the AgisoftTM Metashape software v. 1.5.2.(Agisoft LLC, St. Petersburg, Russia) resulting in a dense point cloud dataset of 67 million points with ground resolution 4.36 mm/pix and a reprojection error of 0.838 pix.Further processing of the point cloud data enabled the production of a detailed DEM with pixel dimensions 8.41 mm.The RPAS flown manually due to the complex geomorphology.The point cloud was georeferenced using the RPAS sensor GPS data and a set of 3 artificial ground control points GCP.Considering the steep geometry of the slope and the relevant hazard to access higher parts, the GCP have been placed in accessible lower parts, with a RMSE of 20.4 cm (XY) and 10.7 cm (Z).
Having georeferenced the 3D point cloud, an additional validation was performed by comparing the orientation of specific joints measured by geological compass with the one extracted by the compass plugin on CloudCompare™.It should be pointed out that the latter procedure, i.e., CloudCompare™ is a "manual", non-automatized one.As it is shown in Figure 4, the difference between the extracted information based on point cloud-based procedure and the geological compass is less than 10%.
accessible lower parts, with a RMSE of 20.4 cm (XY) and 10.7 cm (Z).
Having georeferenced the 3D point cloud, an additional validation was performed by comparing the orientation of specific joints measured by geological compass with the one extracted by the compass plugin on CloudCompare™.It should be pointed out that the latter procedure, i.e., CloudCompare™ is a "manual", non-automatized one.As it is shown in Figure 4, the difference between the extracted information based on point cloud-based procedure and the geological compass is less than 10%.

Structural Analysis Based on Point Cloud-Oriented Approaches
In order to detect the discontinuity sets and measured their orientation at the studied area, the generated point cloud was analyzed by applying recently developed SfM-based methodologies.Supervised analysis using the open source software Discontinuity Set Extractor (DSE) (developed by [30]), and manually oriented (CloudCompare™) approaches were applied in the developed 3D model.Considering the former approach, a supervised classification of discontinuity sets was performed to the obtained 3D point cloud by applying the methodology proposed by [30].This approach calculates for every point the normal vector of a set of "knn" neighbors and its corresponding pole in a stereonet plot.As the normal is calculated using the PCA algorithm, a test is performed using the concept of tolerance η: the value of the third eigenvalue over the overall sum of eigenvalues.
According to [30], if η is greater than a threshold, the set of points could be considered as non-planar and discarded, but if not, it could be considered as a candidate of a planar area.Then, using the Kernel Density Estimation (3) [41], the density function is estimated and finding the peaks of the most represented orientations, the user determines the number of principal planes and their orientations.The method assigns a set to every point if the angle between the point's normal and the principal plane is less than a threshold value (stablished by the user).Finally, the points members of a set are extracted, and a spatial cluster analysis is performed using the DBSCAN method [42].Each cluster of points is member or a plane, so planes' equations are calculated for every outcrop.This result enables subsequent studies of normal spacing [30] and persistence [31].In particular, the following parameters were employed on DSE; number of nearest neighbors (knn) = 30, η = 0.2, while the minimum angle between the normal vector of a discontinuity set and the normal vector of the point was 30 (cone = 30).Moreover, once planes are identified and located on the slope face, then the normal spacing can be automatically estimated.
Furthermore, the manual procedure was followed by using the Compass plugin of the freeware CloudCompare™ software, aiming to measure discontinuity orientations directly on the point cloud.Having extracted the orientation of the joints, the developed stereonet was analyzed on the Dips v7 software developed by Rocscience™.
In addition, based on the information provided by the point cloud, i.e., joint spacing, it was feasible to calculate both the mean block volume at the rock slope based on Equation ( 6), while the volume of specific blocks that are likely to be detached form the rock face was manually measured on 3D model.The results obtained from the application of the precedent approaches at site 1 and 2 are described in the following subsections.

Structural Analysis at Site 1
The dimensions of the 3D studied area of site 1 is approximately 2 m length and 2 m height covering an area of 4 m 2 .The number of points of the original point cloud is 94,502, while the relevant number of points of the classified point cloud is 70,109.The discontinuity sets at site 1, as they were semi-automatically extracted based on DSE (Figure 5), are three (3) and their characteristics are listed in the following Table 2.
Geosciences 2020, 10, x FOR PEER REVIEW 8 of 19 Where % is the number of assigned points to a DS over the total number of points.Furthermore, the discontinuity spacing was assessed based on the procedure that can be followed on DSE.More specifically, having extracted the discontinuity sets (DS), [30] developed a methodology where every single point is labeled with its corresponding DS and every exposed planar surface is calculated.Then this method computes for every DS the normal spacing between an exposed plane and its nearest point member of the same DS.As a result, the normal spacing is calculated as the mean value of all the normal spacings for each DS [30].In this case study, the spacing for the discontinuity set  is 0.22 m, for  is 0.23 m and for  is 0.19 m (Figure 6), while the relevant values resulted from the field survey are 0.3 m, 0.18 m, and 0.3 m, respectively.Based on geological judgement and having estimated on the field the linear persistence of  ,  , and  as less than 1 m (very low persistence according to the classification suggested by [32]), the DS have been characterized as non-persistent.
Thus, considering the point cloud-based values of spacing, the indirectly measured value of RQD is 69% (R2 = 13) based on the proposed equation by [40], while the relevant parameter R3 of RMR classification is assessed as 10 and 8. Adopting that the other parameters required for the  Where % is the number of assigned points to a DS over the total number of points.
In addition, the same area has been investigated regarding the joint orientations based on Compass Plugin on CloudCompare™ software.More specifically, 44 measurements of discontinuities orientation were performed, and analyzed in stereographic projection on Dips v7.Afterwards, the determination of the representative orientation (dip, dip direction) of the discontinuity sets was extracted (Figure 5) and listed in Table 2. Furthermore, the discontinuity spacing was assessed based on the procedure that can be followed on DSE.More specifically, having extracted the discontinuity sets (DS), [30] developed a methodology where every single point is labeled with its corresponding DS and every exposed planar surface is calculated.Then this method computes for every DS the normal spacing between an exposed plane and its nearest point member of the same DS.As a result, the normal spacing is calculated as the mean value of all the normal spacings for each DS [30].In this case study, the spacing for the discontinuity set J 1 is 0.22 m, for J 2 is 0.23 m and for J 3 is 0.19 m (Figure 6), while the relevant values resulted from the field survey are 0.3 m, 0.18 m, and 0.3 m, respectively.Based on geological judgement and having estimated on the field the linear persistence of J 1 , J 2 , and J 3 as less than 1 m (very low persistence according to the classification suggested by [32]), the DS have been characterized as non-persistent.

Comparing Traditional and 3D-Based Characterization of Rock Mass at Site 1
The outcome arisen from the comparison of the information regarding the discontinuity sets obtained by the DSE and CloudCompare™, as well as from the field survey based on geological compass is shown in Figure 4 and listed in Table 2.It can be primarily stated that the measurements of discontinuity set  are almost in agreement at all the approaches.The  discontinuity set was also identified in the three systems, as a high angle joint.However, it should be pointed out that although the fact that the measured strike of  is similar among the approaches, the main dip direction, as it was assessed based on the data provided by the geological compass, is opposite to the ones extracted by the point cloud-based approaches; as it is shown in Figure 5, poles of similar dip direction have also been identified both on DSE and CloudCompare™.This difference could be preliminary justified by the fact that by analyzing the developed point cloud, either based on semiautomated (DSE) or manually approaches (CloudCompare™), it is feasible to obtain a larger set of data comparing to the one provided by the traditional field survey.Therefore, it is believed that the point cloud-based approaches can assess the orientation of the joints in a more representative way regarding the rock mass.The orientation of the extracted  DS is almost in agreement between the DSE and CloudCompare™ approaches.The relevant joint system ( ) that was assessed based on the analysis of the data provided by the geological compass, has similar dip direction with the ones extracted by the point-cloud based approaches but it is considered as a lower angle joint.

Structural Analysis at Site 2
The studied area of site 2 has maximum length and width 8 m and 13 m, respectively.It is considered as an area free of human activities and a source zone of future detachments of rock blocks.Considering the high elevation of the site, the orientation of the discontinuity sets was retrieved solely based on the approaches, i.e., DSE and CloudCompare™ that were applied on the developed point cloud.
More specifically, concerning the former approach (DSE), the number of points of the original point cloud is 113,797, while the relevant number of points of the classified point cloud is 57,454.The discontinuity sets at site 2, as they were semi-automatically extracted based on DSE (Figure 7), are three and their characteristics are provided on the following Table 3.
In addition, the same area has been investigated regarding the joint orientations based on Compass Plugin on CloudCompare™ software.More specifically, 58 measurements of discontinuities orientation were performed, which have been further analyzed by means of a density analysis of attitudes in stereographic projection on Dips v7 of Rocscience (Figure 7).Then, the determination of the representative orientation (dip and dip direction) of the discontinuity sets was extracted and listed in Table 3.Thus, considering the point cloud-based values of spacing, the indirectly measured value of RQD is 69% (R2 = 13) based on the proposed equation by [40], while the relevant parameter R3 of RMR classification is assessed as 10 and 8. Adopting that the other parameters required for the evaluation of RMR basic remain the same, the average value of RMR basic is 78 and the rock mass is classified as good.

Comparing Traditional and 3D-Based Characterization of Rock Mass at Site 1
The outcome arisen from the comparison of the information regarding the discontinuity sets obtained by the DSE and CloudCompare™, as well as from the field survey based on geological compass is shown in Figure 4 and listed in Table 2.It can be primarily stated that the measurements of discontinuity set J 1 are almost in agreement at all the approaches.The J 2 discontinuity set was also identified in the three systems, as a high angle joint.However, it should be pointed out that although the fact that the measured strike of J 2 is similar among the approaches, the main dip direction, as it was assessed based on the data provided by the geological compass, is opposite to the ones extracted by the point cloud-based approaches; as it is shown in Figure 5, poles of similar dip direction have also been identified both on DSE and CloudCompare™.This difference could be preliminary justified by the fact that by analyzing the developed point cloud, either based on semi-automated (DSE) or manually approaches (CloudCompare™), it is feasible to obtain a larger set of data comparing to the one provided by the traditional field survey.Therefore, it is believed that the point cloud-based approaches can assess the orientation of the joints in a more representative way regarding the rock mass.The orientation of the extracted J 3 DS is almost in agreement between the DSE and CloudCompare™ approaches.The relevant joint system (J 3 ) that was assessed based on the analysis of the data provided by the geological compass, has similar dip direction with the ones extracted by the point-cloud based approaches but it is considered as a lower angle joint.

Structural Analysis at Site 2
The studied area of site 2 has maximum length and width 8 m and 13 m, respectively.It is considered as an area free of human activities and a source zone of future detachments of rock blocks.Considering the high elevation of the site, the orientation of the discontinuity sets was retrieved solely based on the approaches, i.e., DSE and CloudCompare™ that were applied on the developed point cloud.
More specifically, concerning the former approach (DSE), the number of points of the original point cloud is 113,797, while the relevant number of points of the classified point cloud is 57,454.The discontinuity sets at site 2, as they were semi-automatically extracted based on DSE (Figure 7), are three and their characteristics are provided on the following Table 3.
Geosciences 2020, 10, x FOR PEER REVIEW 10 of 19 Where % is the number of assigned points to a DS over the total number of points.Having obtained the orientations of the discontinuity sets (DS) by the DSE and CloudCompare™, it can be seen in Figure 7 that the discontinuity sets that formed the rock mass are three ( ,  , and  ) and they are almost identically matching regarding the orientation (Dip and Dip direction).For this area that has not been affected by blast techniques and excavation, we could state  Where % is the number of assigned points to a DS over the total number of points.
In addition, the same area has been investigated regarding the joint orientations based on Compass Plugin on CloudCompare™ software.More specifically, 58 measurements of discontinuities orientation were performed, which have been further analyzed by means of a density analysis of attitudes in stereographic projection on Dips v7 of Rocscience (Figure 7).Then, the determination of the representative orientation (dip and dip direction) of the discontinuity sets was extracted and listed in Table 3.
Having obtained the orientations of the discontinuity sets (DS) by the DSE and CloudCompare™, it can be seen in Figure 7 that the discontinuity sets that formed the rock mass are three (J 1 , J 2 , and J 3 ) and they are almost identically matching regarding the orientation (Dip and Dip direction).For this area that has not been affected by blast techniques and excavation, we could state that these two approaches match.Approximately 43% of the points fall on surfaces belonging to J 1 while the second prevailing joint system is J 2 (36%).The gently-dipping discrete J 3 joint set is also present (21%).
Having defined the DS, we proceed to the estimation of the normal spacing following the procedure proposed by [30] by considering the cluster analysis on DSE software.Having tested several combinations, at the relevant field on DSE, we finally estimated the number of clusters by applying a value 200 as the minimum of number of points per cluster for J 1 , J 2 , and J 3 .As an outcome, the density function of the measured normal spacing of the DS is computed (Figure 8).that these two approaches match.Approximately 43% of the points fall on surfaces belonging to  while the second prevailing joint system is  (36%).The gently-dipping discrete  joint set is also present (21%).
Having defined the DS, we proceed to the estimation of the normal spacing following the procedure proposed by [30] by considering the cluster analysis on DSE software.Having tested several combinations, at the relevant field on DSE, we finally estimated the number of clusters by applying a value 200 as the minimum of number of points per cluster for  ,  , and  .As an outcome, the density function of the measured normal spacing of the DS is computed (Figure 8).
Based on geological judgment and measurements on the 3D model, the discontinuity sets  and  are considered as a full persistent (linear persistence > 1 m) with a spacing of 0.59 m and 0.5 m, respectively.The third discontinuity set ( ) is considered as non-persistent ((linear persistence < 1 m) and the relevant number of spacing that was considered for the calculation of Jv is 0.82 m.Thus, taking into account the proposed by [40,43] correlation between RQD and the "volumetric joint count" (number of joints per cubic meter), when drill cores are not available, we estimated that Jv is equal to five and the RQD of the rock mass is computed as 98%.
Considering that the block size is a volumetric expression of joint density [44] and that in our case the extracted discontinuity sets occur approximately at right angles to each other, the mean volume Vblocks of the blocks that could be detached from the rock slope was estimated based on the Equation ( 6) proposed by [43].In particular, the mean block volume at the rock slope is computed as where  ,  , and  is the estimated normal spacing of discontinuity sets  ,  , and  , respectively.Adopting that the parameters regarding the joint roughness, weathering and filling are similar with the ones obtained by the conventional field survey on the lower part of the rock face, the mean value of RMRbasic is computed as 82 (83 for  and 81 for  and  ), and the rock mass is classified as very good.The values of RMRbasic, assigned at the lower (site 1) and the higher (site 2) part of the rock slope, are similar and some minor changes are due to the differences in joint spacing.

SMR Assessment
For the purposes of this study, the slope mass rating (SMR) classification was assessed using the open source software SMRTool [29,37].This software guides the user through the calculation of the SMR index introducing the orientations of the slope and discontinuity set (i.e., dip and dip direction), the RMRbasic value and the excavation method.The software calculates the failure mode, the auxiliary angles, the adjustment factors and the RMR index value [37].In addition to the calculations, it shows a sketch of the orientations of the slope and the discontinuity set to aid the user to understand the Based on geological judgment and measurements on the 3D model, the discontinuity sets J 1 and J 2 are considered as a full persistent (linear persistence > 1 m) with a spacing of 0.59 m and 0.5 m, respectively.The third discontinuity set (J 3 ) is considered as non-persistent ((linear persistence < 1 m) and the relevant number of spacing that was considered for the calculation of J v is 0.82 m.Thus, taking into account the proposed by [40,43] correlation between RQD and the "volumetric joint count" (number of joints per cubic meter), when drill cores are not available, we estimated that J v is equal to five and the RQD of the rock mass is computed as 98%.
Considering that the block size is a volumetric expression of joint density [44] and that in our case the extracted discontinuity sets occur approximately at right angles to each other, the mean volume V blocks of the blocks that could be detached from the rock slope was estimated based on the Equation (6) proposed by [43].In particular, the mean block volume at the rock slope is computed as where S 1 , S 2 , and S 3 is the estimated normal spacing of discontinuity sets J 1 , J 2 , and J 3 , respectively.Adopting that the parameters regarding the joint roughness, weathering and filling are similar with the ones obtained by the conventional field survey on the lower part of the rock face, the mean value of RMR basic is computed as 82 (83 for J 1 and 81 for J 2 and J 3 ), and the rock mass is classified as very good.The values of RMR basic , assigned at the lower (site 1) and the higher (site 2) part of the rock slope, are similar and some minor changes are due to the differences in joint spacing.

SMR Assessment
For the purposes of this study, the slope mass rating (SMR) classification was assessed using the open source software SMRTool [29,37].This software guides the user through the calculation of the SMR index introducing the orientations of the slope and discontinuity set (i.e., dip and dip direction), the RMR basic value and the excavation method.The software calculates the failure mode, the auxiliary angles, the adjustment factors and the RMR index value [37].In addition to the calculations, it shows a sketch of the orientations of the slope and the discontinuity set to aid the user to understand the calculation.The software is programmed in MATLAB and available on the Internet (https://personal.ua.es/en/ariquelme/smrtool.html).The SMR provides a preliminary assessment of the slope stability and it is obtained from RMR by subtracting a factorial adjustment factor depending on the joint-slope relationship and adding a factor depending on the method of excavation.As a classification system, it was recommended by [45] as an additional material to RMR concept, which is suitable for slopes, and has been worldwide used during last thirty years as a geomechanics classification for rating rocky slopes [46].
Through the widely application of SMR, several issues have been identified according to [46].In particular, the SMR classification is slightly conservative while the slope stability is highly influenced by the excavation method justifying its inclusion.In addition, the failure modes derived from SMR are validated in practice and the slope height is not taken into consideration for the classification of the rock mass based on SMR.
The proposed five different rock classes by [45] are class V (SMR < 20) described as very bad where the type of failures that should be expected are big planar or soil-like, class IV (20 < SMR < 40) described as bad with planar and big wedges, class III (40 < SMR <60) described as normal that requires systematic support for joints and wedges, class II (60 < SMR < 80) described as good requiring occasional support at some blocks, and class I (80 < SMR < 100) for completely stable slopes where none support is required.
In this study, having employed on SMRtool the requested information, the SMR per discontinuity set and per likely to developed wedges is calculated and the possible slope failure type is displayed based on both approaches proposed by [47,48].It is additionally pointed out that the SMRtool provides two graphics where the dip direction and the dip of the slope and discontinuity set are plotted, helping users to understand the expected type of failure [49].The proposed slope mass rating (SMR) is obtained from RMR by using the following formula: where F1 depends on parallelism between joints and slope face strikes and particularly between the angular relation between the dip direction of the considered joint (αj) and slope dip direction (αs).F2 refers to joint dip angle (βj) for the planar mode of failure.F3 reflects the relationship between the dip of slope (βs) and the dip of joints (βj).F4 is the adjustment factor for the method of excavation used for the studied slope; for natural slopes is equal to 15, for presplitting F4 = 10, for smooth blasting F4 = 8, and for normal blasting or mechanical excavation F4 = 0.
The orientation of the natural slope is derived during fieldwork by using a geological compass.As it is shown in Figure 9, the slope is divided in two sectors, named S1 and S2, regarding its direction.In particular, the dip direction of the sector 1 is defined as 160 • (S1) and 180 • for the sector 2 (S2) with a dip of 85; the direction of the latter slope plane (S2) is parallel to the one of the railway networks.Considering that the sector S1 is a natural slope while the sector 2 was excavated based on blasting and mechanical methods, the relevant employed values of F4 are 15 and 0, respectively.For the natural slope (S1), the SMR has been calculated using SfM data sets from the site 2 as they have been analyzed by means of DSE software.As an outcome, the SMR values that were calculated for the three delineated discontinuity sets, based on [47] approach, are 94, 92, and 87 for J 1 , J 2 , and J 3 , respectively, while for the wedges W12, W13, and W23 the relevant SMR values are 86, 96, and 87, respectively.The relevant values based on the [48] approach are 92, 89, and 84 for the discontinuity sets J 1 , J 2 , and J 3 , respectively, and for the developed wedges are 86, 95, and 86, respectively.Based on these results, the discontinuity set J 3 and the wedges formed by J 2 -J 3 and J 1 -J 2 are related to the lower rating of the slope (Table 4).
For the slope that had been excavated (85/180) for the construction of the railway, the relevant values of SMR that were calculated for the three delineated discontinuity sets are 79, 71, and 72 for J 1 , J 2 , and J 3 , respectively, while for the wedges W 12 , W 13 , and W 23 the relevant SMR values are 77, 81, and 72, respectively, following the procedure suggested by [47].Based on these results, the slope is considered as stable while some blocks may be detached from the slope; the failure probability is approximated as 20% according to [50].The discontinuity set J 2 and the wedge formed by J 2 and J 3 are related to the lower rating of the slope (Table 5).

Evaluation of Rockfall Hazard
Following the visual inspection during the field survey, a 3D point cloud analysis on CloudCompare™ took place aiming to determine the location and the geometry of potentially likely to fail blocks.As it is shown in Figure 10, three source areas, i.e., blocks, have been clearly identified and the shape and size (volume) of each area were measured based on the tools provided by the CloudCompare™.Two blocks (points a and b in Figure 10) are considered as cases that are likely to fail while the third point (c in Figure 10) is considered as the source area of a block that has been already detached and fallen.
Geosciences 2020, 10, x FOR PEER REVIEW 14 of 19 to fail blocks.As it is shown in Figure 10, three source areas, i.e., blocks, have been clearly identified and the shape and size (volume) of each area were measured based on the tools provided by the CloudCompare™.Two blocks (points a and b in Figure 10) are considered as cases that are likely to fail while the third point (c in Figure 10) is considered as the source area of a block that has been already detached and fallen.
The volume of the block shown in the following Figure 10b is manually measured as V = 0.26m 3 based on the tools provided on CloudCompare™, a value that is in agreement with the estimated mean volume V of the blocks (V = 0.24m 3 ), calculated based on the extracted information regarding the joint spacing.Considering the third point (point c), the maximum volume of the detached block that can indirectly be estimated from the 3D model is equal to V = 1.12m 3 , a value that is four times larger than the estimated mean block volume.This may be justified in the case that this block was not detached as a unique boulder from the slope face, but it had been previously fragmented on two or more smaller blocks.Nevertheless, this case indicates the potential for larger blocks to being formed, detached, and failed on the railway.The next step for the evaluation of rockfall hazard is to simulate the most critical profiles as they have been identified with the 3D approach [33].The determination of these fall tracks in this study was achieved by considering the geometry of the slope as it was extracted based on the information provided by the developed point cloud.The goal of this run-out analysis is to evaluate critical parameters of the rockfall such as kinetic energy and bounce height, and to evaluate the run-out distance in order to qualitatively assess the possibility of a rockfall to reach the railway, which is located at the base of the rock face.In our case, considering the estimated values of kinetic energy of rockfalls that were simulated at this sector, the fall tracks that showed the higher values are the ones The volume of the block shown in the following Figure 10b is manually measured as V = 0.26 m 3 based on the tools provided on CloudCompare™, a value that is in agreement with the estimated mean volume V of the blocks (V = 0.24 m 3 ), calculated based on the extracted information regarding the joint spacing.Considering the third point (point c), the maximum volume of the detached block that can indirectly be estimated from the 3D model is equal to V = 1.12 m 3 , a value that is four times larger than the estimated mean block volume.This may be justified in the case that this block was not detached as a unique boulder from the slope face, but it had been previously fragmented on two or more smaller blocks.Nevertheless, this case indicates the potential for larger blocks to being formed, detached, and failed on the railway.
The next step for the evaluation of rockfall hazard is to simulate the most critical profiles as they have been identified with the 3D approach [33].The determination of these fall tracks in this study was achieved by considering the geometry of the slope as it was extracted based on the information provided by the developed point cloud.The goal of this run-out analysis is to evaluate critical parameters of the rockfall such as kinetic energy and bounce height, and to evaluate the run-out distance in order to qualitatively assess the possibility of a rockfall to reach the railway, which is located at the base of the rock face.In our case, considering the estimated values of kinetic energy of rockfalls that were simulated at this sector, the fall tracks that showed the higher values are the ones related to the blocks a and b that are located at higher elevation.
The 2D fall tracks of the rockfall was simulated using the Rocfall v.7 software, which is a robust, easy to use computer program that is available from Rocscience.This software performs a probabilistic simulation of rockfalls and can be used to design remedial measures and test their effectiveness [51].In general, a detached block descends very rapidly through the air by falling, bouncing, and rolling [52], while its trajectory mainly depends on the geometry of the slope.Considering the type of rockfalls in this study, it is shown in Figure 11 that both cases represent a free fall.The initial points in the models have been defined as a single point in Rocfall v.7 software (Rocscience, Toront, ON, Canada) based on the information provided by the point cloud.
In particular, the profile_1 concerns the fall track of the block detached by point a (Figure 11a).The result of this simulation indicates that this failure cannot threaten the railway and consequently the rockfall hazard can be neglected.This is because after the rebound on the ground, the direction of the boulder is towards the base of the rock slope and not towards the railway.On the other hand, the profile_2 which represents the trajectory of the block detached from point b, indicates that the rockfall will end up at the railway after the rebound on the ground (Figure 11b,c).
The maximum total kinetic energy of the rockfall (profile_2) at the location of the impact with the ground is 134 kJ while the maximum height of the rockfall after the rebound occurs in a distance of 5 m from the base of the rock slope and was estimated as 3.6 m.At this location, the maximum value of the kinetic energy was estimated as 9.3 kJ, while the probability of a boulder to exceed the 5 kJ of the total kinetic energy is approximately 15% (8 of 50 blocks used as a sample).The relevant value of a boulder to reach the height of 2 m is less than 30% (14 of 50 blocks).
In particular, the profile_1 concerns the fall track of the block detached by point a (Figure 11a).The result of this simulation indicates that this failure cannot threaten the railway and consequently the rockfall hazard can be neglected.This is because after the rebound on the ground, the direction of the boulder is towards the base of the rock slope and not towards the railway.On the other hand, the profile_2 which represents the trajectory of the block detached from point b, indicates that the rockfall will end up at the railway after the rebound on the ground (Figure 11b,c).The maximum total kinetic energy of the rockfall (profile_2) at the location of the impact with the ground is 134 kJ while the maximum height of the rockfall after the rebound occurs in a distance of 5 m from the base of the rock slope and was estimated as 3.6 m.At this location, the maximum value of the kinetic energy was estimated as 9.3 kJ, while the probability of a boulder to exceed the 5 kJ of the total kinetic energy is approximately 15% (8 of 50 blocks used as a sample).The relevant value of a boulder to reach the height of 2 m is less than 30% (14 of 50 blocks).

Discussion and Conclusions
This study primarily aims to characterize the rock mass and define the orientation of the discontinuity sets using conventional methods and by applying recently developed SfM-based methodologies on a rock slope located at the area of Nestos, Greece.In particular, the latter ones are focused on the development of a 3D model of the investigated area and the extraction of information regarding the orientation and spacing of discontinuities based on semi-automated and manually procedures.

Discussion and Conclusions
This study primarily aims to characterize the rock mass and define the orientation of the discontinuity sets using conventional methods and by applying recently developed SfM-based methodologies on a rock slope located at the area of Nestos, Greece.In particular, the latter ones are focused on the development of a 3D model of the investigated area and the extraction of information regarding the orientation and spacing of discontinuities based on semi-automated and manually procedures.
Initially, traditional field survey methods have been applied at the lower part of the rock face aiming to define the geometrical characteristics of the joints.It is pointed out that this site is affected by the blasting techniques applied during the construction of the railway network.For the same area, this information was additionally extracted based on open source software, i.e., Discontinuity Set Extractor (DSE) [30] and CloudCompare™.The DSE-based approach aids the discontinuity sets identification and gathers points in planar clusters.Once planes are identified and located on a slope face, then the normal spacing can be automatically estimated.On the other hand, the procedure required by CloudCompare™ for extracting the orientation of joints, is totally subjective and can be influenced by the quality of the generated point cloud and the experience of the researcher.
This study shows that on heavily fragmented rock mass, i.e., site 1, the CloudCompare™ approach could not be considered as appropriate due to the small size facets and is not recommended for the identification of joints' orientation.In addition, at site 1, a discrepancy exists between the estimated spacing of joints resulted from the traditional and the RPAS-based methods.This outcome agrees with the conclusion of [30] that more clusters can be detected based on the point cloud-based approach, resulting to a smaller spacing between the discontinuity sets, compared to the ones measured during traditional investigation.Regarding the classification of the rock mass, the DSE-based characterization slightly differs from the one performed based on the information provided by the field survey.This difference is due to the fact that the structural analysis performed by DSE is more detailed in comparison to the coarser analysis that is conducted during field survey in accessible zones.
Therefore, at reachable heavily fragmented rock masses, it is suggested to initially perform a detailed field survey that can be used as core data and afterwards conduct a structural analysis based on the developed 3D point-cloud, using open source software, i.e., Discontinuity Set Extractor, DSE, in order to define the orientation of the discontinuity sets (DS) and to estimate the joint spacing.
Considering the second site, it was concluded that both DSE and CloudCompare™ extracted similar information regarding the orientation of discontinuity sets.This probably happens because this area is free of human activities, and consequently not so fragmented as the one at the lower part of the slope.Thus, at the areas that are not reachable due to steep morphology and where conventional methods cannot be applied, the point-cloud approaches should be used for the rock mass classification and the identification of the orientation of discontinuity sets.At the zones where planar surfaces can clearly be defined, CloudCompare™ could be accurately define the dominant discontinuity sets and can be additionally used for validating the orientation of joints.However, comparing to DSE, the extraction of this information is considered as time consuming.
Combining the information collected during traditional field survey and the data extracted from the generated 3D model, rock mass characterization was realized based on classification developed by [2], i.e., RMR and afterwards the slope mass rating (SMR) was assigned.
The secondary goal of this research was the assessment of rockfall hazard and the evaluation of critical parameters of the fall.In order to achieve this, we used the developed 3D model of the slope face for detecting and estimating the volume of potentially unstable blocks.This procedure can be performed using the relevant tool on CloudCompare™ while the accuracy of the obtained block dimensions strongly depends on the quality of the point cloud.
In this case study, two source areas were identified and the characteristics of the likely to occur rockfalls were evaluated.Considered the results of the simulation, it is concluded that only one of the two examined fall tracks can threat the railway line.The proximity of the rail network to the rock slope and the type of rockfall, i.e., free fall, indicate that construction of a rockfall barrier or a rock shed are the recommended protection measures.However, scaling is suggested as the primary mitigation measure that should be realized in advance in order to minimize the risk.

Figure 1 .
Figure 1.(a) Black circle indicates the location of the study area at Northeastern Greece; (b) simplified geological map (1:50,000 scale) of the wider area modified by [39], where the study area is defined with a black square; (c) red square indicates the location of the study area on aerial imagery (image source: © Google Earth; Digital Globe 2020); (d) perspective view of the slope, (e) generated 3D point cloud of site 2, (f) and generated 3D point view of site 1.

Figure 1 .
Figure 1.(a) Black circle indicates the location of the study area at Northeastern Greece; (b) simplified geological map (1:50,000 scale) of the wider area modified by [39], where the study area is defined with a black square; (c) red square indicates the location of the study area on aerial imagery (image source: © Google Earth; Digital Globe 2020); (d) perspective view of the slope, (e) generated 3D point cloud of site 2, (f) and generated 3D point view of site 1.

Figure 2 .
Figure 2. Flowchart of the methodology followed in order to evaluate the rock mass rating (RMR) and assess the slope mass rating (SMR) of the rock mass in the area of Nestos, Greece (DS: Discontinuity set).

Figure 2 .
Figure 2. Flowchart of the methodology followed in order to evaluate the rock mass rating (RMR) and assess the slope mass rating (SMR) of the rock mass in the area of Nestos, Greece (DS: Discontinuity set).

Figure 3 .
Figure 3.The investigated area during the field survey.The blue colored Schmidt hammer (length 28 cm) is indicated as a scale.

Figure 3 .
Figure 3.The investigated area during the field survey.The blue colored Schmidt hammer (length 28 cm) is indicated as a scale.

Figure 4 .
Figure 4. Validation of the structure from motion (SfM)-based developed point cloud through the comparison of the orientation of discontinuity sets defined based on geological compass and the generated 3D model.

Figure 5 .
Figure 5. (a) Stereonet of the density of the normal vector's poles and its corresponding planes based on Discontinuity Set Extractor (DSE); (b) discontinuity sets at site 1, colored as blue ( 1 ), green ( 2 ) and red (  3 ); (c) Stereographic projection of concentration lines of discontinuity poles from CloudCompare™ approach; (d) Stereographic projection of concentration lines of discontinuity poles from field survey based on geological compass

Figure 5 .
Figure 5. (a) Stereonet of the density of the normal vector's poles and its corresponding planes based on Discontinuity Set Extractor (DSE); (b) discontinuity sets at site 1, colored as blue (J 1 ), green (J 2 ) and red (J 3 ); (c) Stereographic projection of concentration lines of discontinuity poles from CloudCompare™ approach; (d) Stereographic projection of concentration lines of discontinuity poles from field survey based on geological compass.
Geosciences 2020, 10, x FOR PEER REVIEW 9 of 19 evaluation of RMRbasic remain the same, the average value of RMRbasic is 78 and the rock mass is classified as good.

Figure 6 .
Figure 6.Statistical analysis of normal spacing of discontinuity sets extracted from the DSE.

Figure 6 .
Figure 6.Statistical analysis of normal spacing of discontinuity sets extracted from the DSE.

Figure 7 .
Figure 7. (a) Stereonet of the density of the normal vector's poles and its corresponding planes based on DSE, (b) Stereographic projection of concentration lines of discontinuity poles from CloudCompare™ approach, (c) discontinuity sets at site 2, colored as blue (J 1 ), green (J 2 ), and red (J 3 ), and (d) J 1 DS, (e) J 2 DS, (f) J 3 DS.

Figure 8 .
Figure 8. Statistical analysis of normal spacing of discontinuity sets extracted from the DSE.

Figure 8 .
Figure 8. Statistical analysis of normal spacing of discontinuity sets extracted from the DSE.

Geosciences 2020 , 19 Figure 9 .
Figure 9. Sectors of the rock face colored based on the slope direction.Site 1 and 2 are also shown.

Figure 9 .
Figure 9. Sectors of the rock face colored based on the slope direction.Site 1 and 2 are also shown.

Figure 10 .
Figure 10.(a) Delineation of likely to failure blocks, (b) size characteristics of the block b, (c) size characteristics of the source area c.

Figure 10 .
Figure 10.(a) Delineation of likely to failure blocks, (b) size characteristics of the block b, (c) size characteristics of the source area c.

Figure 11 .
Figure 11.2D rockfall trajectories.(a) profile_1, (b) profile_2, (c) focus on the lower part of rock face showing the fall track reaching the railway (profile_2).Light green color refers to the marble bedrock and green color to the debris layer.Location of source areas (seeders) are indicated by blue crosses.

Figure 11 .
Figure 11.2D rockfall trajectories.(a) profile_1, (b) profile_2, (c) focus on the lower part of rock face showing the fall track reaching the railway (profile_2).Light green color refers to the marble bedrock and green color to the debris layer.Location of source areas (seeders) are indicated by blue crosses.

Table 1 .
Assigned parameters of the discontinuity sets used for the estimation of the RMR basic.

Table 2 .
Discontinuity sets and their characteristics in terms of dip and dip direction resulted from field-based and point cloud-based structural analysis.

Table 2 .
Discontinuity sets and their characteristics in terms of dip and dip direction resulted from field-based and point cloud-based structural analysis.

Table 3 .
Discontinuity sets and their characteristics in terms of dip and dip direction resulted from point cloud-based structural analysis.

Table 3 .
Discontinuity sets and their characteristics in terms of dip and dip direction resulted from point cloud-based structural analysis.

Table 4 .
Evaluation of SMR parameters of the discontinuity sets and formed wedges at sector 1.

Table 5 .
Evaluation of SMR parameters of the discontinuity sets and formed wedges at sector 2.

Table 4 .
Evaluation of SMR parameters of the discontinuity sets and formed wedges at sector 1.

Table 5 .
Evaluation of SMR parameters of the discontinuity sets and formed wedges at sector 2.