A Comprehensive Approach to Quantitative Risk Assessment of Rockfalls on Buildings Using 3D Model of Rockfall Runout

Abstract

Rockfalls are incidents of nature that take place when rocks or boulders break from a steep slope and fall to the ground. They can pose considerable threats to buildings placed in high-risk zones. Despite the fact that the impact of a rockfall on a building can cause structural and non-structural damage, few studies have been undertaken to investigate the danger associated with this event. Most of these studies indicated that the risk resulting from rockfall hazards is hard to forecast and assess. A comprehensive quantitative risk assessment approach for rockfalls on buildings is developed and described in this paper and applied for the Mtein village in Mount Lebanon. This method employs a 3D model to simulate the rockfall trajectories using a combination of digital elevation data, field surveys, and orthorectified aerial photographs. The spatial and temporal probability of rockfalls were evaluated using the analysis of historical data in two triggering-factor scenarios: earthquake and precipitation. The findings show that, during the period of 1472 years between the years 551 (the first observed large earthquake in Lebanon) and the current year of the study (2023), the temporal probability will potentially be equal to 0.002 and 0.105 in the cases of earthquake- and rainfall-triggered rockfalls, respectively, while the maximal damage values are expected to be 232 USD and 10,511 USD per year, respectively. The end result is a final map presenting the risk values assigned to each building that could be damaged by rockfalls.

1. Introduction

1.1. Definition and Generalities

Rockfalls (RFs) are geological hazards that occur when rocks break free from a cliff or slope and fall downward. They are complex events that involve rocks detaching from a slope and then falling, rolling, and bouncing down the slope. The main thing that makes RFs move is the interaction between the rocks and the slope, but there are other factors that can affect their motion, such as natural and man-made obstacles and the way in which the rocks interact with each other. A multitude of reasons, including severe rain, earthquakes, and human activity, can cause them.

For instance, and because heavy rain events can trigger RFs, climate change could represent an additional risk to people, buildings, and infrastructure in terms of RFs [1] Earthquakes can cause RFs, as the ground shaking caused by an earthquake can destabilize slopes and trigger an RF. Earthquakes can also cause cracks in rocks, which can make them more likely to fall [2]. Human activities can trigger RFs in a variety of ways, including mining and quarrying, construction, deforestation and road building. All of these activities can destabilize slopes by removing rock and soil, changing the weight distribution on the slope, or creating cracks in the rock. This can make the slope more susceptible to collapse and trigger an RF [3]. Thus, RFs can happen quickly and unexpectedly with no previous signs, and it is difficult to understand how rocks will behave during an RF because of different unexpected changes in all of these factors [4]. RFs can be catastrophic as they can cause extensive damage to buildings and infrastructure [5,6] and represent a menace to people’s lives. The falling rocks contain an enormous amount of energy, which might cause the structure (buildings/infrastructure) to fall apart or impose significant damage [7] and also cause injury/death to persons affected by their impacts [8]. To limit the potential harm to structures, developers and land planning authorities must have an in-depth understanding of the causes and trajectories of RF accidents [9]. This knowledge allows them to determine places at high risk of RF, allowing them to make informed judgments about permissible and non-accepted building placements. When choosing and planning where to build it is important to identify potential RF hazards, as well as their impact areas, with the help of the research on rock slope stability [10,11]. Rock slope stability is the ability of the slope of a cliff to resist the RF hazard of a specific rock at a particular location. Thus, a stable rock slope is less likely to experience RFs; however, even a stable rock slope can experience rockslide if it is subjected to a sufficiently large or intense triggering event, such as an earthquake, heavy rain events, or freeze–thaw cycle.

1.2. Risk Assessment Methods

To assess the risk of RFs, scientists must study the features of the rock slope, such as the type of rock, the slope angle of the cliff, and the environmental conditions surrounding the studied cliffs and their rocks that are prone to RFs. There are many ways to conduct this, including field observations, laboratory tests, computer models, and the analysis of satellite images [12]. By analyzing the data that are obtained within a simulation framework, it is possible to forecast and evaluate the probability of future RF incidents and their potential impacts.

The information generated from these simulations can be used by developers and land planning authorities to conduct an RF risk assessment. This latter is then employed to make educated judgments about building placement and to create urban master plans for a region. They can, for example, avoid building structures in high-risk regions or use protective measures such as RF barriers or slope stabilization techniques to reduce the potential harm from RFs [4,13,14]. RF risk assessment is a difficult process due to several factors, including incomplete historical records, the specific RF nature of each location, the difficulty of accurately predicting where RFs will occur, and the lack of knowledge about how much damage RFs will cause [15,16]. In order to address these issues, probabilistic approaches for assessing the risk have been progressively developed through scientific research. There are numerous approaches for assessing the risk of RF. These approaches can be categorized into three general types: qualitative, semi-quantitative, and quantitative. The use of expert judgment to assess the risk is involved in qualitative approaches.

Qualitative approaches are frequently employed for preliminary assessments or assessments with little data. Semi-quantitative methods involve combining expert opinion with quantitative data to determine the hazard and risk of RF. Semi-quantitative approaches are frequently employed for more complex assessments based on more data than those of qualitative methods. Quantitative methods entail the use of mathematical models to quantify the hazard and risk of RF. Quantitative methods are often used to assess complex hazards and need more data than both aforementioned methods. Quantitative methods provide more objective and accurate results than both qualitative and semi-quantitative methods; however, they are also more complex and require more data [17].

