Evaluation of rockfalls at a historical settlement area in the Ihlara valley (Cappadocia, Turkey) using different methods

Rockfalls are one of the most dangerous natural events in hilly terrains. This study presents the results of an investigation program to analyze the possibility of a rockfall from a slope to nearby residential buildings in a historical settlement area. Various rockfall analysis techniques were implemented in the study for this purpose. The kinematical analysis revealed the potential of different structurally controlled modes of failure in the slope, especially wedge type and block toppling were the most significant ones. Finite element analysis suggested a stable slope considering the safety factor of 2.19 for the existing geological and geotechnical conditions of the studied slope case. A possible rockfall trajectory was determined and located as an input in the 2D rockfall program based on the field measurements. Different shapes and sizes of blocks were used in the rigid body model for a more realistic numerical simulation of rockfall events. According to the 2D model results, there was no danger of rockfall for the investigated downslope buildings. However, to stay on a safe side, a suitable control measure with a specified dimension was proposed to manage rockfalls in the study area.


Introduction
Rockfalls are among the most dangerous natural events in mountainous terrains, road cuts, quarry faces, and coastal cliffs (Spadari et al. 2012). Rockfall is defined as "the free-falling or precipitous movement of newly detached segments of the bedrock of any size from a cliff or other very steep slope" (UN Glossary 2020). They occur when a rock fragment, isolated rock, or boulder is detached from a steep rock wall, cliff, or slope and then travels some distance by free-falling, bouncing, rolling, and sliding ( Fig. 1-a). Since they occur suddenly, it is not possible to forecast rock falls in advance. Therefore, they pose a severe threat to humans, buildings, and groundworks, particularly if they are not accurately documented for the risk they source ( Fig. 1-b). It is necessary to assess the danger induced by rockfall to secure endangered human lives, residential areas, and infrastructure in closer vicinity (Dorren 2003). Favorable geology and climate are the main factors controlling rockfall occurrence. They include intact behavior of the bedrock, fractures in the rock mass, exposure to weathering, presence of water, freezing-thawing processes, root-wedging, and exterior forces (Attewell and Farmer 1976;Bozollo and Pamini 1986;Giani 1992;Dorren et al. 2006;Binal and Ercanoğlu 2010;Asteriou et al. 2012). Some biological, mechanical, and environmental actions (i.e., seismic activity, ground vibrations, pore-water pressure changes due to heavy rainfall, erosion of surrounding soil during heavy storms, root-growth, or leverage by roots moving in strong winds) can also trigger rockfalls (Hoek 2007).
In this study, an investigation program was undertaken to figure out rock falls to houses in a historical settlement area due to a nearby rock slope (Fig. 2). The study area includes the historical Belisirma village (Guzelyurt-Aksaray), which belongs to a special protection area covering 14 km length Ihlara valley and its surroundings. The site has only been studied in terms of geotechnical and rock mechanics perspectives by a few researchers. Binal (1996) investigated the instability mechanisms observed in volcano-sedimentary rocks, especially in the Kizilkaya ignimbrite in the valley. According to the kinematic analysis results, the study found that the toppling failure was the most probable mechanism in a single block with an eroded base and the system of several blocks in the Kizilkaya ignimbrites. The secondary toppling occurs when an individual block's base angle exceeds 8 o and should exceed 11 o due to differential settlement in the block groups. Sarı and Çömlekçiler (2007) characterized the Kizilkaya ignimbrite using the rock mass rating (RMR) (Bieniawski 1989) classification system. They described the rock mass as "good" with a total score of 62. Tunusluoglu and Zorlu (2009) investigated the rockfall hazard at one of the natural and historical monuments in the Cappadocia region. In the fieldwork, they described the location and dimension of the moved-in-place blocks having the potential to fall. After determining the necessary input parameters, they analyzed the threat caused by these blocks using rockfall software. A potential risk map showing the areas under rockfall danger around the castle was prepared in the study. Sari (2009) developed a stochastic model to estimate the Kizilkaya ignimbrite properties using the Monte Carlo simulation method. He defined the best-fitting frequency distributions for each input parameter of the generalized Hoek-Brown criterion during the simulation process. He found that the strength and deformability parameters of the Kizilkaya ignimbrite could be approximated by asymptotic distributions skewing to larger values like intact rocks. Zorlu et al. (2011) considered the effect of the landforms on the rockfall hazard assessment in the Cappadocia region. They found that differences in the durability between the horizontal layers of the slope were the significant causes of severe rockfall events. In the region, the uppermost levels were composed of wellcemented limestones and welded ignimbrites, whereas the lower parts and slopes were composed of soft tuffs and non-welded ignimbrites. Taşpınar (2015) developed various twodimensional (2D) computer models to assess rockfall hazards in the different parts of the Ihlara valley. She found that it was likely to observe different falling scenarios depending on the block sizes in the Kizilkaya ignimbrite. Ozturk et al. (2019) proposed a low-cost and useful approach for determining orientation (strike and dip) of fractures in a computerized environment using mobile phones and photogrammetric methods. They applied the proposed methodology to a part of the Ihlara valley where the rock units' columnar structure did not allow the discontinuity measurements with a conventional surveying method. A recent numerical study by Sari (2021) in a different part of the Ihlara valley has shown that the valley's cliffs were highly susceptible to secondary toppling failure. In this process, the soft basal rock has been eroded due to water and wind erosion over time (i.e., undercutting), while the block falls from above the hard ignimbritic rock were observed as a result of differential settlement and deformation at the base. Nearly vertical orthogonal joint sets in the upper blocky rock mass have a primary role in promoting the valley's rockfall activity.
It is possible to observe different failure mechanisms in a rock slope with the same geological features depending on the joints' geometric properties and the slope face. To better reflect actual rockfall behavior, it requires a careful characterization of parameters in the field for numerical modeling. Firstly, discontinuity measurements were performed on the rock mass source leading to rockfall on the northeast of the buildings. Different structurally controlled failure modes were investigated for the kinematical analysis depending on the orientation of slope and friction angle of the discontinuities. Secondly, the minimum safety factor was searched under the geological and geotechnical conditions of the slope using the finite element method (FEM). Lastly, different shapes and sizes of blocks, possibly fallen from the slope due to weak structural planes, were simulated using a two-dimensional (2D) rockfall analysis program. Maximum runout distances, bounce height, and total kinetic energy of the rock blocks were determined along the assumed trajectory. Based on the rockfall danger zone defined by the 2D model, a retention wall with a specified location and size was implemented to secure the houses in the study area.

