From rockfall source area identification to susceptibility zonation: a proposed workflow tested on El Hierro (Canary Islands, Spain)

Rockfall modelling is a common topic in the landslide literature, but a comprehensive workflow for rockfall susceptibility zonation remains a challenge. Several aspects of the modelling, such as rockfall runout simulation, are consolidated, but others still show inconsistencies and ambiguities, such as the source area identification or the criteria to obtain probabilistic susceptibility zonation. This study proposes a workflow for rockfall susceptibility zonation at the regional scale that integrates (i) source area identification criteria, (ii) deterministic runout modelling, (iii) approaches for the runout classification, and (iv) robust procedures for validation and comparison.

The workflow is tested on El Hierro Island (Canary Islands, Spain) and considers the effect of different methods to identify the rockfall source areas that are used as input data for rockfall runout modelling. The runout outputs are classified to derive rockfall susceptibility zonation considering different types of classification (i.e. supervised versus unsupervised methods). The source area identification reflects situations with limited data availability and scenarios with a large amount of topographic, geological, and geomorphological information. The first approach is based on slope angle thresholding, the second uses a statistical method based on empirical cumulative distribution functions (ECDFs) of slope angle values, and the third involves the combination of multiple multivariate statistical classification models where the source area is the dependent variable and thematic information is the independent variables. The source area maps obtained from the three methods are utilized as inputs for a rockfall runout model (STONE) to derive rockfall trajectory count maps. Two classification approaches are applied to generate probabilistic susceptibility maps from the trajectory counts: unsupervised and supervised statistical methods using distribution functions. The unsupervised classification employs only the rockfall trajectory counts as input, whereas the supervised classification requires additional data on the areas already affected by rockfalls. To complement the workflow, statistical methods and metrics are proposed to verify, validate, and compare the susceptibility zonation.

1 Introduction

Rockfalls are dangerous natural hazards with a relevant socioeconomic impact worldwide (Borella et al., 2019; Mateos et al., 2020). Changes in environmental conditions related to the growth of the population, land use intensification, and industrial development have the potential to increase the impact of rockfalls in many different regions (Farvacque et al., 2019; Othman et al., 2021; Santangelo et al., 2020). In addition, climate change is expected to modify precipitation patterns, which will have the effect of increasing the frequency and extension of rockfalls (Gariano et al., 2015; Sarro et al., 2021). As a consequence, there is an increasing interest in improving the reliability and accuracy of tools and products able to support rockfall management and mitigate its impact (Noël et al., 2021; Omran et al., 2021; Santos et al., 2024).

A rockfall susceptibility estimate shows where rockfalls are likely to occur. A susceptibility zonation is a map that classifies the territory based on the likelihood of rockfall occurrence. Several techniques have been developed for rockfall susceptibility zonation, which can be broadly classified into qualitative and quantitative approaches. Qualitative approaches may include geomorphological analysis and heuristic techniques, whereas quantitative approaches are different statistical or deterministic-based analyses (Reichenbach et al., 2018).

Rockfall runout modelling allows us to obtain information on the spatial distribution of boulder trajectories and their velocity, energy, and heights (Carlà et al., 2019; Gallo et al., 2021) and plays a relevant role in rockfall assessment, supporting the identification of rockfall-prone areas and the characterization of boulder behaviour (Crosta et al., 2015; Pfeiffer, 2019). A wide range of algorithms or models is currently available for the calculation of runout zones. Nevertheless, deterministic modelling remains inherently uncertain due to insufficient data on the mechanical and geometrical properties of the terrain, limiting the reliability of trajectory identification. To address this limitation, stochastic approaches have been introduced to account for the variability in the input parameters. In the literature, different modelling approaches have been proposed based on data availability, the environmental setting, and the type of analysis. An incomplete list consists of STONE (Guzzetti et al., 2003; Sarro et al., 2020), RocPro3d (Sarro et al., 2014, 2018), Hy-Stone (Dinçer et al., 2016; Lanfranconi et al., 2020), RAMMS (Dhiman and Thakur, 2021), RockyFor3D (Francioni et al., 2020; Robiati et al., 2019), and RocFall (Kakavas et al., 2023; Pérez-Rey et al., 2019).

The output of runout models is commonly used to estimate the rockfall susceptibility degree by classifying rockfall trajectory counts (Dorren et al., 2023; Nanehkaran et al., 2022; Noël et al., 2023). This susceptibility measures the degree to which terrain can be affected by future slope movements. In other words, it is an estimate of “where” landslides are likely to occur and, in mathematical language, can be defined as the probability of the spatial occurrence of slope failures, given a set of geo-environmental conditions (Reichenbach et al., 2018).

Rockfall modelling, both stochastic and deterministic, presents errors associated with the input data employed to replicate the rockfall process (Straub and Schubert, 2008). The inaccuracy in defining rockfall source areas is highly relevant in modelling, since source areas provide the starting state for rockfall trajectories (Frattini et al., 2013; Rossi et al., 2020).

The location of source areas depends on several characteristics, such as slope morphology, lithology, and geological discontinuities (Alvioli et al., 2021; Sarro et al., 2018; Yan et al., 2023). At the local scale, in situ analyses commonly involve discontinuity characterization and escarpment recognition. Frequently, logistical and safety issues in the field constrain these methods. Remote sensing techniques, such as laser scanners and uncrewed-aerial-vehicle-based (UAV-based) photogrammetry, are nowadays widely used to address these limitations and obtain detailed observations of slopes (Gallo et al., 2021; Giordan et al., 2020; Sarro et al., 2018). Although both fieldwork and remote sensing methods are successful at the local scale, their utility at the regional scale is limited.