In order to apply a specific risk assessment technique, the suitability of all of these techniques must be determined and judged according to the specific context of the study area, such as the spatial scale and complexity of the RF risk, the availability of data, and the budget of the study. The spatial scale of the hazard is the size of the area that is potentially affected by RFs. This can range from a small area at the base of a rock slope (as the accumulation zone) to a large area, such as a valley or an entire mountainside. The spatial scale of the hazard will influence the selection of the most appropriate method of assessment and management. The complexity of the hazard is the degree to which it is influenced by a variety of factor combinations. These factors can include the geology of the rock slope, the slope geometry, the presence of discontinuities in the rock, the groundwater conditions, and the external loading on the slope. The budget available for RF risk assessment and management will also influence the method that is chosen. More sophisticated methods are typically more expensive. Therefore, it is important to choose a method that is appropriate for the budget available and for the level of risk that needs to be managed.

In general, most of the studies made on RFs are limited to hazard assessments. These are evaluations of the potential of an RF to occur which do not take into consideration the consequences of such an event [18,19]. This is due to the difficulty in creating quantitative estimates of risk [20] that arises through the uncertainty about the system being assessed and the uncertainty about the information that is available [16,21].

For, instance, quantitative risk assessments (QRAs) are rarely carried out in areas at risk of RFs [22] and, because they are complex and because there is often a lack of resources in these areas, are site specific [19]. As a result, QRAs are undertaken only in places where there is a historical record of RF events [20,23].

As used in the present study, the quantitative risk assessment (QRA) method will be discussed in more detail, as follows. The QRA is a method for assessing the risk of RF hazards. It was originally developed by Corominas et al. (2005) [23] to assess the risk of landslides [22,24], it was then modified to incorporate the unique features of RF processes [20,23,25,26]. This method allows stakeholders to compare risks across different locations and to make informed decisions about how to manage risk [22]. It provides useful information when comparing the costs and benefits of different management plans in different scenarios toward the end of making optimal decisions and enacting effective policies.

1.3. Quantitative Risk Assessment

Agliardi et al. (2009) [20] has outlined the elements involved in evaluating RF risk. These elements include the likelihood (hazard probability) of RF events, the extent to which an RF will damage structures, the susceptibility of structures to damage (their vulnerability), and the estimated cost of damage.

Thus, to effectively develop RF mitigation measures, it is necessary to analyze the details of the RF risk, including the probabilities of RF incidents, the spatial probability, the intensity of impacts on structures, the vulnerability of the structures, and the total costs associated with different circumstances. This knowledge can help stakeholders and guide their decisions when planning and implementing cost-effective actions that protect people, infrastructure, and property.

QRAs use probabilities to quantify the likelihood of RFs occurring, the presence of structures in their paths, and the vulnerability of those structures to damage as well as the cost of such impact on structures [15].

Risk to be assessed (damage) = hazard probability (occurrence and spatial reach) of a rockfall × probability of element being at a particular point in the RF pathway at the moment of the RF × vulnerability × value of the element at risk

(1)

where hazard is evaluated by combining the annual probability of occurrence, which is typically obtained from RF inventories [27] and from the spatial probability (reach probability) of impact on structures using the simulations of three-dimensional numerical models [20].

The probability of an element being at point X at the moment of an RF depends on a number of factors, including the location of the element relative to the RF source, the type of element, the time of day, and the weather conditions. The probability can be estimated using a variety of methods, including historical data, numerical modeling, and expert judgment.

Furthermore, vulnerability considers the intensity of the RF as well as the resistance against the impact of an RF of a given at-risk element that is present at the impact point.

The value of the element at risk is an important factor that is used when assessing the risk of damage from an RF. This value can be defined as the economic, social, nor environmental value of the damage and can be estimated using a variety of methods, including the market value, replacement cost, social value, or environmental value of the element at risk.

The geology of the area, slope angle, and the presence of fractures or other discontinuities in the rock all influence the likelihood of an RF occurrence. The size and shape of the RF, the distance it travels, and the sort of structures in the vicinity all influence the level of damage. The susceptibility of structures to damage is determined by the structure’s strength and construction, as well as the materials employed in its creation. Finally, the estimated value of damage is computed by multiplying the likelihood of falling rocks by the cost of damage from that size of RF and the vulnerability of the assets.

As RF accidents are uncommon and only occur in certain places, and as QRA are complex and interact with structures in a variety of ways, this method is applied and validated in this paper using a sophisticated modeling approach toward RF processes [22,24].

The following section of the paper will detail the methodology used to assess the RF risk, including the selection of the area of interest, the presentation of data collection and their employed techniques. Section 3 presents and discusses the risk assessment results and analyzes the implications of the findings. The fourth and final part summarize the study principal findings, highlight the study limitations, and make suggestions for future research.

2. Methodology

According to the QRA equation, the methodology for risk calculation consists of calculating the values of the following parameters: reach probabilities; vulnerability and the economic value (cost) of the element being at risk, in this case buildings; and the temporal probability of RF occurrence. Three-dimensional RF simulation was employed to obtain the reach probability and the intensity (which is a part of the vulnerability assessment) of RFs. The 3D simulation will be performed using Rockyfor3D. The simulation requires input data, such as (a) digital elevation model; (b) geological and soil properties, derived from geological and soil maps and from the data obtained from site visits; and (c) RF properties including size, lithology and the spatial distribution of rocks (determined by RF map** based on aerial orthorectified photos).