Study Area
The historical Belisırma village belonging to Guzelyurt town is located 35 km southeast of Aksaray city. This village is part of the touristic Ihlara valley, which is well-known due to its historical and cultural heritage and natural beauty. The length of the valley is 14 km with a 52 km 2 area (Şimşek 1997). Hasan-Melendiz mountain range, which lies along the E-W direction and the highest altitude belongs to Hasandağ (3257 m), is bordered by flat-lying plains with an average altitude of between 1100 and 1400 m. The most significant relief structure in the region is the drainage system of the Melendiz river on the north flank of the Hasandağ (Doğan et al. 2019). This system cuts thick ignimbrite-intercalated fluviolacustrine sediments. As a result, canyon-like deep valleys were generally formed in the region. The study area's stratigraphy is mainly composed of late Miocene-Quaternary age young volcano-sedimentary units, including the pyroclastic rocks of Hasandagi ashes, Selime tuff, and Kizilkaya ignimbrite.

Hasandağı Ashes (Th)
Hadandağı ashes were formed by the sequential deposition of lacustrine sediments and volcanic products in the area. Rock types in this formation showed vertical and horizontal transitions. It has a flat topography with a gentle slope. Since it filled the previous surface topography, Hasandağı ashes offer highly variable thicknesses in the region. In general, this formation consists of grayish-white ashes and lapilli in which the white-colored vitric ash matrix includes obsidian, pumice, and lava fragments in sizes from coarse-grained sand to 5-6 cm in diameter. Hasandağı ashes cover the Kizilkaya ignimbrite on the west flanks of the Ihlara valley.

Selime Tuff (Ts)
Selime tuff, which covers a considerable part of the study area with its light purple to yellowish-white color, was first named after Beekman (1966). It lies in a narrower band along the valley starting from the Ihlara town and outcrops on a broader domain around the Selime town. Selime tuff is recognized with Neogene-aged fairy chimneys, a spectacular badland landscape for Turkey's Cappadocia region. Fairy chimneys are specific landforms that originated from differential eroding of pyroclastic-flow deposits consisting of a very poorly sorted mixture of volcanic ash or tuff during pumice lithification with rock fragments in fluvio-lacustrine environments (Sarikaya et al. 2015). When it was eroded together with overlying hard Kizilkaya ignimbrite, an impressive mesa capping a badland landscape on the valley's right bank was observed in the region. Selime tuff is rich in basalt, spilite, obsidian, tuffite, pumice, and andesite fragments showing locally color change facies in a light purple to whitish-yellow. In these facies, the radius of the pumice fragments also increases in the tuff.