Many methods with different degrees of complexity have been proposed to identify rockfall source areas at the regional scale based on deterministic, probabilistic, or statistical approaches (Muzzillo et al., 2018). Deterministic methods identify rockfall source or detachment locations using models based on mechanical principles, while statistical methods are based on the analyses of historical catalogues of past rockfall events. For the probabilistic identification of source areas, supervised multivariate classification or machine learning models are employed to predict rockfall detachment locations (i.e. the dependent or grouping variable) based on a set of explanatory variables (i.e. independent variables).

Most of the approaches are based on the numerical analysis of digital elevation models (DEMs) and additional environmental datasets. Source areas can be identified by analysing local topography using surface slope thresholds, which denote the area with favourable conditions for boulder detachment. Larcher et al. (2012) proposed an equation to define rockfall source areas by linking the slope angle threshold and the resolution of a DEM. Rockfall source areas can also be identified empirically or derived from the decomposition of slope frequency distributions using morphometric methods based on the slope angle thresholds. Several studies determined a correlation between this threshold and the angle of internal friction of the rock massif (Loye et al., 2008; Paredes et al., 2015). Thus, the evaluation of slope frequency distributions can determine the angle of internal friction associated with each lithological unit of the rock massif, and it is used as the threshold beyond which the massif becomes unstable. In the same way, Loye et al. (2009) developed a model based on the Gaussian distribution of the slope angle values. According to the results of this slope angle distribution, for each morphological unit, the steepest slopes are selected as potential source areas (Zhan et al., 2022). Additionally, Wang et al. (2021) identify rockfall source areas controlled by rock mass strength using relief–slope angle relationships.

Other identification techniques at the regional scale are based on the analysis of remote sensing multi-temporal imagery, such as interpretation of orthophotos from optical aerial or satellite data. The use of distinctive imaging features or signs, such as scars or deposits, has been shown to be feasible in several different research projects (Liu et al., 2020; Mateos et al., 2016; Scavia et al., 2020). However, this technique is limited by the availability of satellite data and the difficulty of analysing some areas (shadowed slopes, steep slopes, and/or vegetation). Moreover, photo interpretation is time-consuming, and this often hampers its application over large areas (Alvioli et al., 2021).

Recently, advanced heuristic methods and statistical tools were proposed to identify the location of source areas with good results. A heuristic method depends on the site characteristics, and its application requires validation and special adaptation processes (Fernandez-Hernández et al., 2012). Conversely, statistical methods can be performed to assess different levels of likelihood based on geomorphological, geological, and geo-environmental factors. These methods, such as multivariate analysis, logistic regressions, or frequency ratios, are more flexible than heuristic methods but require training with representative data samples. Hybrid methods combine statistical and experimental methods, such as neural networks or machine learning decision analysis, to reduce the amount of data required and improve the accuracy of the results (Fanos and Pradhan, 2019; Rossi et al., 2020).

Rockfall modelling has seen considerable improvements in recent years; however, a significant challenge remains concerning the seamless integration of different modelling components spanning the source area definition to the susceptibility zonation. In addition, numerous studies have explored methods to identify source areas, but there are no specific studies that analyse how they influence the rockfall modelling results.

This study proposes a workflow that combines existing knowledge and new concepts and techniques to estimate rockfall susceptibility and to enhance its robustness and applicability, moving beyond case-specific solutions to provide a robust and transferable framework. The workflow integrates (i) three source area identification criteria, (ii) deterministic runout modelling, (iii) supervised and unsupervised approaches for the runout classification, and (iv) statistical methods and metrics to verify, validate, and compare the susceptibility zonation. The application at a regional scale on El Hierro Island (Canary Islands, Spain) allowed a preliminary evaluation of the workflow.

The article is organized as follows: Sect. 2 describes the test area, Sect. 3 presents the variety of methodologies employed, Sect. 4 presents the results, and Sect. 5 discusses the results and highlights the main conclusions.

2 Test site and data

2.1 Geographical and geological setting

The Canary Islands are a volcanic archipelago located in the Atlantic Ocean, within the African plate. The archipelago is made up of seven major islands (Fig. 1) and some smaller ones, which, together with underwater reliefs, form an extensive volcanic domain. The islands are the result of a long magmatic history that started 70 million years ago and has continued to the present with the recent volcanic eruption on La Palma (September 2021).

Figure 1 Areas used to classify and validate the simulated rockfall runout.

El Hierro is the westernmost and youngest island, with an area of 268.71 km² and a population of 11 147 inhabitants (Instituto Canario de Estadística (ISTAC), 2021). The climate is subtropical oceanic along the coast and is very mild and sunny for most of the year, with rainfall concentrated in October to March. Heavy storms are frequent, associated with intense rainfall and strong winds that often trigger landslides. The average temperature ranges between 19 and 25 °C, with maximum values in August.