Historical records of large earthquakes and heavy rainfall events are used to estimate thresholds for initiating RFs. Based on these thresholds, temporal probabilities in each case are calculated.

In line with this and in the present study, two scenarios are taken into consideration, the first considers an earthquake as the triggering factor, and the second considers rainfall.

Furthermore, aerial photos are used to check and evaluate the distribution of buildings and their characteristics, to allow for the determination of their resistance against RFs (this phase represents the second part of the vulnerability assessment) and the cost of repair in case a building is hit by an RF. The methodology used in this study is presented by the flowchart in Figure 1.

2.1. Study Area

The study region within the Mtein village, located in western Mount Lebanon, was chosen because of its elevated RF potential.

The elevations of the study area range from 900 to 1200 m, and the average yearly precipitation ranges between 1200 to 1400 mm. According to the geological map (scale 1:50,000) of Dubertret [28], the study area is mostly made of limestone and dolomite, with sandstone in some places. This information is important for understanding the RF potential of the area.

In order to simulate the trajectories of RFs using 3D modeling software, for both earthquake-triggered and rainfall-triggered cases, it is essential to use the RF map.

Four aerial photos, with a scale of 1:10,000 and taken in 1999, were used in this study to generate the RF map using digital elevation model (with 5 m spatial resolution and 10 m vertical resolution) ArcGIS software (v10.8). The focal length of the used camera (of 153.401 mm) was used to correct the distortions in these photos using ERDAS Imagine software (v16.7.0). The corrected photos, of the study area covering 12.7 square kilometers, were then combined and digitized in ERDAS imagine. The photos cover large spatial extents of five villages: Mtein, Mchikha, Bzebdine, Hasabaya El Matn, and Qaaqour. The lengths of most slopes that have experienced previous RF hazards range between 250 and 350 m.

For instance, 209 deposited RFs (from previous RF events), of different rock sizes, were detected and mapped using the findings and information acquired throughout field trips (in June 2021) and aerial imaging [29]. The rock, already deposited from previous rockfall events, of sizes varying from 0.7 to 2284 cubic meters, were dispersed throughout the study area. The majority of the RFs had a volume of 25 cubic meters. Observations also revealed that certain built structures were in the trajectory of the potential falling rocks. These structures were built between the location of the RF source and their accumulation places.

The map (Figure 2) shows the locations of previous RFs, their rock volumes, and their maximal reaching distance (observed travel distance) as well as their release areas.

The maximum reaching points of RFs, represented by the lines in red in the map, varied in different areas because of the different characteristics of the landscape, such as the slope steepness, the rock type on the slope, and any obstacles that rocks may have encountered. In addition, all buildings that were observed in the orthophotos (existing before the year 2000) were digitized using the GEOAI Urban tool [30] and presented in the map. New buildings that were built during and after the year 2000 were digitized manually using the satellite basemap in ArcMap software (Version 10.8).

2.2. Runout Simulations

RF runout simulations are based on the shape; material properties; volume and density of the rocks; the slope angle and length (slope profile); the surface roughness of the slope, represented by the tangential coefficient of restitution; and the soil type, represented by normal coefficient of restitution [31,32]. The RF map, digital elevation model (with 5 meters spatial resolution and 10 meters vertical resolution), rock properties (rock density, block dimensions and shape), and the soil types and slope surface roughness, the first of which is directly related to normal (nCOR) and the latter of which is directly linked to the tangential (tCOR) coefficients of restitution, are inserted as inputs to run the RF runout simulations in Rockyfor3D software (Version 5.2.15) [33], a 3D modeling software for probabilistic and process-based simulations of RFs.

The mass of the RF blocks is calculated using the density and the volume of the rock at the source. The aerodynamic drag and rolling ability of the blocks are calculated based on their size and shape. The friction between each block and the slope is calculated using the roughness of the slope surface determined by the soil properties. The tangential coefficient of restitution measures rebound in a direction parallel to the contact surface, whereas the normal coefficient of restitution measures rebound perpendicular to the contact surface. These coefficients describe the ratio of rebound to incidence velocities during the collision.

The surface roughness is also determined by the heights of obstacles encountered, in a scattered pattern, along the slope. The heights are measured from the downward perspective of the slope [33]. The roughness in Rockyfor3D is represented by three classes: rg70, rg20, and rg10, which are used to calculate the tangential coefficient of restitution. These classes represent the probability of encountering different roughness levels, with rg70 being the most common (70% probability) and rg10 being the least common (10% probability). Each class is assigned to a roughness value, which represents the height of the obstacle that a falling rock would encounter. For example, if rg70 is equivalent to 0.05 meters, then 70% of the obstacles on the slope have a height of 5 centimeters. Two sets of simulations were conducted in Rockyfor3D. The first set simulated earthquake-triggered RFs, while the second set simulated RFs that occurred in fully moisture-saturated rocks and wet soils (representing the rainfall triggered RFs case). Among several output files generated by Rockyfor3D, only two outputs will be used in the analysis: (a) mean kinetic energy and (b) probability of reach. Using the said outputs of these simulations, the intersection impact points of the RF runout trajectories and the existing buildings prone to RF risks could be determined.

It is worth noting here that the directions by which rocks would fall from the source area to the accumulation zone were taken into consideration.

The following subsection describes the initial conditions (earthquake) that were used to simulate RFs in the software. The simulation parameters are derived from the collected data, such as the size and distribution of boulders observed in the aerial photos and during field visits, and these parameters are then calibrated to ensure that the simulated RFs match the real-world data as closely as possible.

