Slope stability and rockfall assessment of volcanic tuffs using RPAS with 2-D FEM slope modelling

Steep, hardly accessible cliffs of rhyolite tuff in NE Hungary are prone to rockfalls, endangering visitors of a castle. Remote sensing techniques were employed to obtain data on terrain morphology and to provide slope geometry for assessing the stability of these rock walls. A RPAS (Remotely Piloted Aircraft System) was used to collect images which were processed by Pix4D mapper (structure from motion technology) to generate a point cloud and mesh. The georeferencing was made by Global Navigation Satellite System (GNSS) with the use of seven ground control points. The obtained digital surface model (DSM) was processed (vegetation removal) and the derived digital terrain model (DTM) allowed cross sections to be drawn and a joint system to be detected. Joint and discontinuity system was also verified by field measurements. On-site tests as well as laboratory tests provided additional engineering geological data for slope modelling. Stability of cliffs was assessed by 2-D FEM (finite element method). Global analyses of cross sections show that weak intercalating tuff layers may serve as potential slip surfaces. However, at present the greatest hazard is related to planar failure along ENE–WSW joints and to wedge failure. The paper demonstrates that RPAS is a rapid and useful tool for generating a reliable terrain model of hardly accessible cliff faces. It also emphasizes the efficiency of RPAS in rockfall hazard assessment in comparison with other remote sensing techniques such as terrestrial laser scanning (TLS).

Most rockfall hazard publications deal with hard, wellcemented rocks, such as limestone (Samodra et al., 2016) or various other types of sedimentary rocks (Michoud et al., 2012), such as igneous or metamorphic rocks.In contrast, very few previous studies have dealt with cliff face stability and rockfall hazard of low-strength rock such as volcanic tuffs (Fanti et al., 2013;Margottini et al., 2015).Volcanic tuffs are very porous rocks and are prone to weathering (Arikan et al., 2007).While the current paper deals with a low-strength pyroclastic rock, it has a slightly different approach to cliff stability analysis.In this study, slope stability is assessed by using a combination of remote sensing techniques, field measurements and laboratory testing of tuffs with 2-D FEM (finite element method) analyses of slopes.In contrast to other case studies, this study operates on a smaller scale and studies the possibilities of wedge and planar fail-ures.More specifically, in our context, the cliff face is unstable as is evidenced by falling blocks.Due to rockfall hazard, the small touristic pathway was closed to avoid causalities.The current paper analyses the cliff faces by condition assessment and stability calculations.Thus, this research provides an assessment of how RPAS-based images and photogrammetric processing can be used to derive a surface model at sites that are difficult to access.The paper also demonstrates the combined use of photogrammetric, surveying and engineering geological methods at difficult ground conditions when assessing rock slope stability.

Study area
The study area is located on a medium-height mountain range in NE-Hungary, where a hardly accessible jointed rhyolite tuff cliff face was studied.On the top of the cliff, a touristic point, the Sirok Castle is located (Fig. 1).The steep rhyolite tuff hill has an elevation of 298 m at the transition area of two mountain ranges, Mátra and Bükk mountains.The tuff is very porous and prone to weathering (Vásárhelyi, 2002;Kleb and Vásárhelyi, 2003;Török et al., 2007).
Although the first castle was first constructed in the 13th century AD, due to war damages and reconstructions, the current structure encompasses wall sections representing different construction periods.In these days, the partially ruined walls have been restored, and the castle is open to tourists, but the southern slopes are closed due to rockfall hazard.The hill represents a rhyolite tuff that was formed during the Miocene volcanism (Badenian-Early Pannonian period).The cliff face was formed due to late Miocene volcanic activity and is a part of the Inner Carpathian volcanic chain.The geological map of the closer area clearly reflects the dominance of pyroclastic rocks, with isolated occurrences of Oligocene and Triassic rocks (Fig. 1).The cliffs are steep and display several joints and discontinuity surfaces.The present study focuses on the southern hillslope of the castle hill, where major rockfalls occurred in the recent past (Fig. 2).The study area is divided into smaller units, where RPAS and rockfall hazard assessment analyses were carried out (Fig. 3).