The morphology of the island is the result of numerous volcanic events associated with important geological features. One of the most characteristic features of El Hierro is the occurrence of large landslides, which correspond to the escarpments of El Golfo, El Julan, and Las Playas, located in the north, southwest, and southeast, respectively (Fig. 1). The northern part, El Golfo, with cliffs that reach an elevation of more than 1100 m, is a hazardous area for rockfalls. During the period of 2011–2012, a submarine eruption took place about 2.5 km off the coast near the village of La Restinga. The highest seismicity was in the El Golfo area, with two earthquakes of magnitude 4.4 and 4.6 in mid-November 2011. The seismic events triggered rockfalls near the Los Roquillos tunnel, one piece of strategic infrastructure that connects the municipalities of Frontera and Valverde, the most populous villages on the island. After the event, the first field observations carried out by technicians of the Geological and Mining Institute of Spain (IGME-CSIC) allowed us to evaluate the cliff stability along the road HI-5, where the Roquillos tunnel is located. The report that was prepared showed a complex scenario for the analysis of rockfall hazard and the definition of source areas. The field surveys revealed that the dikes that outcrop on the escarpments of El Golfo and Las Playas are preferential rockfall source areas. Recently, on 14 March 2021, a large rockfall along the El Golfo escarpment alerted the population and caused a general alarm.

2.2 Available data and products

For El Hierro Island, the following data are available: (1) a digital elevation model (DEM) at a 5 m × 5 m resolution (Centro de Descargas del CNIG (IGN), 2024) that was used to compute morphometric parameters (elevation, slope, curvature, landform classification, etc.) and (2) lithological information derived from the geological map provided by IGME-CSIC at a scale of 1:25 000. The map was reclassified into seven geotechnical classes (Sarro et al., 2020; Rossi et al., 2020), ranging from class 1, which includes soft soils (such as lapilli and sand), to class 7, which includes extremely hard rocks (dikes and volcanic breccias).

In this paper, we have used different thematic data to identify source areas and to perform rockfall modelling and susceptibility zonation. The methods to identify source areas require diverse types of information. (i) Unsupervised slope thresholding (ST_RSA) and slope angle empirical cumulative distribution functions (ECDFs) (CDF_RSA) require only slope data; (ii) supervised ST_RSA and CDF_RSA require slope data and the location of source areas (i.e. normally mapped in the field; Rossi et al., 2020); and (iii) probabilistic identification (PROB_RSA), together with the location of source areas, exploits the following additional geo-environmental information as conditioning factors: topography parameters (i.e. slope, curvature, and aspect derived from the DEM), lithology, and the presence of dikes (Rossi et al., 2020).

For the runout modelling, the following additional data were used: (i) a sample of mapped rockfall deposits in polygon format for the supervised cumulative distribution function (CDF) analyses of rockfall trajectories, (ii) a sample of areas affected by or with no evidence of rockfalls for the receiver-operating-characteristic-based (ROC-based) model performance evaluation, and (iii) a sample of the rockfall boulder locations (i.e. silent witnesses) for the violin plot and boxplot susceptibility zonation.

The rockfall information used in the runout simulation classification, and validation was derived using diversified techniques and sources of information. Using field investigations conducted from 2012 to 2018 (47 records), aerial image interpretation (84 records), and data from the MOVES database (BDMoves, 2024) (78 records), we have identified rockfall deposits (red polygons in Fig. 1) that include single detached boulders (i.e. mapped as points; black dots in Fig. 1c) and talus deposits (i.e. mapped as polygons; blue polygons in Fig. 1d). Additionally, areas with no evidence of rockfall activity were recognized in the field by experts, with the support of geomorphological and topographical maps (i.e. green polygons in Fig. 1).

3 Methodology

The methodology proposed in this study can be formalized in a workflow that considers different steps (Fig. 2).

• 1.

  The first step is typical of any rockfall study, where available relevant data are collected (field surveys, photo interpretation, etc.).

• 2.

  The second step focuses on the identification of rockfall source areas, a critical input for the subsequent analyses performed using different approaches.

• 3.

  The third step is the deterministic rockfall runout modelling using the source areas from step (2) as input. The main output is a map of the cumulative count of rockfall trajectories.

• 4.

  The fourth step derives probabilistic susceptibility zonation through the classification of the trajectory count values from step (3). Unsupervised and supervised classification approaches based on the empirical cumulative distribution function (ECDF) are applied for this purpose.

• 5.

  The fifth step validates and verifies the susceptibility maps and assesses the reliability of susceptibility zonation using quantitative multi-criterion evaluation techniques and statistical metrics.

The five methodological steps, their application to the study area, and the results are illustrated in the following sections.

Figure 2 The rockfall modelling workflow.

3.1 Identification of rockfall source areas

The maps of the source areas are crucial data for the rockfall analysis. In our study area, we have used three different approaches: (i) a morphometric schema based on the slope thresholding, (ii) cumulative distribution functions (CDFs) that consider slope information and geology, and (iii) a probabilistic model.

3.1.1 Slope thresholding

This method (hereafter referred to as ST_RSA) relies on a simple morphometric approach, which identifies potential rockfall detachment zones as those areas with a slope angle above a given threshold. Even though rockfalls initialize mainly on steep slopes and the steepness of the hillslope surface can be used to identify potential source areas, it is more realistic to determine a slope threshold using distinctive evidence (e.g. deposits, inventory) rather than arbitrarily establishing one (Michoud et al., 2012).

