UAV-enabled reconnaissance and trajectory modeling 1 of a co-seismic rockfall in Lefkada 2

UAV-enabled reconnaissance and trajectory modeling 1 of a co-seismic rockfall in Lefkada 2 Charalampos Saroglou, 3 Pavlos Asteriou 4 Dimitris Zekkos 5 George Tsiambaos 6 Marin Clark 7 John Manousakis 8 Department of Geotechnical Engineering, School of Civil Engineering, National Technical 9 University of Athens 10 Department of Civil and Environmental Engineering, University of Michigan, USA 11 Department of Earth and Environmental Science, University of Michigan, USA 12 Elxis Group, S.A, Athens, Greece 13


Introduction
Active faulting, rock fracturing and high rates of seismicity contribute to common rockfall hazards in Greece.Rockfalls primarily damage roadways and houses (Saroglou, 2013) and are most often triggered by rainfall and secondly seismic loading.Additionally in recent years, some rockfalls have impacted archaeological sites (Marinos & Tsiambaos, 2002, Saroglou et al., 2012).The Ionian Islands, which includes Lefkada Island, experience frequent M w 5-6.5 earthquakes, as well as less frequent larger (up to 7.5) earthquakes.The historical seismological record is particularly well constrained with reliable detailed information for at least 23 events since 1612, which induced ground failures at the island of Lefkada (Papathanasiou et al. 2005) and an average of a damaging earthquake every 18 years.In the recent past, a M w 6.2 earthquake occurred on August 14 2003 and was located offshore the NW coast of Lefkada.Significant damage was reported, particularly in the town of Lefkada, where a PGA of 0.42g was recorded.Landslides, rockslides and rockfalls occurred along the western coast of the island (Karakostas et al. 2004, Papathanasiou et al., 2012).
On November 17 th 2015, an M w 6.5 earthquake struck the island of Lefkada and triggered a number of landslides, rockfalls and some structural damage.The most affected area by large rockslides was the west coast of the island, especially along its central and south portion which are popular tourist destinations (Zekkos et al., 2017).
The landslides completely covered the majority of the west coast beaches and damaged access roads.A rockfall in Ponti village, which lies in the southeast side of Lefkada near the Gulf of Vassiliki, was triggered during this earthquake and was responsible for one of two deaths caused by the earthquake.Of particular interest is the very long travel path of the rock block, which was about 800 m in plan view from the point of detachment to the end of its path.Near the end of the rock fall path, the block impacted a family residence, penetrated two brick walls and killed a person in the house.The block exited through the back of the house and came to rest in the property's backyard.
The Ponti village rockfall site is characteristic of earthquake induced rockfall and an example of how seismically-induced rockfall impacts human activities.It also exemplifies the limitations of common 2D rockfall analysis to predict specific aspects of the rockfall trajectory as measured by field evidence.In order to create a highly accurate model of the rockfall propagation in 2D and 3D space, the rock path and the impact point on the slope was identified by a field survey.The study was performed using an Unmanned Aerial Vehicle (UAV) with an ultra-high definition (UHD) camera, which produced a high-resolution orthophoto and a Digital Surface Model (DSM) of the terrain.The orthophoto was used to identify the rolling section and the bouncing points of the rock along its trajectory.The high-resolution DSM made it possible to perform kinematical rebound analysis and a 3D rockfall analysis.

Ponti rockfall -site conditions
The slope overhanging Ponti village is formed in limestone and has a maximum height of 600 m and an average slope angle of 35 0 to 40 0 (Figure 1).The geological formations at the Ponti rockfall site are limestones covered by moderately cemented talus materials.The thickness of the talus materials ranges between 0. and bottom of the surveyed area, the 3D model is georeferenced to a specific coordinate system.