Materials and methods
The research uses two major methods: (i) RPAS and (ii) engineering geology.The applied methods are summarized in a flow diagram displaying the combination and links between the two methods (Fig. 4).The flow chart has four major realms that have both RPAS and engineering geological units: (i) preparation, (ii) field survey, (iii) data processing and calculation and (iv) interpretation.The RPAS line is described in detail in the next part of the paper, but is also linked to previous publications providing overview of image acquisition, image processing and interpretation (Civera et al., 2012;Westoby et al., 2012;Remondino et al., 2014).The engineering geological part of the flow chart is also explained below and has strong links to publications describing the application of RPAS to landslide characterization and rock slope stability assessment (Niethammer et al., 2012;Tannant, 2015).the exposure time was set to 1/1400 s and the images were compressed at a rate of 4.5 bits pixel −1 .There were three imaging flights: two around noon and one at about 17:00.The flying times were 13, 12 and 13 min at which 390, 365 and 419 images were captured.All 1174 images were involved in the photogrammetric object reconstruction (Fig. 5).

RPAS data acquisition and terrain modelling
The photogrammetric reconstruction has been done using Pix4Dmapper (Pix4D, 2017), which is based on structure from motion (SfM) technology (Lowe, 2004;Westoby et al., 2012;Danzi et al., 2013).SfM automatically identifies tie points considering initial requirements (e.g.preliminary image centre positions, time stamps) (Table 1).Camera calibra-tion was executed during post-processing, and no prior calibration was needed (Pix4D, 2017).After the image alignment, the image projection centres and attitudes can be observed (Fig. 5).Then 19.3 million points were obtained by the photogrammetric reconstruction, which was appropriate for the engineering geological application.The technology allows a higher resolution to be obtained, but it was not necessary.The average point density is about 670 points m −3 , but there are areas with double point density.
For georeferencing, particular tie objects were measured by the Global Navigation Satellite System (GNSS).The used GNSS receiver was a Leica CS10 with a Gs08plus antenna (GS08, 2014;CS10, 2014).The measurement was done in RTK mode supported by the Hungarian RTK network (RTKnet, 2013).There were seven measured ground control points (Fig. 6) (GCPs); the mean 3-D measurement accuracy was 4.9 cm (minimal value was 2 cm, maximal value 9 cm).The RPAS technology has produced a considerable number of data points (observations).Since this point cloud is difficult to manage due to its size and heterogeneous point spacing, it requires a sophisticated resampling step, which was done using CloudCompare, and the spatial resolution of the point cloud was set to 1 cm.
The RPAS data collection was validated by the use of terrestrial laser scanning.The necessary data were captured by two scanners: a Faro Focus S 120 3D (Faro, 2016) and a Z + F Imager 5010C (Z + F, 2014).The terrestrial laser scanning was executed on the same day as the RPAS mission.The raw point cloud measured by the Faro scanner contained 1.9 billion points, while the Z + F point cloud had 0.8 billion points.Both point clouds included X, Y and Z coordinates, intensity and RGB colour values.RPAS-and TLSbased point clouds were compared with CloudCompare software (CloudCompare, 2014) (Fig. 7).
As can be noted in Fig. 7, the largest difference between the two sources is almost less than 10 cm and the majority of the differences are about 1 cm.The point cloud was then imported into Geomagic Studio 2013 (GeomagicStudio, 2013) and meshed.The triangle side length was 5-7 cm.To support the engineering geological survey, several horizontal and vertical sections were derived in Geomagic DesignX 2016 (Ge-omagicDesignX, 2016); these profiles were exported in CAD format.
The next step was to make cut-offs focusing only on the cliffs; this was done by CloudCompare, followed by the points being exported in LAS format (LAS, 2012).The exported points could then be imported into SAGA GIS 2.1.2(Conrad et al., 2015), where the necessary DTMs were created using an inverse distance weighting (IDW) algorithm (IDW, 2013).The derived DSM-grids have 5 cm spatial resolution, which is adequate for morphologic analyses (Fig. 6) and suits to slope stability analysis.The morphology analysis has concentrated on the catchment area (CA) (Costa-Cabral and Burges, 1994;Haas et al., 2016) (Fig. 4), although several other morphological indices (e.g.topographic wetness index, stream power index) were derived.These indices express the potential relationship between surface geometry and geological parameters.