2.2.1. Case of Earthquake

Three major faults in Lebanon (Yammouneh, Serghaya, and the Mount Lebanon thrust) are the most likely to generate earthquakes with magnitudes greater than 7. Additionally, secondary faults may cause earthquakes with magnitudes up to 6.5. Accordingly, Lebanon is at risk of experiencing moderate to severe earthquakes [34,35]. Harajli et al. 2002 [35] created a peak ground acceleration (PGA) map for Lebanon based on a 0.1 probability of exceedance in 50 years with an ultimate earthquake magnitude of 7. This map was used to calculate the initial velocity of RFs at the time of detachment, which is a parameter that is added to the simulation in the case of an earthquake. According to Harajli et al. (2002) [35] the local PGA is 0.20 m/s² (g).

Mavrouli et al. (2009) [36] have indicated that the PGA along the slope is obtained by linear interpolation between the acceleration at the base (PGAb) and that at the slope crest (PGAcr). The authors have indicated that the acceleration at the base is 0.20 g and the acceleration at the crust is 1.5 × 0.20 = 0.3 g. Therefore, the value of the PGA for the slope face (PGAsf) is 0.27 g, which was interpolated using the values of the PGA at the base and crest and based on the average altitude of the study area (1100 m).

In other words, the PGA map was used to determine the initial velocity of RFs at the time of detachment in the event of an earthquake. The PGA along the slope (PGAsf) is 0.27 g, which was calculated using linear interpolation.

According to Saroglou et al. (2018) [37], the initial horizontal velocity of the block at the time of detachment could be determined by taking into consideration the equilibrium of the produced work and the kinetic energy, as indicated in the following equation:

$${\mathrm{V}}_{\mathrm{x}}=\sqrt{2\times \mathrm{P}\mathrm{G}\mathrm{A}\mathrm{s}\mathrm{f}\times \mathrm{s}}$$

(2)

To start moving down the slope, the rock must be displaced by a small distance, denoted by Saroglou et al. (2018) [37], who set the value of s to 0.05 m in the calculations. Accordingly, the block’s initial horizontal velocity is 0.519 m/s, and its initial vertical velocity is considered to be negligible [37]. A value of 0.5 m was added to the initial height of the falling rock in the simulations, to account for the effect of the earthquake as Rockyfor3D software does not allow users to input the initial velocity or energy of a falling rock as a parameter. This value (0.5 m) is equivalent to the initial horizontal velocity of 0.519 m/s.

2.2.2. Case of Rainfall

The study area is characterized by a heavy rainfall pattern [38]. Most of the rocks in the area are karstified limestone, dolomitic limestone, and sandstone [39]. These rocks are porous and fractured, which allows water to easily seep into the pores and cracks. This increases the rock weight, the water pressure in the pores, and the width of the cracks during rainfall events. The wider the joints and fractures in a rock, the weaker it is [40]. In addition, water erodes the soil around rocks [41] and reduces the shear strength of soil materials [42]. Thus, rainfall is assumed to significantly increase the mass of a rock block without changing its volume, leading to an increase in the rock’s density. Heavy rainfall affects not only the rocks on a slope but also the materials that make up the slope itself. The high moisture content caused by heavy rainfall slows down rocks as they move down the slope [42]. In the simulations, the density of the rocks and soil was increased to account for the weight of the water they absorbed. The normal restitution coefficient, which represents rebound perpendicular to the contact surface when a rock bounces off the ground, was reduced to account for the dam** effect of wet soil during impact [42].

The rainfall-triggered RF simulation was used to study how the moisture content of the soil and rocks affects the rebound velocity of rocks. The normal coefficient of restitution (nCOR) was lowered because the moisture content of the soil is inversely proportional to the coefficient of restitution. This means that the wetter the soil, the slower rocks will bounce off of it [42]. The coefficients of restitution are indicated in the following equations [43]:

$$\mathrm{nCOR}=\frac{{\mathrm{v}}_{\mathrm{r},\mathrm{n}}}{{\mathrm{v}}_{\mathrm{i},\mathrm{n}}}$$

(3)

$$\mathrm{tCOR}=\frac{{\mathrm{v}}_{\mathrm{r},\mathrm{t}}}{{\mathrm{v}}_{\mathrm{i},\mathrm{t}}}$$

(4)

where v_i,n is the normal incident velocity, v_i,t is the tangential incident velocity, v_r,n is the normal rebound velocity and v_r,t is the tangential rebound velocity.

Hence, in the simulation software, the original soil type value [39] was lowered to reduce the normal restitution coefficient, since Rockyfor3D software has seven standard soil types, each with a different normal restitution coefficient. The average density of the rocks was chosen based on their nature, with limestone having a density of 2550 Kg/m³ and sandstone having a density of 2350 Kg/m³. Nonetheless, the rocks in the simulation are assumed to be fully moisture saturated. To match the densities of fully moisture-saturated rocks [44], the densities of sandstone and limestone were increased from 2350 kg/m³ and 2550 kg/m³ to 2450 kg/m³ and 2650 kg/m³, respectively. The modeling parameters (rock density, boulder volume, boulder shape, surface roughness, digital elevation model) needed to run the simulations of 3D rockfall (RF) runout models in Rockyfor3D are included in each raster map in ESRI ASCII grid format. When creating the input data, the worst-case scenario was taken into account. This scenario considers the spherical shapes of rocks as well as the absence of forests/trees that could block the rock movement along the slope during RF.