According to Fu et al. (2021), more than 80 % of 2238 rockfall records collected in Sichuan (China) over the past 30 years occurred on hillslopes with slopes ranging between 30 and 50°, with most of them around 40°. As the result of a historical rockfall study in the Yosemite Valley (California, USA), Guzzetti et al. (2003) identified slopes above 60° as potential release points. In the region of Vaud county (Switzerland), Jaboyedoff and Labiouse (2011) determined slope thresholds between 47 and 54°. Based on the experience of the Trentino Geological Survey, Frattini et al. (2008) selected as source area cells with slope angles over 37° in Val di Fassa (Dolomites, eastern Italian Alps). Overall, most of the previous studies cited reveal slope thresholds over 30°.

Sarro et al. (2020) proposed a slope threshold over 40° on Gran Canaria (Canary Islands), an island with similar topographical and geological conditions to El Hierro. Detailed evaluations revealed that the source areas in Gran Canaria are primarily associated with hard, very hard, and extremely hard rocks, corresponding to geological types such as dikes and breccia, phonolite, massive basalt, trachyte, and ignimbrite. Considering that the geological context of El Hierro where rockfalls are observed is similar to that of Gran Canaria, we have defined the threshold above 40°. The map obtained using the slope thresholding method is binary, with 0 corresponding to stable areas and 1 to rockfall-prone detachment zones.

3.1.2 Statistical identification of rockfall source areas using slope angle ECDFs

For the second identification of rockfall source areas, we utilized the empirical cumulative distribution functions (ECDFs) of slope angle values (hereafter referred to as CDF_RSA).

An ECDF function returns the probability that a random variable is less than or equal to a given value (Lee et al., 2022). In mathematical terms, this is expressed by Eq. (1):

$$\begin{array}{}\text{(1)}& F_x\left(x\right)=P\left(X\le x\right)={\sum }_{t\le x}f\left(t\right),\end{array}$$

where F_X(x) denotes the ECDF of a random variable X, whose probability distribution is f(x).

ECDF has lower and upper limits of 0 and 1, respectively, and gives a cumulated probability that increases with the x value. Equation (2) shows the values taken by ECDF or F_X(x) for infinite boundaries of the random variable and Eq. (3) the relation between F_X(x) values for successive values of x.

$$\begin{array}{}\text{(2)}& F_x\left(-\mathrm{\infty }\right)=\mathrm{0},F_x\left(\mathrm{\infty }\right)=\mathrm{1},\end{array}$$

$$\begin{array}{}\text{(3)}& \forall x_{n+\mathrm{1}}\ge x_n,F_x\left(x_{n+\mathrm{1}}\right)\ge F_x\left(x_n\right).\end{array}$$

In our study, we selected the slope value as the random variable X, and using a supervised approach, we analysed only the slope values corresponding to the mapped rockfall detachment areas (source area inventory in Rossi et al., 2020) to derive CDF_RSA. Thus, CDF_RSA gives the probability that the slope in rockfall source areas is less than or equal to a given value. This function represents the cumulative probability of the slope causing rockfalls and can be used as a quantitative probabilistic estimation of rockfall detachment for given slope values. The map of source areas obtained using the CDF_RSA approach is a probabilistic map, with values ranging from 0 to 1 for a 0 % probability or a 100 % probability of being a potential rockfall detachment area, respectively. The slope values corresponding to a classification of 1 in the CDF_RSA approach range from 62 to 85°, with a mean slope of 77°. In contrast, the slope values associated with a classification of 0 do not exceed 47.27°, exhibiting a mean slope of 16°.

3.1.3 Probabilistic identification of rockfall source areas using LAND-SUITE

The third method for the source area identification (hereafter referred to as PROB_RSA) proposes a probabilistic modelling framework that applies a combination of multiple multivariate statistical classification models, using the source area locations mapped in the field as the dependent variable and a set of thematic data as independent variables (i.e. morphometric data derived from DEMs and lithological data). The model uses morphometric parameters derived from the digital elevation model and lithological data as an expression of the mechanical behaviour of the rocks as input.