Engineering geology and slope stability analysis
Geological data and written resources (Balogh, 1964;Haas, 2013, Lukács et al., 2015) provided background information for the planning of engineering geological field surveys (Fig. 4).Major lithotypes were identified and described and geological profiles were recorded during the engineering geological field surveys (Fig. 4).Rock joints, discontinuity surfaces and fault systems were measured by using compass and structural geological software applied in mobile phone.The structural geological data were analysed by Dips software.Strength parameters were assessed on-site by using  a Schmidt hammer.Ten rebound values were measured on each surface and mean values and SDs were also calculated.This method has been used previously to gather rapid data on the rock strength of cliff faces (Margottini et al., 2015).The data set was compared to rock mechanical laboratory tests.Samples for laboratory analyses were collected on-site (Fig. 4).Major rock mechanical parameters were measured under laboratory conditions on cylindrical specimens.These were drilled from blocks and cut into appropriate sizes using cutting disc.The sizes of the tested specimen were made according to EN (European Norm) on air-dried and watersaturated samples.The specimens were grouped according to their bulk density and the propagation speed of the ultrasonic pulse wave.Strength parameters such as uniaxial compressive strength and indirect tensile strength (Brasilian) were measured according to relevant EN standards and ISRMsuggested methods (International Association of Rock Mechanics).Modulus of elasticity was also calculated (Table 2).
The generalized Hoek-Brown failure criterion (Hoek et al., 2002) was used to determine strength parameters of the rock mass.Altogether, 53 cylindrical test specimens were used for the tests.
The falling blocks can endanger tourists on the footpath below the castle on the southern slopes; therefore the stability analysis of the rocky slopes was focused on this part of the cliff (Fig. 3).First, the rock mass failure was analysed with by the RocFall FEM software of the Rocscience (RS2).The steepest sections were determined from the terrain model obtained from RPAS data.The GSI (Geological Strength Index) values of the rock masses were determined according to Marinos et al. (2005).The global stability of the hillslope of selected sections was calculated with RS2 software.Since the rhyolite tuff is a weak rock with few joints, the rock mass failure and the failure along discontinuities were also analysed.This kinematic analysis had been done with a stereographic tool.The orientations of main joint sets were obtained from DSM model-based morphological analyses with the use of catchment area tool.It assumes that major flow paths are related to joints; i.e. the fracture system controls the drainage pattern (Costa-Cabaral et al., 1994).At accessible areas joints and fractures were also measured onsite on the southern and south-eastern parts of the hillslope.Additional control field measurements were also made in the underground cellar system of the castle, where the tuff is also exposed.The Dips software was used for the kinematic analysis.The direction of the hillslopes and the direction of the discontinuities were compared to determine the location of the potential hazardous failure zones on the hillside.Stereographic plots were generated showing the possible failure planes for all slope directions, and the safety factor of the possible planar failure was calculated by RockPlane software.Wedge failure was modelled by Swedge software.Toppling failure due to geological and geomorphological conditions cannot occur.Slope stability calculations and stability assessment formed the last part of the engineering geological analyses (Fig. 4).