2.2.3. Collection and Calculation Risk Equation Parameters

The annual cost of damage caused by RFs is referred to as RF risk. It is calculated by factoring in the likelihood of an RF (occurrence probability), the vulnerability of a structure to RF damage, and the economic value of the structure. The geology of the area, the slope angle, and the presence of fissures or other discontinuities in the rock all influence the possibility of an RF.

The vulnerability of a structure to RF damage is determined by its strength (to resist the RF impacting the structure), construction, and the materials employed in its construction.

The RF economic risk value, in the present study, is represented by the cost of repairing or replacing the damaged structure [10], without taking into account any human or any additional loss. Thus, the risk value of an RF, in a particular area, could be calculated by multiplying the said three factors using the following equation [22,24]:

$$\mathrm{R}={\displaystyle {\sum }_{\mathrm{M}\mathrm{i}}}\mathrm{P}\left(\mathrm{M}\mathrm{i}\right)·\mathrm{P}\left(\mathrm{X}\mathrm{J}|\mathrm{M}\mathrm{i}\right)·\mathrm{P}\left(\mathrm{T}|\mathrm{X}\mathrm{j}\right)·\mathrm{V}\mathrm{i}\mathrm{j}·\mathrm{C}$$

(5)

where R is the risk value resulting from rockfalls of magnitude Mi and intensity j on a vulnerable element positioned at a distance X from the RF source (the magnitude Mi is typically quantified by the volume of rock that falls), P(Mi) is the probability of rockfall occurrence of magnitude Mi, P(Xj|Mi) is the reach probability of a rockfall with intensity j reaching a spot at a distance X from the RF source (the intensity j is typically quantified by the kinetic energy created by the RF at a distance X from the RF source), P(T|Xj) is the probability of the element being at spatial point X at moment T of the rockfall, Vij is the element’s vulnerability to a rockfall of magnitude Mi and intensity j, and C is the value of the element at risk.

2.2.4. Probability of Occurrence—P(Mi)

According to Jaiswal and van Westen (2009) [45], the temporal probability is obtained by multiplying the annual exceedance probability (AEP) by the probability of landslide initiation, in the case of threshold exceedance. In this study, it is assumed that, whenever the threshold is exceeded, the RF is initiated, so the probability of landslide initiation is equal to one. It is worth noting that the threshold capable of inducing a rockfall is considered in this study to be equal to 5 according to the Richter scale in the case of earthquake, and to 100 mm cumulative precipitation in the case of rainfall over a period of 3 days.

The annual exceedance probability (AEP) is the chance that a general event of a certain magnitude/intensity/both will happen in any given year [24]. To determine the probability of RF, triggered by a factor whose intensity exceeds a certain threshold, within a certain time period, a Poisson probability model can be used [45,46,47]. The Poisson distribution model that can be used to calculate the probability of having a particular number of RFs in a time period t (N(t)) is given in the following:

$$\mathrm{P}\left[\mathrm{N}\right(\mathrm{t}) \ge 1]=1-{\mathrm{e}}^{-\mathrm{t}/\mathrm{\mu }}$$

(6)

The mean recurrence interval between successive RFs (μ), can be obtained from the multi-temporal RF inventory data.

Case of earthquake:

Previous research has shown that only earthquakes with a magnitude of 5 or greater can trigger slope failures [48,49], as strong ground motions have been known to trigger RFs in the past [50].

As large magnitude earthquakes (magnitude above 5 according to the Richter scale) have occurred in previous centuries in Lebanon, Brax et al. 2019 [51] assessed the published studies on earthquakes in Lebanon and compiled a list of events with estimated micro-seismic epicenters. This list includes the earthquakes that occurred in the years 551 (with magnitude 7.3), 1170 (magnitude 7.3), 1202 (magnitude 7.6), the two earthquakes in the year 1759 (on 30 October with magnitude 6.6 and on 25 November with magnitude 7.4) and lastly the earthquake of the year 1837 (magnitude 7.0) [52].

The return period of large earthquakes hazards in Lebanon is estimated to be 475 years [52,53]. In the present study, the return period is considered the same (475 years) according to aforementioned references. According to Equation (6), and during the time period (t) of 1 year, the temporal probability (annual probability) is equal to 0.002.

Case of heavy rainfall:

Intense rainfall can alter hydrogeological conditions and water levels, triggering surface movements of soils and altering materials. These movements can initiate RFs, landslides or debris flows. Knowing that there is a direct relationship between the amount of accumulated rainfall and the occurrence of instabilities, one can say that, generally, the greatest probability of these instability events occurs during the periods of heavy precipitations [54].

A recent study conducted by Leyva et al. (2022) [54] has investigated the relationship between rainfall and RFs over the period 2010–2016 in the Anaga Nature Reserve, a mountainous region with steep slopes in Tenerife. The cumulative rainfall over a period of 3 days was the amount of rainfall considered in this study.

The results show that (a) there is a strong correlation between rainfall and RFs, as the probability of an RF occurrence increases significantly after periods of heavy rainfall, and that (b) the amount of rainfall required to trigger an RF varies depending on the slope angle. In the latter context, the authors found that, on slopes with an angle of 45 degrees, the amount of rainfall required to trigger an RF is about 100 mm over a period of 3 days, as well as 200 mm accumulated rainfall for slope angle of 60 degrees.