Kizilkaya Ignimbrite (Tk)
The welded Kizilkaya ignimbrite outcropping of various sizes in the region has the widest distribution, which is mesa hills with a badland landscape. The most beautiful location in which the contact relationship of this mesa unit and the Selime tuff underneath is seen in the north of Kızılkaya village. Beekman (1966) named ignimbrites as Kızılkaya ignimbrite about this area. It is generally whitish-gray, with slightly weathered surfaces turning a pinkish color. The ignimbrite, which has a flat-looking topography, was divided into pieces by the Melendiz river and its tributaries forming mesa hills at various heights. In the Ihlara valley, which is in the form of a deep narrow canyon, Kizilkaya ignimbrite shows pervasive columnar jointing due to cracks developed in sub-vertical direction during the cooling. The thickness of ignimbrite reaches approximately 60 m in the valley. It contains pumice fragments with weakly textured at the upper levels, densely textured at the middle classes, and rough, densely textured at the bottom up to 30 cm in size. A red-colored cooking zone is noticed between Kizilkaya ignimbrite and the underlying Selime tuff. The thickness of the Kizilkaya unit varies between 2-75 m in the study area. Kızılkaya ignimbrite's age was 4.9-5.5 million years found by the K/Ar method in biotites (Ayhan and Papak 1988).

Kinematical analysis
The kinematical analysis is a practical preliminary method to investigate the stability of rock slopes in which the failure is mostly controlled by the predominant discontinuity sets in the rock mass (Hoek and Bray 1991). According to the structural and engineering geology studies related to the field, discontinuity systems can create kinematically different rock mass instabilities in the studied slope. For initiating a rockfall, a block should be formed by two intersecting joints, and it is kinematically capable of dislodging from the slope face. In this process, kinematical analysis is the first step in identifying slope conditions for potential structurally controlled failures, such as plane, wedge, and toppling (Hoek and Bray 1991). Therefore, before doing any rockfall analysis, it is crucial to recognize structural discontinuities with respect to the orientation of the slope face. Fig. 4 shows the typical failure modes and conditions for the kinematic analysis. These may take the form of a plane failure ( Fig. 4-a) in which rock mass contains continuous joints dipping out of the slope face, and they strike parallel to the slope face; a wedge failure ( Fig. 4-b) formed by two intersecting planes of discontinuities; toppling failure  in which rock mass contains fractures steeply dipping into the face; and a circular failure  in the rock fill, very weak/weathered rock or closely jointed rock masses containing many randomly oriented discontinuities (Wyllie and Mah 2004).
Given that discontinuities play a critical role in the kinematical analysis, 53 discontinuity measurements were performed on the hard ignimbritic unit covering the slope ( Fig. 5-a).
Two sub-vertical orthogonal joint sets were delineated from the rock mass forming the slope crest. A third joint set lying horizontally is almost perpendicular to these two sets (Fig. 6).
Sub-vertical joint sets (Set 1 and Set 2) were developed as a result of thermal cracks formed while the cooling of ignimbrite massif took place. Near horizontal joint set (Set 3) was formed from different outflow facies during ignimbritic deposition (Sarı and Çömlekçiler 2007). Table 1 presents some characteristics of the discontinuity sets in the rock mass and other input parameters necessary for the kinematic analysis. Tilting tests were conducted in the field to estimate the friction angle of the discontinuity surfaces ( Fig. 5-a). This test is based on measuring the friction angle between two blocks placed on top of each other in the field directly when the upper one starts to slide on the lower one. It is an in-situ test for better reflecting the actual conditions in the field and considering the size effect. The average of three tilting tests was calculated as φ=45 o .
The rock mass lies on the upper slope of the houses subjected to rockfall (see Fig. 2-a). The analysis considers the SW face of the slope looking towards the houses ( Fig. 2-b).
Accordingly, the slope's dip/dip direction was recorded as 85 o /205 o in the field. For the rock mass in the slope, three joint sets had been identified with a dip/dip direction of joint set 1 (86 o /275 o ), joint set 2 (88 o /194 o ), and bedding set 3 (2 o /201 o ). Goodman (1989) and Hoek and Bray (1991) described the method employed in kinematic analysis using the lower hemisphere stereonet projection. The data in Table 1 was used in Dips v.8.0 software (Rocscience Inc. 2020a) for analyzing modes of the planar, wedge, and toppling failures.