High-resolution Orthophoto
A 5cm pixel size orthophoto was generaated based on thhe methodology outlined earlier.As showin Fig. 3, the rolling section and the bouncing locations of the rock block throughou its course were identified.The rolling section is discerned as a continuous and largely linear mark left in the densely vegetated terrain that is indicative of the damage caused.Impact points that are part of the bouncing section of the rock, are identified as circular to ellipsoidal bare earth craters with no disturbance in between.The last bouncing point before impacting the house is clearly identified on the paved road.The plan view ortho-imagery, along with the original footage of the video collected was crucial to the qualitative identification of these features.The alternative, i.e., land-based, conventional field reconnaissance was practically impossible tto perform in the the densely vegetated and steep terrain.

Digital Surface Model
A profile section and a 10 cm Digital Surface Model (DSM) paired with the plan view orthophoto were first developed (Manousakis et al., 2016) and made possible the identificaiton of terrain features such as structures, slope benches or high trees, that could affect the rock's path downhill.However, this resolution of the DSM proved to be not only unnecessarily high and thus difficult to manipulate in subsequent rockfall analyses, but also resulted in numerical instabilities during the rockfall analyses.
Therefore, a downscaled 2 m DSM was produced for the rockfall analysis.This was implemented through an aggregate generalization scheme where each output cell is assigned the minimum of the input cells that are encompassed by that cell.In

Seismic acceleration
The epicenter of the earthquake according to the National Observatory of Athens, Institute of Geodynamics (NOA) is located onhore near the west coast of Lefkada.
The causative fault is estimated to be a near-vertical strike-slip fault with dextral sense of motion (Ganas et al., 2015(Ganas et al., , 2016)).Based on the focal mechanism study of the earthquake it was determined that the earthquake was related to the right lateral Kefalonia-Lefkada Transform Fault (KLTF), which runs nearly parallel to the west coasts of both Lefkada and Kefalonia island, in two segments (Papazachos et al. 1998, Rondoyanni et al. 2012).The previous earthquake in this zone occurred in August 2003 with a magnitude of 6.2.
A strong motion station recorded the ground motions in the village of Vasiliki located at a distance of 2.5 km from the site.The ground motion characteristics of the recording are summarized in Table 1 and are presented in Figure 4, according to an ITSAK preliminary report (ITSAK, 2016).In comparison with the recordings at other locations in Central Ionian, it was evident that the strongest acceleration was encountered in Vasiliki area.

Topography effect
Peak ground acceleration along the rock slope is the intensity of base shaking modified by site and topographic effects (Mavrouli et al., 2009).In the present case,

Assessment of rock block's initial velocity
The initial horizontal velocity of the block, at the time of detachment, was calculated considering equilibrium of the produced work and the kinetic energy according to equation 1.
where PGA sf is the acceleration on the slope at the location of detachment and s the initial displacement of the block in order to initiate its downslope movement.
The initial horizontal velocity was calculated equal to 0.67 m/sec, considering a displacement in the order of s = 0.05 m.The vertical component of the initial velocity is assumed to be zero.

Trajectory analysis
In order to estimate the possible rock paths and design remedial measures, simulation programs are used in design practice, which are mostly based on the lumped-mass analysis model.The trajectory of a block is modelled as a combination of four motion types; free falling, bouncing, rolling and sliding (Descoeudres and Zimmermann, 1987).