With respect to the amount rainfall that could initiate slope instability, the present study uses the same method adopted by the aforementioned research [54], which considered the cumulative rainfall amount over a 3-day period.

Further, and to be more conservative while forecasting the slope instability risks, the worst-case scenario, which considers the least amount of rainfall (100 mm) corresponding to the shallowest slope angle (45 degrees), was taken into account in the present study.

In the present article, daily precipitation data recorded from the year 1981 through to the present (2023) was obtained from the Earthdata webpage of NASA [55]. By analyzing the obtained data, the cumulative heavy rainfall (≥100 mm) over three consecutive days was detected in the years 2010 and 2019, in the months of December and January respectively. Thus, in the present research, the return period for rainfall-triggered RFs is assumed to be 9 years, while, according to Equation (6), which considers t = 1, the annual probability of RF occurrence is calculated and is then equal to 0.105.

2.2.5. Reach Probability—P (Xj|Mi)

Reach probability is the percentage of blocks that reach a distance from the RF source cliff. A simple calculation can be conducted by dividing the number of blocks that passed through a given distance by the total number of falling blocks [56]. This is typically used to determine the RF runout zone, which represents the area where blocks are likely to stabilize after they have fallen from the cliff. For this probability, Rockyfor3D calculates the reach probability in each location of the zone of interest.

The reach probabilities are then calculated for each scenario, i.e., of earthquake and rainfall cases as presented in Figure 3 and Figure 4, respectively. Finally, intersection points are created between the runout zones of simulated RFs and the existing buildings already rasterized. Accordingly, the reach probabilities are attributed to these points.

2.2.6. Probability of the Element Being at Point X at the Time of the Rockfall Occurrence P (T|Xj)

The spatial probability of the element being at point X at time T of the rockfall occurrence is assumed to be equal to the point of intersection created between simulated RFs and buildings. This is due to the fact that buildings (containing the intersection points) are fixed in position and are located on the trajectory or fallen rocks.

2.2.7. Physical Vulnerability—(Vij)

In the context of RF hazards, many authors have defined two types of corresponding vulnerabilities: societal vulnerability [57] and physical vulnerability [58]. The first of these is related generally to human injuries and deaths that are the result of RF events.

RF physical vulnerability is the degree to which a set of assets can be damaged or lost due to RFs [59]. Both types of vulnerability represent a key factor in RF risk assessment and management, as their estimation is essential for making reliable predictions and mitigation plans for RF impacts. Nonetheless, the present study is limited to an examination of RF risk from a physical vulnerability perspective.

According to UNDP/UNESCO 1982 [60], four main groups of methods exist for assessing the level of vulnerability of a structure: (i) categorization methods, (ii) inspection and empirical (rating) methods, (iii) analytical methods, and (iv) experimental methods.

The categorization methods classify structures into typological classes based on their structural type, construction materials, and construction techniques. They are employed because different structural typologies have different vulnerabilities to natural hazards.

The inspection and empirical (rating) methods assign numerical values to each structure based on its condition. Generally, this is done by inspecting the structure and looking for signs of damage, such as cracks, deterioration, and poor construction practices.

The analytical methods calculate the resistance of a structure against natural hazards damage by using civil engineering equations and computer models.

The experimental methods involve testing the structural components and entire structures to determine their properties. This type of vulnerability assessment for structures is considered the most accurate; however, it is also the most expensive and time-consuming method.

In the present study, a combination of the first three of these four methods is applied. Accordingly, the vulnerability was assessed using the combined methodology developed by Li et al. (2010) [61], and was calculated according to the following:

$$\mathrm{V}=f\left(\mathrm{I},\mathrm{R}\right)=\left\{\begin{array}{cc}2{\left(\frac{\mathrm{I}}{\mathrm{R}}\right)}^2& \mathrm{if} \frac{\mathrm{I}}{\mathrm{R}}\le 0.5\\ 1-2{(1-\frac{\mathrm{I}}{\mathrm{R}})}^2& \mathrm{if} 0.5<\frac{\mathrm{I}}{\mathrm{R}}\le 1.0\\ 1.0& \mathrm{if} \frac{\mathrm{I}}{\mathrm{R}}>1.0\end{array}\right.$$

where V represents the vulnerability, I represents the RF intensity, and R represents the resistance of buildings.

The RF intensity I in this study is scaled from 0 to 1, based on the previous research by Lateltin et al. (2005) [62]. The authors indicated that I is determined mainly by the RF kinetic energy and that, therefore, the intensity of RFs can be estimated according to

Table 1. The estimation of RF intensity (Source: Singh et al. (2021) [8]).

Rockfall Intensity	Kinetic Energy: E (KJ)	Proposed I Value
Low	Less than 30	0.2
Moderate	30–300	0.6
High	More than 300	1

.

The resistance of buildings, R, is calculated as the following equation using resistance factors associated with the following four variables: structure typology (ε1), state of maintenance (ε2), construction quality (ε3) and number of floors (ε4) as determined by Li et al. (2010) [61].

$$\mathrm{R}=({\mathrm{\epsilon }1\times \mathrm{\epsilon }2\times \mathrm{\epsilon }3\times \mathrm{\epsilon }4)}^{1/4}$$

(7)

With respect the first variable, buildings must be categorized into different structural typologies where each topology has a specific resistance factor ($\mathrm{\epsilon }$1). The resistance factors for the three structural typologies are defined as follows: brick masonry ($\mathrm{\epsilon }1$ = 0.3), stone masonry ($\mathrm{\epsilon }1$ = 0.6), and reinforced concrete ($\mathrm{\epsilon }1$ = 1).

Once the buildings have been categorized, they must be inspected and assigned numerical values according to rating (empirical) methods in order to represent the state of maintenance ($\mathrm{\epsilon }$2). This involves rating the structure by looking for signs of damage and poor construction practices, as follows: damaged ($\mathrm{\epsilon }2$ = 0.2), distressed ($\mathrm{\epsilon }2$ = 0.4), good ($\mathrm{\epsilon }2$ = 0.8), and excellent ($\mathrm{\epsilon }2$ = 1).

Furthermore, the construction quality evaluation is determined by several factors, such as the type of construction, the materials used, and the age of the structure. The higher the quality, the higher the resistance and the lower the vulnerability rating. The resistance values corresponding to the quality are defined as follows: low ($\mathrm{\epsilon }3$ = 0.2), moderate ($\mathrm{\epsilon }3$ = 0.5), and high ($\mathrm{\epsilon }3$ = 1).

The fourth resistance factor corresponds with the number of floors, where the higher the building, the lower the vulnerability rating. This assumption, made by Singh et al. (2021) [8], is based on the fact that taller buildings typically have more depth and solid foundations as well as core walls and shear walls. These resistance factors are defined as follows: 1 floor ($\mathrm{\epsilon }4$ = 0.2), 2 floors ($\mathrm{\epsilon }4$ = 0.5), 3 and more floors ($\mathrm{\epsilon }4$ = 1).

Consequently, the most critical parameters to be determined are I and R. The first parameter, I, was obtained through the simulation results for each scenario. By examining the simulation output, specifically the mean kinetic energy map, the kinetic energy value was extracted from the raster map at the point where the rockfall impacts the side of the building. This point represents the potential impact site on the building structure.

The coefficient R was calculated using Equation (7), wherein distinct values were assigned to each parameter (ε1, ε2, ε3, and ε4) in the equation for every building. These assignments were based on the assessment of the buildings’ characteristics and conditions situated along the trajectory of rockfall propagation.

2.2.8. Value of the Element at Risk

The cost of repairing building damage caused by RFs is difficult to estimate because it depends on many factors, such as the severity of the damage, the type of building, and the cost of materials and labor. Additionally, the cost of repairing damage can vary significantly depending on the location of the building. For example, it may be more expensive to repair building damage in a peripheral area than in a more centrally urbanized area. The area of a building that is exposed to RFs could be employed to accurately measure the value of the element at risk [63]. Considering the area of buildings that are exposed to RF risk leads to a better understanding of the potential damage caused by the RFs. Based on expert opinion and according to experience in assessing building damages [64], the cost of damage is estimated, in this study, to be 400 USD per square meter. This cost is expensive considering that the zone of interest is a peripheral area and requires additional cost to transport materials.

3. Results and Discussion

Recalling the previous section, 1000 simulations for each scenario (earthquake and rainfall scenarios), were performed. An additional value of the horizontal kinetic energy was added to represent the effect of earthquakes on rock blocks in the release area. This is expressed by lifting up each block by 0.5 m, as suggested in Section 2.2.1. However, in the case of the rainfall scenario, two simulation parameters were changed differently. The change of the first is represented by lowering the normal coefficient of restitution (soil type), which leads to a dam** effect on the rock block while moving down along the slope. The second parameter is the density of rock block itself, and this is increased due to water absorption during or right after a rainfall period. The change in the second parameter leads to an increase of moisture content and therefore an increase of the density of rock block. The densities of sandstone and limestone were increased from 2350 kg/m³ and 2550 kg/m³ to 2450 kg/m³ and 2650 kg/m³ (see Section 2.2.2).

The results reveal that 48 buildings in the study area could potentially be hit by RFs in both the earthquake and rainfall scenarios. However, in the earthquake scenario, an additional 5 buildings could potentially be hit, for a total 53 damaged buildings. This suggests that earthquakes can dislodge more rocks and cause them to fall over a wider area.

Regarding the damage size, the study also shows that earthquakes will potentially present more damage to buildings in the RF impact area than in the case of rainfall-triggered RF. Findings show that the risk of damage to the majority of buildings was higher in the earthquake scenario compared with the rainfall scenario. This is likely because the earthquake vibrations can tremble buildings and make them more vulnerable to damage from falling rocks. This is also related to the observation that the highest impact energy resulting from the case of earthquake was 10,673.7 KJ, while, for that of rainfall, it was 10,451 KJ.

The mean kinetic energies for RFs in both scenarios are presented in Figure 5 and Figure 6. It is worth noting that the potential highest kinetic energies in both scenarios are observed in the same spots, and that this is due to the steepness of slope, which is is around 60 to 65 degrees. Additionally, and with respect to building conditions that affect their vulnerability, the field survey of buildings that could potentially be hit by RFs has revealed that 85% of the buildings are made of concrete, 11% made of stone masonry and 4% made of brick masonry. In addition, 53% of these buildings are well-maintained and in excellent condition, 32% are in good condition and 15% are in damaged/distressed condition. Notably, nearly one-third (28%) of the buildings are poorly constructed.

Regarding vulnerabilities, the calculation results are equal to 1 for all of the buildings except for one of them (0.68 for earthquake and 0.62 in rainfall case), this is explained by (1) the high intensity of RFs, as mean kinetic energy was higher than 300 KJ, and (2) the low resistance of buildings against RF impact.

These results are combined with the impact probability and the value of elements at risk according to equation 5, allowing one to compute the annual cost of RF risk on buildings (as shown in Figure 7 and Figure 8). Results show that the highest cost of RF risk is of 232 USD in case of earthquake, and 10,511 USD in the case of rainfall, per annum.

The proposed approach for risk assessment (quantitative method) in this study, is limited by the absence of statistically representative historical catalogues of RF events as well as their impacts and damages in Lebanon. This makes it difficult to accurately characterize the frequency distribution of RF volumes for individual sites, which is essential for using the approach in practice.

The observed size of deposited rockfalls varies significantly, ranging from 0.7 to 2284 cubic meters. However, for the purpose of this study, a volume of 25 cubic meters was selected as the initial estimate for the volume of the rockfalls that originated from the release areas. This choice was based on the high frequency of this size within the study area, without taking into consideration the fragmentation. Given the broad spectrum of sizes observed, it was assumed that this volume could be the better representative of typical rockfalls from release areas. Larger rockfalls are expected to yield increased runout distances and higher kinetic energies on the slope, while smaller rockfalls are anticipated to have shorter runout distances and lower kinetic energies.

Risk assessments are also affected by the quality of the database that is used to calculate the vulnerability functions. Additionally, the data needed to fully quantify the value of elements at risk (e.g., buildings, infrastructure) is not always available. This introduces further uncertainty into the estimated costs.

Despite these limitations with respect to RF events, dating the RFs can improve the estimated return periods, as it can provide information on how many boulders of a given size could potentially fall in a specific time period.

The locations where RFs commence (release areas) are difficult to detect and to correctly map. Nonetheless, with the help of aerial photographs and field investigation, the release zones on the unstable slope were located and mapped.

Additionally, partially detached rock blocks and rock masses overlying joints, which are orientated in a way that makes them more likely to break, are indicators of possible RF release zones. Finally, the simulations were carried out where boundaries of the release areas were chosen, in order to cover the locations of all unstable blocks on hills.

4. Concluding Remarks and Future Research

Although no historical information on RFs is available, the presented case study demonstrates that advanced numerical modelling techniques can be used to establish a valuable RF risk evaluation method that can be implemented to improve the quantitative risk assessment at various levels of complexity based on the available data. The simulation was validated by comparing the results of each simulation to the actual observed rockfall deposits. It was noted that the locations of the deposited rockfalls in the simulation aligned with the actual deposition sites.

This study developed a new methodology with which to calculate the risk of RFs for individual buildings in a rural area. The study first used a computer model to simulate, using a computational approach, the RF runout zone and other characteristics, such as the speed, energy and reach probability of the RFs. Then, the study considered both the intensity of the RFs and the resistance of the buildings to calculate the vulnerability of each building. The resistance of each building encompasses factors such as the construction structure, number of floors, building type, and maintenance state. Moreover, the economic impacts of RF events are estimated. Finally, the study created a risk map for the study area by combining the risk maps for different buildings.

The proposed risk analysis method for RFs has several advantages over existing methods. First, it is particularly well-suited for risk problems that are spread out over a large area. Second, it can combine all stages of risk assessment into one process, which can make the assessment more efficient and precise. Third, it decreases the number of assumptions that need to be made about the physical components of risk, such as the hazard, probability, and energy of impacts on structures.

However, there are still some concerns about the method because there is a lack of statistically reliable data on both RF events, damage, and cost of assets. This lack of data can make it difficult to validate the results of the risk assessment. Nevertheless, the proposed method provides a good estimate of risk, which is essential for evaluating the effectiveness of RF mitigation measures.

Generally, the proposed risk analysis method is a promising new approach for assessing the risk of RFs, as it is more comprehensive and accurate than existing methods. However, that method is still under development, requiring more data to validate the results. Nevertheless, and despite some limitations in the risk assessment, the risk map can still be used to identify areas at high risk of RFs and to develop mitigation measures. Mitigation measures can reduce the hazard or susceptibility of the community to RFs, and they are especially important in high-risk areas. Given the high number of inhabitants in this study area, upcoming studies are needed to carry out an assessment of the RF risk on human lives.

The application of quantitative risk assessment methods is not limited to rockfalls, it can also be applied to other types of landslides. Identifying hazards, assessing the probability of occurrence, and estimating the potential impact are necessary across different types of landslides. However, specific models and techniques may need to be adapted or developed based on the unique characteristics of each landslide type. It is essential to emphasize the modern significance of machine learning algorithms, which have enormous potential in advancing the study of rockfalls and can be fundamental in further research. For instance, the study by Rafiei Sardooi et al. (2021) [65] on landslide risk assessment provides a good example of this potential. Their methodology not only offers a framework for assessing the risk of various landslide types but also delivers remarkably valuable results. Machine learning (ML) and artificial intelligence (AI) can significantly contribute to quantifying risk in several ways, such as automated geological feature map** and landslide simulations, etc. AI can integrate all of the information and data into a single interface in order to analyze landslide risks and develop risk maps using various models and methods like the finite element method (FEM), limit equilibrium method (LEM), and quantitative risk assessment (QRA).

Moreover, the integration of both physical and social science aspects in comprehending and evaluating landslide mitigation measures by community perception represents a groundbreaking approach. This entire perspective not only raises innovation but also raises public trust in the efficacy of such measures [66].