Results
The rhyolite tuff faces consist of moderately bedded ignimbritic horizons and also brecciated lapilli tuffs according to our field observations (Fig. 6).The topmost 10 m of the cliff face, which was modelled from slope stability, com-prises three main horizons and can be modelled as a "sandwich structure".The lower and the upper parts are formed by thick pumice containing lapilli tuffs.These beds enclose nearly 2 m of well-bedded less welded fine tuff and brecciated horizons (Fig. 8).
By combining and comparing all measured data of discontinuities and joints using DSM and its derivative (Fig. 6.) and morphological index (Fig. 9), the joint orientation was outlined.The filed survey validated the obtained structural geological conditions, and six main joint sets (with dip angle/dip directions of 85/156, 88/312, 79/110, 81/089, 82/064, 61/299) were identified with prevailing a NE-SW direction.
The laboratory tests of tuffs provided the input data for stability analysis for the two main lithologies: upper and lower unit of lapilli tuff and middle unit of less welded tuff (Table 3).In the model calculations GSI = 50 was used.
The results of RS2 FEM analyses suggest that the global factor of safety is SRF = 1.27-1.71 in the studied sections (some of the selected sections are shown in Fig. 3).The aim of the analysis is to determine the critical strength reduction factor (SRF) that can be considered the safety factor of the slope (Fig. 10).The SRF factor is influenced by the weak tuff layer (marked by B-D in Fig. 8), which has very low shear strength compared to the lapilli tuff.Colours in Fig. 10.represent the total displacements as a result of the shear strength reduction analysis.Thus these are not real displacements of the hillslope.The figure demonstrates only how the failure of the slope would occur with reduced shear strength parameters.Our failure analyses have proved that the bottom of the slip surface would be in this weaker layer and could lead to a larger mass movement.
Other failure modes that were studied are planar failure and wedge failure, which are often controlled by joints and discontinuities.According to data obtained from remote sensing and according to field measurements there was no uniform spacing of the discontinuities.Stereographic plots with possible failure planes for all slope directions (Fig. 11) indicate that the most hazardous part of the slope is the one where the plane orientation is 75/75.The calculated factor of safety (FS = 1.15) implies a high probability of planar failure.
Three possible wedge failure modes were identified as being the most hazardous according to our calculations by Swedge software (Fig. 12).In these three cases the factor of safety was in the range of 1.38-1.94,representing the hazards of rockfalls along wedges delineated by different joint systems.