Modelling the response to an impact
The most critical input parameters are the coefficients of restitution (COR), which control the bouncing of the block.In general, the coefficient of restitution (COR) is defined as the decimal fractional value representing the ratio of velocities (or impulses or energies; depending on the definition used) before and after an impact of two colliding entities (or a body and a rigid surface).When in contact with the slope, the block's magnitude of velocity changes according to the COR value.Hence, COR is assumed to be an overall value that takes into account all the characteristics of the impact; including deformation, sliding upon contact point, transformation of rotational moments into translational and vice versa (Giani, 1992).
The most widely used definitions originate from the theory of inelastic collision as described by Newtonian mechanics.For an object impacting a rocky slope (Figure 5), which is considered as a steadfast object, the kinematic COR (v COR ) is defined according to Eq. 2.
where v is the velocity magnitude and the subscripts i and r denote the trajectory stage; incident (before impact) and rebound (after impact) respectively.
Two different mechanisms participate in the energy dissipation process; energy loss normal to the slope is attributed to the deformation of the colliding entities, and in the tangential direction is due to friction between them.Therefore kinematic COR has been analyzed to the normal and tangential component with respect to the slope surface, defining the normal (n COR ) and the tangential (t COR ) coefficient of restitution (Eq. 3 and 4 respectively).
where the first subscript, n or t denotes the normal or the tangential components of the velocity respectively.Usage of the lump-mass model has some key limitations; the block is described as rigid and dimensionless with an idealized shape (sphere); therefore the model neglects the block's actual shape and configuration at impact, even though it is evident that they both affect the resulting motion.

Rockfall path characteristics
23 impact points were identified on the slope surface (Figure 6).Their coordinates are presented in Table 2, along block's path starting from the detachment point (where x=0).
The apparent dip of the slope at impact positions was measured from the topographic map; on each impact point a line was set with a length twice the block's mean dimension, oriented according to preceding trajectory direction.Moreover, the impact point was expanded on the topographic map to a rectangular plane with a side twice as much the mean dimension of the block (Figure 7).This plane was then oriented so that one side coincides with the strike direction and its' vertical side toward to the dip direction.Thus, direction difference, ∆φ, was measured by the strike direction and the preceding path and deviation, e, was measured as the angle between pre and post impact planes (Asteriou & Tsiambaos, 2016).Having a detailed field survey of the trajectory path, a back analysis according to the fundamental kinematic principles was performed in order to back-calculate the actual COR values.