Planar Failure
Compared to other modes of failure, planar failure usually occurs seldom in nature. In the kinematic analysis of a planar instability, the critical zone is defined as the area encircled by the outside of the pole friction cone and inside the daylight envelope. The poles located within this zone are likely to slide. Accordingly, the slip plane angle must be greater than the friction angle and smaller than the slope's inclination. Also, the direction of the slip plane and the slope must be approximately parallel to each other. In other words, the difference between the dip direction of the slip plane and the slope should be at most ±20° (Goodman 1989). Fig. 7 presents suitable conditions for a planar sliding on the stereonet. Each discontinuity that falls into the red shaded region has a planar sliding potential. Accordingly, out of 53 discontinuity surfaces measured in the field, only 3 of them are located in the critical region, and these three discontinuity planes are members of Set 2. As a result, the potential for planar sliding on this slope has been calculated as 5.66 %.

Wedge Failure
In case two intersecting planes form a rock wedge, and if it slides along the intersection direction, it is called a wedge failure. In the stereonet, the area inside the friction cone plane and outside the slope's plane is the critical zone. Wedge instability is characterized by the joint planes intersected in the critical zone (Wyllie and Mah 2004). The stereoplot program calculates the trend and the plunge of all intersecting planes. At least two different pole concentrations must be observed on the stereonet for a wedge-type sliding. The red shaded region given in Fig. 8 is considered as the critical zone. The intersections falling in this region have the potential of wedge-type failure. Accordingly, out of the 1377 intersections formed by the three sets on the slope, only 229 were plotted in the critical zone. As a result, the potential for wedge-type failure on this slope was found to be 16.63%.

Toppling Failure
Toppling is a common failure type in slopes with sub-vertical joints. In this failure type, rock columns or blocks have to rotate along a fixed base. Toppling instability is classified broadly as direct or block toppling and flexural toppling. In block toppling, rock columns are formed by one set of joints steeply dipping into the slope face, and the length of broadly spaced cross joints defines the column height. However, in flexural toppling, the rock columns separated by the continuous joints can bend forward before instability (Wyllie and Mah 2004). The slip limit plane in the steorenet plot defines the critical zone in flexural toppling. The deduction of friction angle from the slope angle gives the dip of the slip limit plane, and its dip direction is the same as the dip direction of the slope face. The other condition for toppling is that the difference between the slope's dip direction and the discontinuity sets needs to be within ±30° (Goodman 1989). The pole plots within the critical zone in Fig. 9 illustrate the block toppling. Accordingly, the block failure was observed in only 112 of 1377 intersections defined on the studied slope. That is, the direct block failure potential on this slope was found to be 8.13%. The potential for oblique (lateral) block failure on this slope was calculated as 10.9%. One of the most critical factors controlling block toppling was the presence of a basal plane on which the blocks would topple. The majority of the discontinuities (58.82%) in Set 3 form the basal plane where the slip would occur. When the flexural type of toppling failure was considered in Fig. 10, only four discontinuity planes showed a potential for this type of failure with a percentage of 7.55 among the evaluated discontinuities.
The kinematic analysis results given in Table 2 indicated that wedge and toppling failures would be the leading cause of block instability in the slope. As shown in Fig. 11, some of the wedges observed on the downslope show consistency with the kinematic analysis results.