As described in detail in Rossi et al. (2020), we applied the probabilistic framework using LAND-SUITE (LANDslide – Susceptibility Inferential Tool Evaluator), an R-based open-source program (Rossi et al., 2022; (Rossi and Bornaetxea, 2021). The software allowed us to obtain a probabilistic map that expresses the probability of a certain area being a potential rockfall source area. A logistic regression model integrated into the tool was used for the preliminary analysis of different training and validation scenarios to determine whether the model was sensitive to the selection of dependent variables and to identify the best model training configuration for application on the island. Four scenarios were evaluated, incorporating variations in training and validation areas, as well as the inclusion of active-source areas (areas with recent geomorphological evidence of rockfall detachments) and prone areas (geologically and geomorphologically susceptible to rockfalls but lacking recent detachment evidence). The optimal scenario involved model training using data from four fieldwork sites (Sabinosa, El Golfo, Las Playas, and La Estaca), with validation applied to the entire island. This configuration achieved the best performance, with an accuracy of 91.28 % in training and a small difference in validation (2.68 %), as well as an area under the curve (AUC_ROC) of 0.954, the highest among all scenarios. Therefore, the source map obtained using this scenario stands out as the most consistent model, delivering the best performance in island-wide validation.

The final source area zonation was prepared applying a combination of different statistical modelling methods, namely a linear discriminant analysis, a quadratic discriminant analysis, and a logistic regression model. Then, different LAND-SUITE tools were used to evaluate probabilistic source area maps that resulted from different model applications and configurations to verify the modelling performance and to estimate the associated uncertainty. The resulting probabilistic source area zonation was evaluated by integrating the output, expressing the variation for a variety of probability thresholds. Specifically, contingency matrices and plots along with model sensitivity, specificity, Cohen's kappa indices, and ROC curves with the corresponding area under the curve (AUC_ROC) values were used to compare the observed and modelled source areas and to quantitatively explore the performance of different model configurations, allowing the selection of the best model and the corresponding probabilistic source area map. See Rossi et al. (2022) for the details on the training, validation, and combination procedure.

Similarly to the previous identification approach, the map of the source areas obtained using the method implemented in LAND-SUITE is a probabilistic map, with values ranging from 0 to 1 for a 0 % probability or a 100 % probability of being a potential rockfall detachment area, respectively.

3.2 Deterministic rockfall runout modelling

The rockfall runout simulation is a core analysis in rockfall modelling. On El Hierro Island, it was performed using a physics-based model employing the maps described above as input source areas (Fig. 3a, b, c). Such a model is based on the fundamental principles of mass and energy conservation and is employed extensively worldwide to study rockfall runout. In this study, we used STONE, a distributed three-dimensional software based on physics-based simulations. The software is raster based and applies a lumped mass approach to simulate boulder movement along the topography described by a digital elevation model (Guzzetti et al., 2002). The software requires four main inputs: (i) a digital elevation model, (ii) three coefficients maps (i.e. dynamic rolling friction, normal energy restitution, and tangential energy restitution) that simulate the energy loss of a boulder when rolling and bouncing at impact points, (iii) a map portraying the location of the rockfall source areas, and (iv) a map of the number of simulations to be run during modelling (Table 1).

Figure 3 The column on the left shows the maps of the source areas identified using the three different approaches (a, ST_RSA; b, CDF_RSA; and c, PROB_RSA), and the column on the right shows the cumulative counts of rockfall trajectories for each map (d, e, f). See Table 2 for the pixel counts of each map of the source areas.

Table 1 The table shows values of the coefficients (i.e. dynamic rolling friction, normal energy restitution, and tangential energy restitution) used in the rockfall modelling.

The three maps of the coefficients were estimated by considering the different lithological and geotechnical categories reported on the geotechnical map of El Hierro and by selecting values reported for similar lithologies in the literature (Alvioli et al., 2021; Guzzetti et al., 2003; Mateos et al., 2016; Sarro et al., 2020).

The number of simulations run for each source area pixel was obtained by multiplying the binary (i.e. 0 or 1) or probabilistic (i.e. from 0 to 1) value of the source area maps by 10, successively rounded to the closest integer value.

The main output of the runout modelling computed for the three source area maps is the cumulative count of rockfall trajectories (Fig. 3d, e, f).

3.3 Classification of rockfall runout for susceptibility zonation, model comparison, and validation

The map of the rockfall trajectory counts estimates the possibility of a specific pixel being impacted by a rockfall. To derive rockfall susceptibility maps, the trajectory values can be classified using different systems, including equal interval, natural break, quantile, standard deviation, head/tail breaks, and landslide percentage (Alqadhi et al., 2022; Baeza et al., 2016; Cantarino et al., 2019; Tehrani et al., 2022; Wang et al., 2016) in order to make a qualitative interpretation of the results.

To generate a probabilistic susceptibility map, we employed two classification approaches based on the ECDF of trajectory counts and considering an unsupervised and a supervised method.

The unsupervised classification technique is based exclusively on the raster map of rockfall trajectory counts. This method classifies the map by utilizing the ECDF derived from the values of counts obtained in the entire study area with the rockfall runout model (i.e. cells with count values greater than or equal to 1). The resulting map presents values ranging from 0 to 1, representing a probabilistic estimate of the likelihood of each pixel being affected by a rockfall event. Consequently, pixels equal to 1 indicate areas where the susceptibility model predicts the highest probability of rockfall occurrence.

The supervised classification method works similarly, but in this case, the ECDF analysis considers only the count of the trajectories together with the rockfall deposits and/or rockfall talus mapped in the study area. The rockfall deposit mapping can be affected by uncertainty and, to be reliable, should be statistically representative of the different geo-environmental settings controlling rockfall occurrence and evolution.

This two-fold classification methodology was applied to the maps of trajectory counts obtained by STONE using the three maps of source areas as input (i.e. ST_RSA, CDF_RSA and PROB_RSA). As a result, we obtained six ECDF graphs and six susceptibility maps that we compared and analysed using different analyses. The six susceptibility maps were evaluated in pairs considering the three source area maps and the two classification methods. To investigate and quantify the diversity, we used maps of the differences and histograms that enabled the identification of the locations where the susceptibility maps show a higher (or a lower) likelihood of rockfall occurrence. Additionally, two-dimensional (2D) hexagonal bin count heat maps derived from the different couplings of susceptibility maps were plotted to show the correlations between the model outcomes. Hexagonal binning for map comparison is a technique used in data visualization, particularly when dealing with large datasets in 2D scatter plots. It groups data points into hexagonal bins (rather than traditional square bins) to provide a more structured view of the data distribution. The hexagonal shape is often preferred because it avoids the visual artefacts that can result from aligning data into rectangular grids and provides a more compact and efficient way of packing data points (Wickham, 2016).

To validate the models, we used two rockfall inventories: (i) a polygon-type inventory with zones reached by rockfall boulders and zones without any significant evidence of potential boulders having reached them (i.e. red and green polygons in Fig. 1) and (ii) a point-type inventory with the locations of isolated rockfall boulders at their final reaches after runout (i.e. silent witnesses; black dots in Fig. 1c). We first used the polygon-type inventory to derive ROC plots (Rossi et al., 2010, 2022; Rossi and Reichenbach, 2016) and the corresponding area under the curve (AUC_ROC) in order to show the differences between the modelled and observed susceptibility values and provide quantitative estimates of the final rockfall susceptibility zonation performance, regardless of the classification approach adopted. Successively, we analysed the distribution of average susceptibility values (i.e. violin plots and boxplots) within circular buffers of different sizes built around the boulder locations reported in the point-type inventory to verify the ability of the models to discern susceptible conditions corresponding to and in the vicinities of mapped rockfall boulders. Different buffer sizes allow us to consider uncertainty due to local conditions and boulder locations. In the proposed approach, the location of the mapped boulders is used to evaluate the rockfall susceptibility zonation. Commonly, this information is used to evaluate runout models, verifying whether simulations reach the boulder locations entirely or partially. The violin plots show the distributions of the susceptibility data and specifically their probability density, and together with boxplots, they help visualize summary data statistics such as median values and interquartile ranges.

4 Results

4.1 Comparison of different maps of the source area

Following the steps of the methodology, we first compared the maps of the source areas prepared using three different approaches (see Sect. 3.1), which cover the entire island with consistent and equal spatial coverage.

For the slope thresholding approach (ST_RSA), we determined a threshold of 40° by combining geomorphological data, geological analysis, and historical rockfall events. In this case, for the entire island, a total of 727 603 px were identified as being prone to rockfall detachments, corresponding to 18.19 km² (6.8 % of the island, Table 2). To carry out the rockfall simulation, the binary map was multiplied by 10, resulting in two distinct values: 10 simulations corresponding to the rockfall source areas and 0 elsewhere.

Table 2 The table shows the spatial extension of the source areas identified by the three approaches (i.e. ST_RSA, CDF_RSA, and PROB_RSA).

In the second approach, we used CDF_RSA to obtain a probabilistic source map with values ranging from 0 to 1, respectively, for a 0 % probability or a 100 % probability of the area being a potential rockfall detachment area. Unlike the binary values in the ST_RSA map, this probabilistic information allows us to identify the source areas with different levels of certainty. The map shows that 1 628 048 px have non-zero probability of being a potential detachment area, twice the number of pixels identified using the slope thresholding approach (ST_RSA). Source areas identified through CDF_RSA cover a total area of 40.70 km², around 15 % of the island's surface. In this case, the map of the number of runout simulations has integer values ranging from 0 to 10.

The third source area map obtained using the PROB_RSA method shows a total of 3 339 686 px with a non-zero probability of being a potential detachment area, which is equivalent to 84.99 km² or approximately 31.6 % of the entire island surface. Similarly to the CDF_RSA case, the resulting map of the number of simulations has integer values ranging from 0 to 10.

The comparison of source areas identified using the three methods was performed using a spatial overlay in raster format and using frequency-based criteria. The three maps show a diversified spatial arrangement: a total of 727 423 px were recognized as source areas through the three different methods, with the matching areas mostly located on steep slopes (Fig. 4). No pixels were identified as a source area by ST_RSA alone, as they were always associated with either CDF_RSA or PROB_RSA. The number of pixels identified by PROB_RSA alone is 1 855 918, corresponding to more than 55 % of the pixels identified using other methods or combinations of methods (Table 3).

Figure 4 The map shows the spatial comparison of the source areas identified using the three different approaches (i.e. ST_RSA, CDF_RSA, and PROB_RSA).

Table 3 This table shows the differences in the spatial distribution of source areas identified by the three approaches (i.e. ST_RSA, CDF_RSA, and PROB_RSA).

The largest RSA match is observed between CDF_RSA and PROB_RSA, 816 278 px (20.40 km²), while the largest mismatch was for ST_RSA and PROB_RSA, with a deviation of 2 672 196 px (66.80 km²) detected by PROB_RSA but not by ST_RSA. This provides evidence that the PROB_RSA tends to identify a larger number of source areas, covering a larger portion of the study area (1 855 918 px and 46.39 km²).

An additional analysis to evaluate the possible relation with the geotechnical classes revealed that only ST_RSA is able to identify source areas in regions with soft and hard soils.

4.2 Comparison of rockfall simulation and susceptibility maps

The output of the runout simulations (Fig. 3d, e, f) shows diverse spatial distributions of rockfall trajectory counts, providing potentially different information on the susceptibility to rockfalls. To obtain comparable rockfall susceptibility maps, we classified the trajectory count maps using unsupervised and supervised ECDF analyses (Figs. 5 and 6). The application of the ECDFs to the relative trajectory count maps allows us to derive the six probabilistic susceptibility maps shown in Fig. 5. The figure reveals evident differences between the maps derived from the unsupervised ECDFs (Fig. 5a, b, c) that are reduced when considering the supervised alternatives (Fig. 5d, e, f).

Figure 5 Probabilistic susceptibility maps derived from the application of unsupervised (a, b, c) and supervised (d, e, f) ECDF functions (Fig. 6).

Different plot representations were used to compare the six maps and to understand their differences. Figure 6 shows the unsupervised and supervised ECDF functions derived from the outputs obtained using the three source area maps. The unsupervised distributions show larger ranges and a higher number of cells with low trajectory counts (i.e. values close to 0). Additionally, the comparison of the unsupervised ECDFs (Fig. 6a, b, c) reveals a larger number of cells with high count values for ST_RSA followed by CDF_RSA and PROB_RSA; this behaviour is reversed when considering the supervised ECDFs (Fig. 6d, e, f).

Figure 6 Unsupervised (a b, c) and supervised (d, e, f) ECDF functions derived for outputs obtained from the different source area identification methods.

Figures 7 and 8 show the pairwise differences in susceptibility maps obtained using different source area maps and the diversified classification method. Specifically, the figure portrays the following six pairs of results: (a) ST_RSA-unsup − CDF_RSA-unsup, (b) ST_RSA-unsup − PROB_RSA-unsup, (c) CDF_RSA-unsup − PROB_RSA-unsup, (d) ST_RSA-sup − CDF_RSA-sup, (e) ST_RSA-sup − PROB_RSA-sup, and (f) CDF_RSA-sup − PROB_RSA-sup. The lighter colours (i.e. lower absolute difference values) between supervised map pairs and the frequency counts of the corresponding histograms highlight lower differences between the susceptibility outputs obtained by applying supervised ECDFs.

Figure 7 Maps of the pairwise differences in susceptibility maps obtained for different source area identification methods (rows) and diversified classification methods (columns). Negative values indicate a higher probability for the second of the two methods compared.

Figure 8 Histograms of the pairwise differences in susceptibility maps obtained for different source area identification methods (rows) and diversified classification methods (columns).

The 2D hexagonal bin count heat maps (Fig. 9) derived for the different pairs of susceptibility maps confirm these results, showing a better alignment along the bisector of the higher frequency counts obtained for supervised susceptibility maps (Fig. 5d, e, f). These plots are divided into hexagonal bins, and each bin is coloured based on the count of the susceptibility map values. Dark-reddish shades indicate a higher frequency of measurements within the corresponding hexagon, while lighter areas may indicate sparse values.

Figure 9 The 2D hexagonal bin count heat maps derived for the different pairs of susceptibility maps obtained by applying unsupervised (a, b, c) and supervised (d, e, f) ECDF functions. Dark-reddish shades indicate a higher frequency of measurements within the corresponding hexagon.

In addition, the comparison of the trajectory maps with the simplified geotechnical classes (Fig. 1 in Rossi et al., 2020) reveals that the trajectories mainly cross over hard and very hard rocks and only occasionally cross over moderately soft rocks. In the unsupervised maps, very hard rocks are affected by rockfall trajectories approximately 19 %, 25 %, and 42 % of the time, corresponding to ST_RSA, CDF_RSA, and PROB_RSA, whereas for hard rocks, the percentages decrease to 7 %, 17 %, and 37 %. These percentages can be explained by the geological and morphological setting. Furthermore, the hard-soil class shows considerable percentages above 70 %. This distribution can be justified by their position on the lower part of the slopes, where trajectory paths commonly stop. Trajectories do not cross over soft soils, which are mainly located in flat areas. In the supervised maps, the very hard and hard rocks are passed over by the majority of the trajectories (i.e. 81 %, 81 %, and 88 % for ST_RSA, CDF_RSA, and PROB_RSA).

4.3 Rockfall susceptibility model validation

Figure 10 shows the results of the ROC analysis comparing the different susceptibility maps (Fig. 5) and field observations. The graphs show that the model with the best performance is obtained using the PROB_RSA source areas (AUC_ROC = 0.88) followed by the CDF_RSA (AUC_ROC = 0.84), with ST_RSA performing the worst (AUC_ROC = 0.78).

Figure 10 ROC plots and corresponding AUC_ROC values for the six susceptibility maps shown in Fig. 5. The points show values of the hit rate (also referred to as the true-positive rate or sensitivity) and the false-alarm rate (also referred to as the false-positive rate equivalent to 1 – specificity) for a set of probability threshold reference values.

For the same maps, Fig. 11 shows the distributions of the average values within circular buffers of 5, 50, and 100 m defined around the observed boulder locations. Susceptibility median and maximum values increase with the decrease in the buffer size. The distributions of values change significantly for different source areas when the susceptibility is classified using the unsupervised EDCF, whereas they tend to be more homogeneous when the supervised ECDF is applied.

Figure 11 Violin plots and boxplots derived for the average values of susceptibility within buffers defined around rockfall boulder locations. Plots correspond to the six susceptibility maps shown in Fig. 5.

5 Discussion and conclusions

Rockfall modelling is complex and requires a set of dedicated methodological choices and assumptions. Despite specific aspects of the modelling having been largely resolved in the literature (Ding et al., 2023; Noël et al., 2023; Yan et al., 2023; Yang et al., 2021; Žabota et al., 2019), a comprehensive methodology to assess susceptibility to rockfalls is still missing. To fill this gap, we have proposed a workflow that includes methods for the source area identification, the deterministic runout modelling, the classification of runout output to derive objective rockfall probabilistic susceptibility zonation, and finally the comparison and validation of the results. The methodology was applied on El Hierro Island (Canary Islands, Spain), where rockfalls pose a significant threat to structures, infrastructure, and the population. We have presented three methods to identify source areas of increasing complexity, namely ST_RSA, CDF_RSA, and PROB_RSA, which require diversified input. Table 3 and Fig. 4 show how these methods may provide different inputs (i.e. source area and number of simulations) for deterministic rockfall runout modelling, impacting the rockfall trajectory simulations and the corresponding susceptibility zonation (Fig. 5).

To derive probabilistic susceptibility maps, we propose the use of unsupervised and supervised ECDFs of the trajectory counts. We demonstrate with quantitative metrics (Figs. 8 and 9) how the use of the supervised ECDF approach helps to reduce differences and homogenize zonation at the expense of a dedicated mapping effort to derive a rockfall inventory (Fig. 1). This is a significant methodological finding of this work and shows that even when using simple source area identification methods, such as ST_RSA or CDF_RSA, the supervised ECDF application guarantees a reliable and unbiased zonation of rockfall susceptibility. Traditionally, information on rockfall deposits is mainly used to validate the rockfall modelling results. In this study, we also show the relevance of mapping areas not affected by rockfalls to improve the reliability and robustness of the susceptibility zonation. This can be as relevant as the source area mapping and identification. In fact, the application of this workflow demonstrated that such data play a key role in susceptibility zonation classification, preventing overestimation of the results and enhancing their utility for decision-makers.

This study also explores the strategies to validate the rockfall susceptibility outputs using different types of inventory, such as (i) polygon-type maps portraying the zones reached by rockfall boulders and zones without any significant evidence of having been reached by potential boulders and (ii) point-type inventories with the locations of isolated rockfall boulders at the end of the runout (i.e. silent witnesses). Metrics comparing modelled and observed values (i.e. ROC plots and the corresponding AUC_ROC) can be used to show the performances of susceptibility models quantitatively, regardless of the classification approach adopted (Fig. 10).

The ROC analysis reveals differences in the performance of the three source area identification methods. However, identical AUC_ROC values are obtained for unsupervised and supervised ECDFs when the same source area identification method is used. This highlights the fact that the method used to classify the maps of trajectory counts and derive the susceptibility zonation is crucial. The ROC analysis is sensitive to methodological choices and helped us to select PROB_RSA (followed by CDF_RSA and ST_RSA) as the preferred method to identify rockfall source areas. Such results can be explained by the greater statistical robustness of this method (Rossi et al., 2020), which requires dedicated mapping, a set of thematic information, and the use of specific statistical software (Rossi et al., 2022). In general, the results of the multi-criterion techniques used to validate the outcomes and assess the reliability of the susceptibility zonation demonstrate that the larger the effort in the identification of source areas, the more reliable and accurate the rockfall susceptibility zonation. They also highlight the importance of selecting appropriate source area identification methods and incorporating supervised classification to improve rockfall susceptibility zonation. Furthermore, the study highlights the fact that supervised approaches provide added value by fine-tuning the modelling outputs.

Rockfall point-type inventories can be used for a basic verification of the ability of models to discern susceptible conditions corresponding to and in the vicinities of the mapped or observed boulders. This can be performed by analysing the distribution of susceptibility values within circular buffers of different sizes built around boulder locations. Such distributions can be visualized with violin plots (Fig. 11) that show the effects of different classification approaches for rockfall susceptibility zonation. Figure 11 reveals that susceptibility zonation values vary largely within buffers and tend to increase in the vicinity of the boulder locations mapped. Significant distribution differences can be observed among the susceptibility classification approaches and the source area identification criteria. Unsupervised ECDFs (Fig. 11a, b, c) show diversified shapes, with PROB_RSA characterized by more uniform distributions and higher susceptibility values. Conversely, supervised ECDFs (Fig. 11d, e, f) minimize these differences, reshaping the distributions and making them more similar. This means that supervised ECDFs should be preferred because they reduce the effect of the criteria used to identify source areas on the final susceptibility zonation considerably.

In the analysis of rockfall susceptibility at the regional scale, the access to comprehensive data is frequently limited. This constraint impacts the choice of the methodologies employed to define source areas. When only a digital elevation model (DEM) and bibliographic resources are available, the slope thresholding method is preferred. Where additional data, such as geological or geomorphological information, are accessible, investing time in mapping source areas enables the application of probabilistic methods that yield more robust results. Furthermore, maps of trajectory counts are often considered the final modelling outputs; nevertheless we propose performing a supervised analysis to classify them for more reliable susceptibility zonation. Combining probabilistic methods for the source area identification with supervised classification of trajectory counts ensures a more accurate and balanced susceptibility zonation, enhancing its utility for decision-making processes in rockfall hazard management.

Despite the availability of various kinds of software and methods for rockfall runout simulation, we have selected STONE due to its previous use, validation, and application in the study area. Nonetheless, we recognize that the methodological framework proposed in this study remains relevant even when employing other rockfall modelling software. The unsupervised and supervised ECDF analyses are applicable to the trajectory counts generated by any software.

The proposed methodology provides possible guidance for objective and reliable rockfall modelling that is able to support civil protection, emergency authorities, and decision-makers in evaluating and assessing potential rockfall impacts, and it can be a potential strategic support for rockfall warning systems.