Kinematic analysis and assumptions
The 23 impact points identified on the slope comprise a rockfall path of 22 parabolic segments.The vertical and horizontal length of each segment is acquired by subtracting consecutive points.Since no external forces act while the block is in the air, each segment lays on a vertical plane and is described by the general equation of motion as: ( where: θ the launch angle from the horizon and v the launch (initial) velocity (Figure 8).
Since no evidence can be collected regarding launch angle and velocity, innumerable parabolas satisfy Eq. 5.However, θ is bound between -β and 90 o , so in order to acquire realistic values for the initial velocity, its sensitivity for that given range was addressed (Figure 9).
For the case presented in Fig. 9 (the first parabolic segment) it is seen that for the majority of the release angles, initial velocity variation is low and ranges between 7.2 and 12ms -1 .Additionally, the relationship between release angle and initial velocity is expressed by a curvilinear function, thus a minimum initial velocity value along with its release angle (denoted hereafter as θ cr ) can be easily acquired.
Given the minimum initial velocity and the critical release angle for each parabolic segment, the impact velocity and angle can be calculated.Afterwards, normal and tangential velocity components according to the apparent dip of the impact area, are calculated in order to evaluate COR values.Results are summarized in Table 3.

Coefficients of restitution
It is observed that v cor (Table 3) is slightly greater than one in 5 out of 22 impacts.
According to Eq. 3, this can only be achieved when impact velocity is less than rebound velocity.However, this indicates that energy was added to the block during contact, which is not possible according to the law of conservation of energy.Thus, impact velocity should be greater, which is possible if the launch velocity of the previous impact was more than the minimum, as assumed.
Omitting the impacts with V cor >1, it is observed that kinematic COR ranges between 0.55 and  Initially, a deterministic 2D rockfall analysis was performed using Rocfall software (RocScience, 2004).Considering an initial velocity of 0.67 m/sec, the falling rock primarily rolls on the slope and stops much earlier than its actual run out distance, approximately 400 m downslope from its starting point (Fig. 6; case 1).The restitution coefficients were n COR =0.35, t COR =0.85, which represent properties of bedrock outcrops according to the suggested values provided in the documentation of the software.The friction angle was set to zero.If the friction angle is set to φ=32 0 (as suggested by the software documentation), the rock travels downslope only 50 m.
A separate analysis was performed, with lower coefficients of restitution, resembling that of talus material on the slope (n COR =0.32, t COR =0.82, φ=30 0 ) as proposed by the suggested values provided in the documentation of the software.In this case, the rock block rolled only a few meters downslope.Therefore, it is evident that the actual rock trajectory cannot be simulated.
In order to simulate the actual trajectory as much as possible, various combinations of restitution coefficients and friction angle were considered.The closest match occurred for n COR =0.60 and t COR =0.85, while the friction angle was set to zero and no velocity scaling was applied.Only in such an analysis, the rock block reaches the house; with a velocity equal to v=18 m/s approximately (Fig. 6; case 2).According to the suggested values, these values for the coefficients correspond to a bedrock material (limestone).
In this case, the modelled trajectory is significantly different from the actual one.The main difference is that the block is rolling up to 200 m downslope while the actual rolling section is 400 m (as shown in Fig. 6).Furthermore the impacts on the ground in the bouncing section of the trajectory are considerably different in number ( 14versus 23) and in location from the actual ones.Finally, the bounce height of some impacts seems unrealistically high.(proposed in the manual), denoting talus with a larger percentage of fallen boulders.

3-D rockfall analysis
The block dimensions were considered equal to 2 m 3 and the shape of the boulder was rectangle.In order to simulate the initial velocity of the falling rock due to the earthquake an additional initial fall height is considered in the analysis, which for this case was equal to 0.5 m.
The energy line angles were recalculated from the simulated trajectories and it was determined that the energy line angle with highest frequency (39%) was 30-31 0 .
Based on the 3D analysis no rock blocks would impact the house, although the rock

Lateral dispersion & Deviation
Lateral dispersion is defined as the ratio between the distance separating the two extreme fall paths (as seen looking at the face of the slope) and the length of the slope (Azzoni and de Freitas 1995).According to Crosta and Agliardi (2004)  Lateral dispersion cannot be defined from the actual rockfall event in Ponti since only one path is available.Using the simulated trajectories from RockyFor3D, which are in the 3d space (Figure 13), a lateral dispersion of approximately 60% is shown in the middle of the distance between detachment point and the house.This is significantly higher compared to the findings of Azzoni and de Freitas (1995) and Agliardi and Crosta (2003).Moreover, based on the actual event and intuition, the lateral dispersion computed by RockyFor3D is extremely pronounced.Examining Figure 13, it is notable that the rock paths are severely affected by the topography factors.
Therefore, assessing lateral dispersion seems to be a case specific task.
Asteriou & Tsiambaos (2016) defined deviation (e) as the dihedral angle between the pre-and post-impact planes that contain the trajectory.They found that deviation is controlled by the direction difference ∆φ, the slope inclination and the shape of the block.For a parallel impact (i.e.∆φ=0 0 ) a spherical block presents significantly less deviation compared to a cubical.Additionally, deviation is equally distributed along the post-impact direction and reduces as the slope's inclination increases.On oblique impacts the block's direction after impact changes towards the aspect of slope and as ∆φ increases this trend becomes more pronounced.
Figure 14 presents deviation as a function of direction difference.It is noted that for parallel impacts deviation is also equally distributed along the post-impact direction.
As direction difference increases, deviation becomes positive, which means that the change of direction is following the direction of slope's aspect.These finding are in Based on the analytical analysis performed, it was found that the coefficients of restitution cannot be directly connected to the material type, nor can be considered as constants.The impact angle seems to pose a consistent effect on normal COR, which has been seen also in other recent relevant studies, but has not been incorporated yet on analyses models.
Performing a 2D rockfall analysis with the set of parameters recommended by the developers, was impossible to replicate the actual trajectory revealing some limitations in the present formulations.In an attempt to match the actual rock path to the analysis output, the friction angle of the limestone slope was considered equal to zero.However, the falling rock rolled on the slope and stopped much earlier than its actual runout distance while the impacts on the ground in the bouncing section of the trajectory were considerably different in number and in location compared to the actual ones.
Using the 3D analysis software, some rock trajectories better approximated the actual trajectory using the suggested values by the software developers, testifying that the 3D analysis can be more accurate than the 2D analysis.
Based on the aforementioned analyses it becomes evident that engineering judgement and experience must accompany the usage of commercial rockfall software in order to acquire realistic paths.One should never rest on the suggested set of parameters since the actual outcome can differ significantly, as demonstrated by this case study.
5 and 4.0 to 5.0 m.A few fallen limestone blocks were identified on the scree slope, with volumes between 0.5 and 2 m 3 .Based on the size distribution of these rocks on the slope, the average expected block volume would be in the order of 1 to 2 m 3 .Nat. Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2017-29,2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 February 2017 c Author(s) 2017.CC-BY 3.0 License.Introduction A quadrotor UAV (Phantom 3 professional) was deployed to reach the uphill terrrain that was practically inaccessible.The UAV was equiped with an Ultra-high definition (UHD) 12 MP camera and has the capacity to collect 4K video.The first objective of the UAV deployment was to find the inititiaion point of the rock and then identify the rockfall path (shown in Fig. 1).A particular focus on that part of the task was the identification of rolling and bouncing sections of the rockfall path.In addition, in order to generate a high-resolution orthophoto of the rockfall trajectory, aerial video imagery was collected, and the resulting digital surface model (DSM) was used to perform rockfall analysis.Structure-from-Motion (SfM) methodologies were implemented to create a 3D point cloud of the terrain and develop a 3D model by identifying matching features in multiple images.Compared to classic photogrametry methodologies, where the location of the observing point is well established, SfM tracks specific discernible features in multiple images, and through non-linear leastsquares minimisation (Westoby et al., 2012), iteratively estimates both camera positions, as well as object coordinates in an arbitrary 3D coordinate system.In this process, sparse bundle adjustment (Snavely et al., 2008) is implemented to transform measured image coordinates to three dimensionl points of the area of interest.The outcome of this process is a sparse 3D point cloud in the same local 3D coordinate system (Micheletti et al., 2015).Paired with GPS measurements of a number of control points (in this case 10 fast-static GPS points) at the top, middle Nat.Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2017-29,2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 February 2017 c Author(s) 2017.CC-BY 3.0 License.
local shaking intensity in terms of horizontal PGA was considered.The E-W component of acceleration was considered for the determination of the initial velocity.The peak ground acceleration (PGA) on the slope face (PGA sf ) was obtained by linear interpolation between the acceleration at the base (PGA b ) and at the slope Nat.Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2017-29,2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 February 2017 c Author(s) 2017.CC-BY 3.0 License.crest (PGA cr ).The acceleration at the base was equal to 0.32g and thus at the crest PGA cr = 1.5 PGA b equal to 0.48.Therefore, the seismic acceleration on the slope at the detachment point was calculated equal to 0.45 g.
Nat. Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2017-29,2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 February 2017 c Author(s) 2017.CC-BY 3.0 License.Normal and tangential COR have prevailed in natural hazard mitigation design via computer simulation due to their simplicity.Values for the coefficients of restitution are acquired from values recommended in the literature (Azzoni et al. 1995; Heidenreich 2004; Richards et al. 2001, RocScience, 2004).Those are mainly related to the surface material type and originate from experience, experimental studies or back analysis of previous rockfall events.This erroneously implies that coefficients of restitution are material constants.However, COR values depend on several parameters that cannot be easily assessed.Moreover, the values suggested by different authors vary considerably and are sometimes contradictory.

Figure 11 and
Figure 11 and seems to describe consistently, but on the unconservative side, the trend and the values acquired by the aforementioned analysis and assumptions.
Nat. Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2017-29,2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 February 2017 c Author(s) 2017.CC-BY 3.0 License.The rockfall trajectory model Rockyfor3D (Dorren, 2012) has been used in order to validate the encountered trajectory and determine the reach probability of the falling rock (from the specific source area) on the impacted house.The 3D analysis was based on the down-scaled 2 m resolution Digital Elevation Model (DEM) that was generated from the 10 cm DSM.The terrain features such as low vegetation (e.g.bushes) and the trees were removed from the DEM as they affected the rock's path downhill.The following raster maps were developed for the 3D analysis: a) rock density of rockfall source, b) height, width, length and shape of block, c) slope surface roughness and d) soil type on the slope, which is directly linked with the normal coefficient of restitution, n COR .The slope roughness was modeled using the mean obstacle height (MOH), which is the typical height of an obstacle that the falling block encounters on the slope at a possibility percentage of 70%, 20% and 10% of the trajectories (according to the suggested procedure in Rockyfor3D).No vegetation was considered in the analysis, which favours a longer trajectory.The parameters considered in the 3D analysis for the different formations are summarised in Table4.The spatial occurrence of each soil type is shown in Figure11and the assigned values of n COR are according to the Rockyfor3D manual.The values for soil type 4.1 are slightly different from soil type 4 Nat. Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2017-29,2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 February 2017 c Author(s) 2017.CC-BY 3.0 License.paths are closer to the actual trajectories compared to RocFall software.The reach probability of the falling rocks, initiating from the source point, is shown in Figure 12.
the factors that control lateral dispersion are classified in three groups: macrotopography factors, factors related to the overall slope geometry; micro-topography factors controlled by the slope local roughness; and dynamic factors, associated with the interaction between slope features and block dynamics during bouncing and rolling.Assessing the results of an experimental investigation, Azzoni and de Freitas (1995) commented that the dispersion is generally in the range of 10% to 20%, regardless of the length of the slope and that steeper slopes present smaller dispersion.Agliardi and Crosta (2003) calculated lateral dispersion to be up to 34%, via high-resolution numerical models on natural rough and geometrically complex slopes.
line with trends described byAsteriou & Tsiambaos (2016), but the deviation of the actual trajectory is significantly lower.This can be attributed to the different conditions (i.e.block shape, slope material, slope roughness, incident velocity and angle, and scale) between the experimental program conducted byAsteriou & Tsiambaos (2016) and the Ponti rockfall event.UAV-enabled reconnaissance was successfully used for the identification of the origin of the detached rock, the rockfall trajectory and the impact points on the slope, emphasizing on the motion types of the trajectory (rolling and bouncing sections).A drone with an ultra-high definition (UHD) camera was deployed to reach the inaccessible, steep and vegetated uphill terrain.A high-resolution orthophoto of the rockfall trajectory and a 10 cm DSM was prepared, which formed the basis for an analytical 2D kinematic analysis and a comparison with the outcomes of 2D and 3D rockfall analysis software.The initial velocity of the detached rock was estimated based on site conditions and amplification of the ground acceleration due to topography.It was found that the Nat.Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2017-29,2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 February 2017 c Author(s) 2017.CC-BY 3.0 License.estimation of the initial velocity of the blocks plays a significant role in the accurate re-production of the rockfall trajectory.

Figure 2 .
Figure 2. Impact of rock on house in Ponti, Lefkada, Greece.

Figure 3 .
Figure 3. Top view orthophoto denoting rolling section, bouncing positions and closeups of impact points.

Figure 6 .
Figure 6.Plan view and cross section along block's path (units in m); 2D rockfall trajectory analysis results are plotted with green and blue line

Figure 7 :
Figure 7 : Out of plane geometry

Figure 13 .Figure 14 .
Figure 13.3D trajectory analysis (from RockyFor3D analysis) addition, noise filtering and smoothing processing were implemented to reduce the effect of construction elements and vegetation in the final rasterized model.Note that Nat.Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2017-29,2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 February 2017 c Author(s) 2017.CC-BY 3.0 License.