Shear Strength of Discontinuities
The rock mass is a geologic structure divided by visible defects into intact rock blocks of various sizes. In engineering excavation built in/on rock structures very close to the surface, where the effective stresses are very low, the failure and deformation are mainly governed by the shearing resistance of the discontinuity surfaces (Barton and Choubey 1977). In planar discontinuity surfaces, since there was not enough resisting force against the shearing, it could damage the geologic structure with a minimum stress change. Filling materials like clay, silt, calcite can be easily worn out, and the discontinuity surface is now in the state of residual shear strength (Barton et al. 1974). On the other hand, roughness and undulation have more pronounced effects on the shear strength of natural joint planes. In general, the shear strength of rock joints increases as the surface roughness increases, and it is essential for the stability of rock excavation (Barton 1973). In this study, the non-linear Barton-Bandis (BB) shear failure criterion (Barton et al. 1974;Barton and Choubey 1977) is used to determine the shear strength of rock joints due to its simplicity and easiness. Barton et al. (1974) had studied the natural joint planes and proposed the following equation: (1) where τp and σn are the shear and normal stresses, respectively; JRC is the joint roughness coefficient; JCS is the joint wall compressive strength; and φb is the basic friction angle of the unweathered rock surface.
Later, Barton and Choubey (1977) changed this formula by using 130 direct shear box experimental results of joint surfaces in weathered rocks.
Where φr is the residual friction angle, Barton and Choubey (1977) determined that the friction angle can now be estimated as: Here, r is the Schmidt value of the wet and weathered surface, R is the Schmidt value of the non-weathered fresh surface. It can be seen that the basic friction angle (φb) plays a crucial role in predicting the shear strength of the discontinuities. The basic friction angle characterizes fresh surfaces. The basic friction angle has been investigated by Barton and Choubey (1977) for different rock types, and they recommended 25°-30° for sedimentary rocks, and 30°-35 ° for magmatic and metamorphic rocks. The basic friction angle can be calculated for fresh, smooth surfaces using the tilt test or the direct shear box test in the laboratory (Alejano et al., 2012).
The joint roughness coefficient (JRC) was extracted from the joint surfaces of hard ignimbrite in the slope crest. As seen in Fig. 12-a, a wooden ruler 1 m in length was placed on the joint surface to measure asperity amplitude using a vernier caliper. The average asperity amplitude value (a=22.2 mm) measured on the discontinuity surfaces was located on the graph in Fig. 13 to find a JRC value of 10. The joint wall compressive strength (JCS) was obtained from the Schmidt hammer test. This test is a practical indirect method for predicting the compressive strength of joint planes in the field. The hammer was applied perpendicular to the discontinuity surfaces during the test. Accordingly, the JCS value of joint surfaces in ignimbrite rock was found using the graph in Fig. 14, which was developed by Barton and Choubey (1977) from the equation proposed by Deere and Miller (1966). The JCS value was found to 55 MPa for the unit weight γd=21.95 kN/m 3 and Schmidt hardness value, r=38.2, measured on the weathered joint surfaces in the field.
In this study, the basic friction angle was obtained by performing a simple tilt test where a manually-operated wooden table was utilized in the experiments. The most widely used method for this test was early proposed by Stimpson (1981), and later this method was revised by Alejano et al. (2012). As shown in Fig. 15, two disc-shaped specimens prepared for the Brazilian tensile test were used in the experiments, which had more conservative values than other sample shapes. Eq. (4) was used for the conversion of β to the basic friction angle.
Here β is the slope of the tilt assembly at the start time of the slide.
The residual friction angle was calculated according to Eq. (3) using the basic friction angle obtained from the tilt test. Accordingly, the peak friction angle for discontinuity planes of hard ignimbrite was determined as φp=52.6 degrees. The BB failure envelope of the joints based on Eq. (2) is depicted in Fig. 16. Using the least-square regression technique, the nonlinear BB failure envelope was then linearized in the range of normal stresses effective in the rock mass to implement Mohr-Coulomb's linear model for the numerical analysis in the next section.