Discussion
There are three critical sets of input data in the modelling of rocky slopes: (i) terrain model and slope geometry, (ii) joint systems and (iii) strengths of rock mass.
To obtain the first, the slope geometry, RPAS-based surveying technique was used similarly to previous studies (Giordan et al., 2015).In many previous studies RPAS mission was performed following a flight plan (Eisenbeiss, 2008;Lindner et al., 2016).In our case no flight plan was made prior to the mission, and the windy conditions were also not favourable for a pre-programmed flight route.The flight was  manually controlled and the skilled personnel with a first person view tool controlled the image acquisition.The crucial points were the necessary overlaps of images and the matching.The overlaps were ensured by three consequent flights around the study area that provided a dense network of image acquisition locations (Fig. 5).The obtained 1174 images covered the entire study area with appropriate overlap.The number of images is reasonable, since in previous studies 400 images were taken for a smaller translational rockslide by a GoPro Hero 3 Black camera (Tannant et al., 2017) or approx.400-900 images with a higher resolution camera (18 MP) for a landslide area that is approximately five times larger than this one (Lindner et al., 2016).The Pix4Dmapper software (SfM) was used to identify key points.Nearly 10 000 key points were found on each image (Table 1), which is considered a sufficient number for appropriate matching (Remondino et al., 2014).A camera selfcalibration tool and rolling shutter effect correction were the other key features of this software that allowed the images of GoPro Hero 3+ camera to be processed.The obliquity of images (Rupnik et al., 2014) and the density of obtained data (Remondino et al., 2011(Remondino et al., , 2014;;Rupnik et al., 2015) are crucial in the applicability and accuracy of these images.These were also managed by Pix4Dmapper.The GNNS system and the ground control points (Fig. 6) allowed the rocky slope to be georeferenced.Our results show that the mean 3-D accuracy was 4.9 cm.It is comparable with the ground resolution of 1-3.5 cm pixel −1 of an Italian rockslide survey (Tannant et al., 2017) or 3.3-4.1 cm (Neugirg et al., 2016).This resolution was appropriate for creating a reliable digital surface model.
RPAS-based data were also validated with TLS measurements.The co-use of these remote sensing tools has been previously well documented for other applications such as soil roughness (Milenkovich et al., 2016) in erosion detection (Neugirg et al., 2016) or in cultural heritage (Eisenbeiss and Zhang, 2006).The RPAS-obtained data validation was performed by comparing the two point clouds obtained by RPAS and TLS.The surfaces were resampled in order to homogenize the spatial resolution.The point densities have been tested in CloudCompare, as a unit sphere of volume of 1 m 3 was defined where the points can be counted and then the sphere could be moved along the whole surface.The differences in point clouds are less than 10 cm (Fig. 7), which is considered a reasonably good match in terrain modelling (Neugirg et al., 2016).This computation proved that the average point density in both point clouds are practically the same, although the RPAS densities are more homogeneous, while the TLS has denser point clouds close to the scanning stations, as it was expected on the basis of previous works (Naumann et al., 2013).
Another aspect causing some differences between the two data sets is that the image-based reconstruction is performed by interest operators, usually SIFT (scale-invariant feature transform) or similar computer vision operators (Lowe, 2004).These operators are generally sensitive to intensity jumps, points or corners, and textural changes in the input images.If the image resolution is not adequate or the object is locally "smooth", these operators do not return with surface points and the output of the reconstruction has some filtered effect.Fortunately, the surface reconstruction quality in RPAS processing resulted in minor, ignorable smoothing effects.Comparing the two data sets, it is clearly proven that the geometric resolution of the RPAS-based digital surface model corresponds to the TLS one, offering very similar quality data (Fig. 7).It is necessary to note that vegetation can hamper TLS measurements (Prokop and Panholzer, 2009;Scaioni et al., 2014;Tannant, 2015) and thus limit the comparison of RPAS-and TLS-obtained data (Milenkovic et al., 2016).In our case most of the study area was bare and if vegetation occurred it was manually removed.
The documentation of joint system and discontinuities are crucial for rock wall stability assessment (Tannant et al., 2017).A field survey can only provide reliable data on joint orientation of accessible areas (Margottini et al., 2015); however, the joint system that was found on inaccessible cliffs was not detectable.To overcome this problem, RPASgenerated images were used; the frequency of joints was observed based on these images like in previous studies by Assali et al. (2014), Martino and Mazzanti (2014) and Margottini et al. (2015).The required resolution for joint frequency is of the order of 10 cm, rarely 1 cm (Tannant, 2015;Tannant et al., 2017).The RPAS technique allows for plane surface geometries; however, many joints are not plane surfaces and there are sets in shadows that are difficult to visualize.Thus, RPAS can be used to outline the strike of major joints, but it might cause problems when it comes to the determination of dip and the displacement along the fault planes (e.g.slickensides).
Our field tests indicate that the application of a Schmidt hammer in rock strength analysis is limited when it comes to the analysis of low-strength rocks, such as volcanic tuff (Aydan and Ulusay, 2003).As a consequence, laboratory analy-ses of samples were also required to obtain reliable strength parameters.To measure the strength and to understand the weathering characteristics, samples were taken representing different stratigraphic positions.Our lab test data (Table 3) clearly indicate that a low-strength unit is found in the studied sections (unit marked by B-D in Fig. 8).Whether the low strength of this zone is related to differential weathering (Török et al., 2007) or if it is associated with inherited weakness (micro-fabric) is not clear.This layer is a potential failure zone as was shown by slope stability calculations.A similar intercalation of pumice-rich layered deposits was modelled by Damiano et al. (2017).They found that a pumice-reach weak zone is prone to rainfall-induced landslides.Our results are in good correlation with these findings since the studied rhyolitic volcanic tuff was also proved to be very prone to weathering.A loss in tensile strength of 60 % was measured under simulated laboratory conditions (Stück et al., 2008).Weathering processes have long been known to induce landslides and cause slope stability problems in various lithologies and especially in pyroclastic rocks (Chigra, 2002;Fanti et al., 2013).At the studied rhyolite tuff cliff face, it was shown that a joint system is responsible for slope instabilities: planar and wedge failures were found (Figs. 11 and 12).These failure modes are common in hard jointed cliff faces (porphyry in Agliardi et al., 2013;mica-schist in Tannant et al., 2017;limestone in Feng et al., 2014).Our study demonstrates that joint systems have significant control on slope stability, not only in hard rock lithologies but also in weak tuffs.It is in line with Fanti et al. (2013) and Margottini (2015), since in Italy and in Giorgia rock walls of volcanic tuffs suffered landslides.The kinematic analysis of tuff rock walls of Tuscany (Fanti et al., 2013) also demonstrated that wedge failure and planar failure are the most common failure mechanisms of tuff cliff faces.

Conclusions
The manually controlled flights of RPAS provided excellent information on slope geometry of highly dissected and inaccessible slopes.The necessary overlap between images was ensured by three flights over the small area and by skilled personnel using a first person view system with a synchronous image transfer.The obtained data were managed by Pix4Dmapper (SfM) software allowing the identification of nearly 10 000 key points per image.The TLS-based point clouds proved to be good tools for validating the accuracy of images and data sets of manually controlled RPAS.In our study the maximum difference between the two point clouds was less than 10 cm but was mostly around 1 cm.RPAS collected images and the point cloud-based digital surface model and the catchment area method especially allows the detection of joint system (mainly strikes and partly dips but not slickensides) but field validation and field measurements of accessible joints and faults are recommended to justify joint orientation.The obtained digital surface model was accurate enough to allow cross sections for rock wall stability calculations.The lithology and physical parameters of the studied steep cliffs are not uniform and intercalations of weak layers of vitric tuff and volcanoclastic breccia were found.According to 2-D FEM modelling the intercalating low-strength layer is the one at which potential slip surface can develop, causing larger-scale mass movements, but at present it has low probability.Joint systems have a crucial role in the stability of the studied rhyolite tuff cliff faces.The greatest hazard is related to planar failure along ENE-WSW joints and to wedge failure.
Data availability.The data can be requested from the corresponding author.It will not be publicly available until the project terminates.
Competing interests.The authors declare that they have no conflict of interest.Special issue statement.This article is part of the special issue "The use of remotely piloted aircraft systems (RPAS) in monitoring applications and management of natural hazards".It is a result of the EGU General Assembly 2016, Vienna, Austria, 17-22 April 2016.

Figure 2 .
Figure 2. Studied southern cliff faces: (a) image of the castle obtained by RPAS with marked details, (b) distant view of the eastern part of the cliff section, (c) steep cliffs dissected by joints, (d) vertical to subvertical cliff face with steep joints and traces of rockfall and (e) weathered rounded cliff with larger tafoni.

Figure 3 .
Figure 3. Location of the illustrations in the paper and the rockfall-affected area.Red dotted and dashed line represents zones affected by rockfall.Yellow dashed lines 1 to 5 mark the sections where slope stability was calculated by using 2-D FEM model (Fig. 8).Dotted lines indicate the areas shown in Figs. 5 and 9.

Figure 4 .Figure 5 .
Figure 4. Flow chart showing the methods and obtained data set of this paper indicating the interrelationship between RPAS and slope stability assessment (see details in the text).

Figure 6 .
Figure 6.Colour digital surface model with 1 m contour line interval of the study area.The solid black dots show the ground control points, while the red dotted and dashed line represents zones affected by rockfall.

Figure 7 .
Figure 7. Differences between RPAS and TLS point clouds by CloudCompare shown in metres (modus of differences is at about 0.01 m).

Figure 8 .
Figure 8. Lithologic column of Sirok Várhegy showing the modelled topmost 10 m section of the hill (letters refer to lithologic units).

Figure 9 .
Figure 9. Top view of the cliff (see the location in Fig.3) obtained with RPAS and the catchment area diagram obtained from DSM analysis by using a catchment area module(Costa-Cabral et al., 1994).The latter was used for joint pattern recognition.Numbers refer to major joint systems marked on the catchment area map and on rose diagram of the field measurements and DSM data set.

Figure 10 .
Figure10.The results of the global stability analysis of the slopes (sections 3 and 4 in Fig.3).Total displacements are marked from blue to red (lithology is indicated by letters A-D; note the weak zone marked by B-D; description of lithologies is given in Fig.8.)

Figure 12 .
Figure 12.Examples for the kinematic analysis of wedge failures (main joint sets are marked by red circles 1-6 m).

Table 1 .
Image processing data.
Mean number of key points per image 22 676 Mean number of matched key points per image 9546 Mean reprojection error (pixel) 0.176 Time for SfM processing 40 m:20 s Time for densification (point cloud) 05 h:30 m:24 s

Table 2 .
Rock mechanical tests and relevant standards.