Finite Element Analysis
The stability of the slopes is one of the most studied and widespread geotechnical analyses in engineering problems. It is essential to understand the processes and mechanism driving the instability in an unstable rock slope to evaluate the potential hazard associated with it (Eberhardt et al. 2002). There are several numerical methods used in solving such instability problems. The finite element method (FEM) is a numerical method mostly utilized in the failure behavior of continuous rock masses under effective in-situ stress conditions. Due to the capabilities of FEM, many studies in recent years were conducted using this method in the modeling of rock mechanics related problems 2008;Fu and Liao 2010;Azami et al. 2012;Pain et al. 2014;Alemdag et al. 2019;Sari 2019). Due to the simplicity of its coding programs, less computation power and time, this method was employed to perform the numerical analysis in this study. This method can analyze the rock mass's failure propagation within defined structural defects such as joints, faults, fissures, schistosity, foliation, folding, and beddings.
It is a well-known fact that the stability of rock slopes is significantly affected by discontinuities in the rock mass since the discontinuities control the rock mass's overall structure and mechanical properties (Hoek et al. 2002). The joint intersections are the potential locations of high stress, deformation counters, shear damage, and instability (Mas Ivars, 2010). Besides, the effect of external forces such as groundwater and earthquake loads can be easily incorporated into built FEM models. It is also probable to define different rock material and rock mass constitutive models for the materials forming the slope. At the end of the analysis, a unique safety factor representing the slope's stability condition can be calculated numerically using the shear strength reduction (SSR) method (Sari 2019). A critical strength reduction factor (SRF) can be easily predicted with total displacement and major shear strain counters of the failed slope. Using this method, the possible location of the unstable rock blocks can be defined.
A reliable slope stability analysis needs to assume the right failure criteria with accurate input parameters in jointed rock masses (Sari 2019). An experimental study was undertaken to obtain the unit weight, uniaxial compressive strength (UCS), Brazilian tensile strength (BTS), and Schmidt hardness value of rock materials based on ISRM (2007). UCS and BTS tests were conducted using NX (diameter 54.00 mm) core sizes (Fig. 17). It was unlikely to have high-quality standard core specimens from the coarse-grained tuff blocks since they were exposed to intense weathering. The test results were imported into RocData v.5.0 software (Rocscience Inc. 2020b). The best-fitting generalized Hoek-Brown (HB) failure criterion envelopes are shown in Fig. 18 and Fig. 19 for the intact rock samples. Hoek et al. (2002) developed the HB failure criterion for upscaling the strength envelopes obtained from the rock material to the field scale rock masses. HB failure criterion provides a dimensionless equation given in Eq. (5) to connect geological observations with calculations.
Here, σ1 and σ3 are the major and minor principal stresses, respectively. To estimate the GSI plays a critical role in the HB failure criterion, and a visual chart is provided for a qualitative description of the rock mass structural elements and discontinuity surface conditions. The GSI values of the studied rock masses were directly determined in the field following the GSI chart provided by Marinos and Hoek (2000), as seen in Fig. 20. According to the field observations, the slope's rock structure had been separated into two distinct zones.
The upper welded zone was named "hard ignimbrite" while the lower unwelded zone bearing rock carved structures was named "soft ignimbrite". The rock formation covering the downslope surface between the slope and the dwellings was called "tuff". The strength envelopes of three rock masses showed significant variations as presented in Fig. 21, both for principal and shear-normal stress spaces. The resulting rock mass properties in Table 3 were employed in the 2D finite element analysis program RS2 v.2019 (Rocscience Inc. observations as seen in the exaggerated inset picture. Furthermore, a relatively high value of SRF= 2.19 was found by the SSR analysis. This value indicates that the overall slope failure would occur in a long period. For an immediate collapse, the calculated SRF value should be 1.0 or below. In the case of the total displacement counters on the slope face that were considered in Fig. 24, one can easily discern that the major movements would take place on the slope crest where the hard ignimbrite was separated into many intersecting blocks by discontunity planes. It clearly shows unstable blocks on the studied slope as a source area of potential rockfalls.

Rockfall Analysis
Rockfall analysis is a function of the source area's location, the geometry, and geomechanical properties of both the block and the slope. In theory, if the initial conditions, the slope geometry, and the energy loss at impact or by rolling can be determined initially, it is probable to measure the position and velocity of a block at any time. In practice, however, it is often impossible to accurately determine the exact location of source points of the rock blocks and their size, shape, and geomechanical properties. Besides, the surface formation's geometrical and mechanical properties can change considerably along a slope (Agliardi and Crosta 2003). Slope geometry, slope roughness, static and dynamic friction, rolling resistance, the density of rocks, and the restitution coefficients are the most important parameters defining particular rockfall trajectories. Among these physical properties, the coefficient of restitution plays a significant role in locating the exact rockfall trajectory (Ansari et al., 2015). It is also the most important input parameter quantifying the energy absorption upon an impact (Ji et al., 2019). It depends on the rock type covering the surface, the vegetation cover, the shape and size of the falling block, and the slope's physical properties (Dorren et al., 2006;Dadashzade et al., 2014). Therefore, a reliable estimation of the coefficient of restitution is of profound importance in rockfall prediction and for designing countermeasures against rockfall (Chau et al. 2002).
In the lumped-mass impact models, the normal and tangential coefficients of restitution (Rn and Rt, respectively) are employed to replace the lack of physics required by the basic equations. The normal coefficient of restitution, Rn, is the amount of energy loss during a falling body's impact in the slope's normal direction. The tangential coefficient of restitution, Rt, is the amount of resisting forces caused due to moving parallel to the slope. Site-specific restitution coefficients can be estimated directly from laboratory or field tests, back analysis of falling blocks, or using theoretical estimation methods (Bozzolo and Pamini 1986;Kobayashi et al. 1990;Evans and Hungr 1993;Budetta and Santo 1994;Robotham et al. 1995;Chau et al. 2002). The most common method is to estimate them from the back analysis of measured rock paths and end locations during the field tests of rockfall trajectory (Bar et al. 2016).
In this study, the coefficients of restitution were determined from a back analysis in RocFall v.5.0 software (Rocscience Inc. 2020d). For this purpose, the location and size of the two fallen blocks observed in the study area were measured during field visits. As seen in Fig.   25, one of the blocks was 10 m away from the slope toe, and the second one was 8 m away.
By considering the weights of the blocks and their distances from the slope toe, the corresponding normal (Rn) and tangent (Rt) restitution coefficients of tuff forming the downslope surface were determined Rn=0.25 and Rt= 0.60 according to back analysis results given in Fig. 26. In this process, the required parameters were adjusted incrementally until modeled runout paths showed a similar spatial distribution to the actual block locations.
For the houses located at downslope, two dimensional (2D) rockfall analysis was executed along the AA-section line on the slope using RocFall v.5.0 software (Fig. 27). The program simulates the rockfall trajectory of the potentially unstable rock blocks by calculating runout distance, bounce height, and total kinetic energy for the given 2D profile. RocFall program additionally allows the user to perform both lumped-mass and rigid body rockfall analysis methods. In the lumped-mass models, a block is considered as a tiny particle with a mass.
Falling rock is defined as a point mass; the actual size and shape are neglected even though they would otherwise affect its trajectory (Rocscience Inc. 2020d). In lumped-mass impact models, the mass of fall rock is only used to calculate energies, and it does not affect the overall body trajectory. Only the sliding motion is replicated while the rotation is characterized by a zero friction angle (Basson 2012). On the other hand, in rigid-body impact models, the essence of fall rock behavior is captured by using the laws of motion and kinematics. It is assumed an instant contact period and that the contact area between the impacting objects is insignificant. The fallen rock's shape and size with four types of movement (i.e., fall, slide, bounce, and roll) are considered in rigid-body impact models (Basson 2012). Modeling all four modes of motion is necessary for a realistic, more consistent risk-based rockfall hazard study (Vick et al., 2019).
In this study, the rigid body impact module is used in the rockfall analysis to obtain more realistic outcomes. Discrete rockfall boulders are modeled with arbitrary and random shapes and size/mass. Various smooth and polygon shapes of sphere, square, triangle, egg, oval, ellipse, rhombus, rectangle, pentagon, hexagon, and octagon were specified for the potentially unstable blocks. The size/mass of a block had been selected based on the varying in-situ block sizes and the average density of hard ignimbrite. The densities of the ignimbrite blocks were determined on the core samples in the laboratory. Accordingly, the dry density of ignimbrites was measured as 2.17 gr/cm 3 , and the water-saturated density was 2.24 gr/cm 3 . In rockfall analysis, generally, considering the worst condition, water-saturated density is more applicable in determining the block weights. Thus, the varying block sizes of 10 kg, 100 kg, 1,000 kg, 5,000 kg, 10,000 kg, and 20,000 kg were executed in rockfall analysis. A total of 1,200 simulations were realized using the blocks in different shapes and sizes to cover possible rockfall trajectory scenarios. Due to the slope geometry and block orientation restrictions, initial velocities of the blocks in vertical and horizontal directions were taken to be 0.5 m/s in the model. Other input parameters selected in the analysis are given in Table 4. With difficulty in data availability, some of these input parameters were estimated from the typical values used in the RocFall program's case histories.
The rockfall analysis results for the AA-section profile are presented in Fig. 28 and Table 5.
According to the results of the 2D analysis, blocks falling from the crest of 18 m high rock slope reached a maximum runout distance of 21.13 m. Most of the blocks were stopped after bouncing and rolling at a small distance from the slope toe. It appears that none of the blocks modeled by the program can reach the houses located at a distance of approximately 55 m (see Fig. 28). In terms of movement mechanism, a block that can detach from the upper elevations of the steep slope primarily exhibits a free-fall at the beginning. Later, the block hits and bounces on the horizontal surface near the slope toe. After very low bouncing heights (60-75 cm), the block rolls down along the slope surface, loses its kinetic energy, and stops when it reaches the flat area. The observations made in the field also confirmed the simulated paths. Previously fallen rock blocks in various sizes were also very close to simulated locations in the model.
The number of blocks and their runout distances, the bouncing heights, and total kinetic energies on the slope profile was also determined during the rockfall analysis. To design a practical mitigation structure, the endangered rockfall zone and the impact energy of falling rocks must be accurately predicted (Yan et al. 2020). Accordingly, most of the blocks were rolled not more than 10 m horizontal distance on the slope profile. Out of 1,200 simulated rock boulders, only twelve could reach up to 20 m distance (Fig. 29). When the bouncing heights given in Fig. 30 were evaluated, it was seen that the blocks jumped up to 8 m mean bounce height at first impact, but as they move away from the slope toe, their bounce heights decrease suddenly. It was found that the block that could travel the farthest distance would jump only 1.1 m. When the total kinetic energy distribution given in Fig. 31 was examined, similarly, the blocks had mean kinetic energy (244 kJ) with the effect of the free fall at the beginning, but they lose this energy as a result of friction during the rolling along the lower section of the slope, and this value decreases up to 12 kJ at rest.

Mitigation Against Rockfalls
For the planning of rockfall mitigation measures such as protective barriers, it is necessary to determine the maximum kinetic energies (for the bearing capacity of the protective barrier to hold the rock block) and the largest jump heights (for the barrier height) of the potential rock blocks. Also, the evaluation of the runout distance for estimating the rockfall hazard zone is another necessity. The 2D rockfall analysis showed that ignimbrite blocks of different sizes and shapes detached from the steep slope with a height of about 18 m could not reach the dwellings subject to investigation. However, there is always a risk of falling rocks for dwellings in the study area, and necessary precautions should be taken to mitigate this danger.
Typically, rockfall events can be mitigated either by active or passive measures. In passive mitigation, the rockfall can still occur, but an effort should be made to prevent the consequence. The negative effects of rockfall events can be diminished by employing rockfall catchment fences, drape nets, diversion dams, rock sheds, and forest belts in the runout or deposition zones. However, in the active mitigation, rockfall event occurrence is prevented from the source zone by rock bolting, slope retention systems, shotcrete, altering slope geometry, dewatering the slope, and revegetation. Several researchers (Ritchie 1963;Pierson et al. 1990;Pantelidis 2009 andBar et al. 2016) have outlined some guidelines for designing passive measures.
It is suggested to construct a protection dam on the slope. It is built from locally available stones suitable for the natural texture of the region. Fig. 32 illustrates the optimum location of the dam with specifications of a width of 1.5 m, a height of 2 m, a length of 15 m, and resilient to a total kinetic energy of 500 kJ. Although the jump heights on the slope profile would not reach significant levels (generally <50 cm), the size of in-situ blocks that can detach from the investigated ignimbrite rock is likely to be greater than 2 m when the discontinuity spacing values are considered in Table 1. For this reason, the height of the rock holding structure should be at least 2 m to prevent the block from passing over the wall along the long axis during the rolling (even if it does not bounce), as seen in Fig. 32.

Conclusions
This study investigated the possibility of rock falls to the buildings in a historical settlement area from a slope. According to the kinematical, finite element, and rockfall analyses, the following conclusions were driven from the study.
• As a result of the investigations and measurements conducted in the field, several discontinuity systems can create kinematically different rock mass instabilities in the said slope. Accordingly, there is a potential for wedge and block toppling type discontinuity-controlled instability in the examined slope.
• Wedge type shear failure shows the highest potential (16.63%), while block toppling has a possibility of 8.13%.
• In the FEM analysis, the overall stability of the slope was found as relatively high with SRF=2.19. This value indicates that it was impossible to observe a complete failure of the slope in a short period.
• For more accurate results, various sizes and shapes of blocks were examined in the 2D rockfall analysis. The results clearly showed that the detached ignimbrite blocks from the steep slope would not reach the houses on the downslope. However, necessary preventive measures should be implemented to cease any negative effects in the future.
• Accordingly, it was recommended to build a protection dam made of locally available stones suitable for natural texture, with specifications in 2 m high, 1.5 m wide, 15 m long, and resistant to 500 kJ total kinetic energy.
• It has been observed that the steel fence/net (rockfall barriers or catch fences), which is widely used in rockfall prevention projects, is not suitable for the studied rockfall danger area due to the gentle slope and the natural texture. • The proposed dam is only for protecting the houses located in parcel no 509 from the danger of rockfall, but not for the other houses in the surroundings.