Geohazard Plugin: A QGIS Plugin for the Preliminary Analysis of Landslides at Medium-Small Scale

1. Introduction

Protecting individuals and the environment from hydrogeological risks represents a formidable challenge within the framework of slope stability, particularly rapid over these years of extreme weather events and climate change. Increasing temperatures, permafrost thaw, sea level rise, extreme precipitation events, glacier retreat, changing snowpack dynamics, and increased frequency of intense storms have significantly increased the likelihood and intensity of catastrophic events such as landslides and floods [1,2,3]. The risks and damages caused by these catastrophic events are exacerbated by urbanization, infrastructure development, and economic growth in high-risk areas (such as unstable slopes, coastal zones, or flood-prone regions), where the presence of valuable assets (e.g., buildings, infrastructure, and economic resources) and human populations increase their vulnerability to these hazards [4,5,6]. Within this complex landscape of hydrogeological risk, the interplay of geological, morphological, hydraulic, and climatic factors generates a variety of instability phenomena, each with distinct dynamics, evolutionary characteristics, and spatial impacts [7,8].

Effective management and mitigation of these hazards require robust tools for danger identification, susceptibility and hazard analysis, and risk assessment [9,10]. Modern advancements in remote sensing technologies and Geographic Information Systems (GIS) have become indispensable in this domain, offering powerful data analysis, spatial modeling, and visualization capabilities. Techniques such as satellite-based radar interferometry (InSAR), GPS systems, and sensor networks enable continuous and detailed observation of slope movements and targets exposed at risk. In particular, InSAR and its variant, PSInSAR (Permanent Scatterer InSAR), have proven invaluable for monitoring slow landslides over extended periods, providing millimeter-level precision in displacement measurements [11,12,13,14].

A preliminary spatial assessment to identify critical zones is demanded for rapid and unpredictable events like rockfalls and earth flows, characterized by sudden instabilities and high-velocity movements. GIS-based models have emerged as effective tools for hazard and susceptibility analysis, enabling the preliminary evaluation of large areas based on equivalent parameters. These assessments provide critical insights for identifying areas where focusing detailed studies [15,16,17,18]. However, integrating specialized algorithms into user-friendly GIS platforms remains challenging, limiting accessibility for non-expert users and hindering the practical application of advanced geohazard assessment techniques.

The Geohazard plugin was developed as an open-source tool within the QGIS environment to address this gap. The Geohazard plugin incorporates algorithms targeting key aspects of slope instability, including SAR data reliability evaluation, susceptibility analysis for shallow landslides, and spatial hazard assessment for rockfall phenomena. Practical applications and case studies demonstrate how the plugin simplifies complex analyses and supports practical decision-making for more resilient and sustainable landscapes in landslide-prone areas. The plugin facilitates efficient and scalable geohazard assessments across varying spatial scales, providing land-use planners, hazard managers, and decision-makers with actionable insights into geohazard identification.

2. The Geohazard Plugin

The Geohazard plugin [19] incorporates several specialized algorithms to evaluate the reliability of InSAR data in monitoring slow landslide phenomena, to assess slope susceptibility to shallow landslides, and to perform a preliminary assessment of rockfall spatial hazard. Key algorithms included in the Geohazard plugin are the following:

• Groundmotion – C index: this algorithm evaluates the reliability of SAR data for monitoring slope movements by calculating the C index [20], a parameter that quantifies satellite observation accuracy. The result is a map of the representativeness of InSAR data for slow landslide monitoring, useful for accurate data interpretation in terms of magnitude and velocity of the slope displacements;
• Landslide – Shalstab: this algorithm evaluates slope susceptibility to shallow landslide triggering by calculating a critical rainfall infiltration, based on the hydro-mechanical model developed by Montgomery & Dietrich (1994) [21]. The result allows the preliminary identification of areas prone to shallow landslide under different infiltration conditions at a medium/small scale;
• Rockfall – Droka_Basic and Rockfall – Droka_Flow: these two algorithms estimate zones potentially impacted by rockfall events originating from distributed source points on a slope. The estimation is based on the energy line concept and employs two approaches: Droka_Basic utilizes the Cone Method to identify areas that can be reached from one or more rockfall source points [22]; Droka_Flow is based on a hydrological model to simulate rockfall trajectories, resembling the path of a water droplet descending along the steepest gradient. Both algorithms provide a spatial representation of rockfall phenomena, including indications of the invasion zone, block velocities, and kinetic energy levels reached in the exposed areas. Such results can provide a preliminary assessment of rockfall susceptibility and relative spatial hazard at a medium-small scale.

3. Groundmotion – C Index

The SAR is an active remote sensing technique wherein an antenna emits radiation upon reaching the ground surface and is reflected and captured by the same sensor. SAR data can enable the creation of 3D terrain models; their processing allows for monitoring ground deformations and studying landslide phenomena with millimetric precision. Within the Copernicus program [23], data acquired by the Sentinel-1 satellites are freely available through the EGMS (European Ground Motion Service), updated annually. These data support applications such as monitoring the structural integrities of dams, bridges, railways, and buildings. It also allows urban planners to make data-driven decisions regarding infrastructure development by assessing the potential risks posed by natural hazards, such as landslides or subsidence [24].

To facilitate the interpretation of SAR data, the first algorithm within the Geohazard plugin calculates the C index, an indicator representing the percentage of ground movement detectable by a satellite. The index is important for evaluating the capability of the satellite to monitor ground movements and identify deformations or terrain changes. Additionally, a visibility percentage map is produced, which is a function of the position of the satellite to the slope, accounting for the radar acquisition angles. The visibility map helps determine the effectiveness of the satellite in collecting accurate and comprehensive data. Finally, the model produces a visibility class map, classifying the terrain into areas according to their visibility (Table 1). This classification provides a detailed overview of areas that can be monitored more precisely than others, aiding in optimizing satellite data usage for planning and monitoring activities.

SAR displacement measurements are one-dimensional, capturing only the movement component along the LOS (Line Of Sight), which is tilted with respect to the vertical at an angle that varies depending on the satellite and the acquisition geometry. In flat terrains and slopes exposed to the north or south, two-dimensional analysis is possible using ascending and descending passes. However, the system exhibits limited sensitivity along the direction of the sensor orbital. Since the orbits of the main satellites (ERS, Sentinel, RADARSAT) are almost polar, any deformation occurring along the north-south direction results in minimal projection along the LOS. To address this limitation, a correction index C can be applied to each LOS measurement to compute the actual velocity component in the slope direction (V_SLOPE). The C index accounts for satellite-dependent parameters, such as incidence angle θ, i.e., the inclination of the LOS with respect to the vertical, and track angle φ, i.e., the angle between the projection of the satellite trajectory on the horizontal plane and geographic north, as well as topographic parameters, including terrain slope α and orientation β [20,26] (Figure 1).

The track angle and incidence angle are provided within the processed SAR images for each satellite, while slope and aspect values are derived from a DTM. The resolution of the DTM significantly influences the accuracy of the C index. High-resolution DTMs improve the representation of terrain profiles but introduce greater variability in slope and aspect values, impacting C index computation and subsequent V_SLOPE calculations. For balanced accuracy and computational efficiency, a cell resolution of at least 20 m is recommended, consistent with InSAR data resolution [25]. The C index represents the component of movement detectable by the SAR sensor, with a value ranging from -1 to 1.

The C index is a function of the incidence and track angles and the characteristics of the slope and can be calculated as follows:

$$C=\left[N·\mathrm{c}\mathrm{o}\mathrm{s}\left(S\right)·\mathit{sin}\left(A-{\displaystyle \frac{\pi }{2}}\right)\right]+\left[E·(-1)·\mathrm{c}\mathrm{o}\mathrm{s}\left(S\right)·\mathit{cos}\left(A-{\displaystyle \frac{\pi }{2}}\right)\right]+\left[H·\mathrm{s}\mathrm{i}\mathrm{n}\left(S\right)\right]$$

(1)

where S is the slope of the surface; A is the aspect of the surface with respect to the north; N, E, and H are direction cosines [20]. This allows the calculation of V_SLOPE (mm/year), the projection of the LOS velocity (V_LOS) along the direction of maximum slope, using the following relationship [27,28]:

$$V_{SLOPE}={\displaystyle \frac{V_{LOS}}{C}}$$

(2)

The relationship also allows for estimating whether the movement is toward or away from the satellite. Indeed, for values of C close to 1, the satellite accurately observes the movement, and V_SLOPE tends to V_LOS. For negative C values, the movement is recorded in the reverse direction relative to the LOS geometries. In this case, if the velocity components are consistent with the movement along the slope, V_SLOPE will have a sign opposite to V_LOS. For values of C close to zero, V_SLOPE tends toward infinity, significantly amplifying errors in the projection [20]. Furthermore, it should be noted that positive V_LOS values can indicate either potential errors or real movements. When V_SLOPE shows a positive value, it may indicate that a landslide is moving upward rather than downward, an unlikely condition. Therefore, such cases should be assessed individually [20].

The processing model in QGIS generates a percentage visibility map that illustrates proportion of the velocity V_LOS along the direction of maximum slope. The percentage of visibility is calculated using the formula:

$$\%V_{REAL}=abs\left(C\right)·100$$

(3)

This thematic map helps in quantifying the proportion of the terrain visible to the satellite, accounting for the radar acquisition angles. The visibility percentage serves as a critical metric for assessing the capability of the satellite to acquire accurate and comprehensive data. Additionally, the algorithm generates a visibility classification map that categorizes the terrain into five distinct classes based on visibility levels. This classification provides a detailed representation of areas with varying monitoring precision, highlighting regions where data acquisition is more effective.

The calculation of the C index has been implemented in the Geohazard plugin, input data and results are shown in Table 1.

3.1. Validation of the Model

For validation purposes, the results generated by the model on a training dataset were compared with externally performed calculations. The DTM of the training dataset represents a high mountain area of 154x172 m² in Piemonte (Italy). The DTM has a cell size of 25 m to balance accuracy and computational efficiency. An Excel spreadsheet was created to facilitate this comparison, incorporating Equation 1 to derive the C index. This approach enabled a direct comparison between the output of the algorithm and the externally calculated values. A validation process was carried out on nine cells, with three cells selected for each pair of angles. The input data for the C index calculation included aspect, slope, incidence angle θ, and track angle φ, and the corresponding values were assigned to each cell being validated. The comparison of results confirmed the validity and accuracy of the model, as outlined in Table 2.

4. Landslide – Shalstab

Shalstab is a physically-based model used to assess the triggering of shallow landslides due to infiltrated rainfall, following the method proposed by Montgomery & Dietrich (1994) [21] and Dietrich and Montgomery (1998) [29]. The method (Figure 2) enables the analysis of topographic, soil, and rainfall data to evaluate slope stability and solves the relevant equations using a gridded representation of the terrain, combining a limit equilibrium stability model for infinite slopes with a steady-state hydrological model. The main assumptions underlying the formulation are as follows [30,31]:

• infinite slope;
• planar failure surface, parallel to the slope, located at the interface between the shallow layer and the underlying layer of lower permeability;
• Mohr-Coulomb shear strength criterion expressed in terms of effective stresses;
• steady flow parallel to the slope;
• absence of deep drainage and flow within the underlying substrate.

According to the Limit Equilibrium Method, the stability condition is determined by calculating the Factor of Safety (Figure 2a and Equation 4):

$$F_S={\displaystyle \frac{c\text{'}+(\gamma ·z·{cos}^2\theta -u)·tan\phi \text{'}}{\gamma ·z·sin\theta ·cos\theta }}$$

(4)

$$u={\gamma }_w·h·{cos}^2\theta$$

(5)

where θ is the slope inclination, z is the landslide thickness, c’ and $\mathsf{\phi }$’ are cohesion and friction angle respectively, $\mathsf{\gamma }$ is soil unit weight, u is water pressure on the failure surface,${\mathsf{\gamma }}_{\mathrm{w}}$ is water unit weight, and h is the water table level above the failure plane.

There are several approaches to consider the effects of rainfall and its infiltration into the soil. In Shalstab, steady-state conditions are assumed, meaning there is no temporal variation in rainfall and soil saturation. The model is based on the observation that shallow landslides tend to occur in topographic hollows, where shallow subsurface flow convergence leads to increased soil saturation, pore pressure, and reduced shear strength, thus triggering potential failure. The hydrological component of the model is simplified by assuming that rainfall infiltrates the soil and is routed downslope, following the topographic gradient. This implies that the contributing area at any given point is defined by the specific catchment area derived from the surface topography (Figure 2b).

For a detailed analytical discussion of the hydrological model employed in Shalstab, reference can be made to [29]. The following section outlines the hydrologic component of this model that, coupled with the stability component (Equations 4 and 5) leads to the implementation of the F_S equation (Equation 6):

$$F_S={\displaystyle \frac{c\text{'}+(-W·{\gamma }_w)·z·{cos}^2\theta ·tan\phi \text{'}}{\gamma ·z·sin\theta ·cos\theta }}$$

(6)

The instability condition (F_S ≤ 1) for each cell in the domain is expressed by (Equation 7):

$${\displaystyle \frac{a}{b}}\ge \left[{\displaystyle \frac{c^{\text{'}}}{{\gamma }_w·z·sin\theta ·cos\theta }}+{\displaystyle \frac{\gamma }{{\gamma }_w}}·\left(1-{\displaystyle \frac{tan\theta }{tan\phi \text{'}}}\right)\right]·{\displaystyle \frac{T}{q}}·sin\theta$$

(7)

Where, according to Figure 2b:

• a is the contributing upslope area draining across the contour length b;
• $T=K_s·z·cos\theta$and K_s is the permeability coefficient of saturated soil;
• $W={\displaystyle \frac{h}{z}}={\displaystyle \frac{q}{T}}·{\displaystyle \frac{a}{b·sin\theta }}$ and q is effective rainfall.
• This condition is evaluated only for cells that do not meet the conditions of absolute stability and absolute instability, given by:

$$tan\theta \le \left(1-{\displaystyle \frac{{\gamma }_w}{\gamma }}\right)·tan{\phi }^{\text{'}} (\mathrm{a}\mathrm{b}\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{e} \mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y})$$

(8)

$$tan\theta \ge {\displaystyle \frac{c\text{'}}{\gamma ·z·{cos}^2\theta }}+tan{\phi }^{\text{'}} (\mathrm{a}\mathrm{b}\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{e} \mathrm{i}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{b}\mathrm{i}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y})$$

(9)

The condition of absolute stability identifies cells classified as stable even when the soil is fully saturated. The condition of absolute instability indicates cells that are unstable even in the absence of infiltration.

Furthermore, by imposing the limit equilibrium condition (F_S = 1) and solving Equation 7 with respect to q, the critical infiltration q_cr threshold can be determined (Equation 10):

$$q_{cr}=\left[{\displaystyle \frac{c^{\text{'}}}{{\gamma }_w·z·sin\theta ·cos\theta }}+{\displaystyle \frac{\gamma }{{\gamma }_w}}·\left(1-{\displaystyle \frac{tan\theta }{tan\phi \text{'}}}\right)\right]·{\displaystyle \frac{a}{b}}·T$$

(10)

If the information on the presence of roots in the shallow soil layers is available, their resistance contribution can be incorporated by increasing the soil cohesion, which will then be composed by the sum of the contribution of soil c_soil and roots c_root.

Based on the aforementioned method, the Shalstab algorithm provides the spatial distribution of the absolute stable and absolute unstable cells, as well as the critical effective (infiltrated) rainfall that leads each cell to a F_S of 1. The required input data are listed in Table 3. Although most parameters vary spatially across the slope, to assist the user in preparing a quick analysis, the Geohazard plugin includes an additional tool called Shalstab input raster creator, which allows for the generation of constant-value input raster files from the DTM.

The results of the analysis are given as 4 raster maps, showing: 1) absolute stable zones; 2) absolute unstable zones; 3) distribution of critical infiltrated rainfall q_cr (mm/day); 4) susceptibility classes. In the q_cr map, absolute stable cells assume no data values and absolute unstable cells assume q_cr = 0. The susceptibility raster map is a reclassification of q_cr in 7 classes (Table 4) as suggested by Dietrich and Montgomery (1998) [29]. This classification provides the value of q_cr/T in logarithmic form. This ratio reflects the magnitude of the water infiltration, represented by q, concerning the ability of the surface to move the water downslope, i.e., the transmissivity (T). The larger the value of q relative to T, the more likely the ground is to become saturated. This is a useful representation for identifying areas of the study domain with different behaviour and it can be interpreted as a susceptibility map since it spatially classifies the distribution of instability conditions. It can be converted into a hazard map if an infiltration model is applied, to compute a critical triggering gross rainfall P_cr from q_cr. Given P_cr, the probability of occurrence of the related meteorological event can be made explicit as a return period and the time dimension of hazard is achieved [32].

4.1. Validation of the Model

For the validation of the model, it was applied to the 10 m resolution DTM of the Dego area, south Piemonte (Italian National Geological Mapping Project, sheet 211), which was used to create the landslide hazard map for the Carg Project [33]. In the Carg Project, absolute hazard was defined, including the temporal aspect through the return period of critical rainfall and the assessment of the destructive potential of the phenomenon (intensity). To validate the model, input data derived from the Carg Project were used (Table 5), and the results obtained in terms of absolute stable, absolute unstable cells, and distribution of q_cr were compared.

In Figure 3, the resulting critical rainfall map (q_cr) is shown after applying the Shalstab model to the Dego area. Twenty representative cells in different areas were considered for comparison (Table 6): the values of critical infiltration vary with negligible differences, due to the use of different hydrology calculation algorithms. However, this comparison highlights the accuracy of the obtained results, and the validation can be considered successful.

In addition, to evaluate the influence of DTM cell size on the analysis results, a calculation was repeated on the same area using the 5 m DTM provided by Regione Piemonte. Again, the parameters listed in Table 5 were used. In Table 7 the critical rainfall q_cr is reported for the 20 representative points referred to in the 10 m and 5 m DTM analyses. The susceptibility classes are also reported for comparison. The results show that the difference in q_cr calculation can be not negligible. This is probably due to the fact that DTM resolution affects both the mechanical (Fs) and the hydrological component of the model, resulting in different inclinations of the slope in each cell and thus in different catchment area calculations. This aspect of the problem, that affects the susceptibility classification of the area, will be further analysed in the future.

5. Rockfall – Droka

The Droka model is devoted to the simplified analysis of rockfall runout from diffuse source points across a slope. The analysis is based on the concept of energy line, the line connecting a rockfall source point to the farthest stopping point of a block detached from that source. The energy line is inclined at an angle $\mathsf{\phi }$_p with respect to the horizontal plane (i.e., the energy angle), which can be interpreted as a global block-slope friction angle. This angle is representative of all the dissipative phenomena occurring throughout the motion phases of the block along the slope. The concept of energy line is implemented in the Geohazard plugin through two independent yet complementary algorithms: Droka_Basic and Droka_Flow. Both simulations provide results in the form of point vectors, distributed across the invasion zone. For each point, the following information is provided: minimum, maximum, and mean block velocity [m/s]; minimum, maximum, and mean kinetic energy [kJ]; and the number of simulations that affected a specific point (frequency).

The two models use different methods to identify the areas involved in the phenomenon and can be used independently or combined to obtain complementary information. A detailed description of the algorithms is provided below along with a discussion of their potential and limitations through practical application examples.

5.1. Rockall - Droka_Basic

Droka_Basic module is based on the Cone method [10,18,22], which enables the definition of the area affected by rockfall events on a slope. This is achieved through a cone representing the envelope of all possible rockfall paths. The apex of the cone is located at the source point and its geometry is defined by three angles (Figure 4a): the dip direction angle θ of the slope at the source point defines the orientation of the cone; the energy angle $\mathsf{\phi }$_p outlines its vertical extension and the lateral spreading angle α delineates the extension of the cone in the horizontal plane. The intersection of the cone with the surface of the slope defines the rockfall-affected area, where each point within this area may potentially be reached by a rockfall path originating from the source. When multiple source points are considered, the module can delineate zones subjected to rockfalls through the overlapping of simulated cones. This approach enables the creation of a frequency map, where areas with higher frequencies are the most prone to be reached by fallen blocks, identifying them as the most susceptible to this type of danger. Using the concept of the energy line, the module enables the calculation of the velocity and kinetic energy of the block at any position within the affected area [18]:

$$v(x,y)=\sqrt{2g·\left[H\left(x_0,y_0\right)-h\left(x,y\right)-\sqrt{{(x-x_0)}^2+{(y-y_0)}^2}·tan{\phi }_p\right]}$$

(11)

where g is gravity acceleration, (x₀,y₀) are the source point coordinates, H(x₀,y₀) is the elevation of the source point, h(x,y) is the elevation of the topographic surface in the point (x,y), $\mathsf{\phi }$_p is the energy angle.

Therefore, the kinetic energy E can be estimated at any position given the block mass m [18]:

$$E(x,y)={\displaystyle \frac{1}{2}}·m·v^2(x,y)$$

(12)

The Cone Method in the Droka_Basic module is implemented by constructing a triangle with its apex at each source point, characterized by specified input parameters. A point layer is then generated, containing one point for each DTM cell covered by the triangle. For each point, the height difference ΔH between the ground and the energy line is computed (ΔH=h_k(x,y), Figure 4a). Points showing a positive ΔH value are identified as part of the invasion area. Starting from ΔH, block velocity and kinetic energy are computed for each point according to Equations 11 and 12.

The input data required by the Droka_Basic module are described in Table 8. The attributes of source points are obtained from the DTM while the cone parameters are constant across all source points and are specific to each case study. Calibrating these parameters is essential to achieve reliable and conservative results. It should be based on information from past events in the area, rockfall inventories, territorial databases, and quick field observations of the slope. The output is a vector layer of points located at the center of the DTM cells within the affected area, containing the data listed in Table 8.

5.2. Rockfall – Droka_Flow

The Droka_Flow module employs a hydrological approach to simulate potential rockfall paths by tracing the line of maximum slope for each source point. The path terminates when the inclination of the line passing through the source point is lower than the inclination of the energy line, defined by the user.

Each analysis includes five simulations from each source point. These simulations adjust DTM elevation using a Monte Carlo approach within a user-defined standard deviation around the original DTM elevation values, assuming a simplified normal distribution of elevation for each cell. The mean corresponds to the original DTM height, while the user-defined standard deviation is constant across cells. This approach provides a simplified representation of the variability in rockfall trajectories due to factors like ground impact, bouncing, and obstacles (e.g., rocks, trees, or structures).

The Droka_Flow module relies on the GRASS GIS function r.drain to calculate five paths from each source point. As with Droka_Basic, the analysis result is a vector point layer is then generated, containing one point for each DTM cell covered by the simulated paths. For each point, the height difference ΔH between the energy line and the ground is calculated, leading to block velocity and kinetic energy calculations when ΔH is positive (Equations 11 and 12; Figure 4b). The paths are then merged for the statistical analysis of block velocity and kinetic energies.

The input data required by the Droka_Flow module are described in Table 9. Again, source points are obtained from the DTM while the energy angle and the standard deviation of the distribution of cells elevation are specific to each case study. In particular, the standard deviation should be related to the DTM size, as discussed in Section 6. The output is a vector layer with point attributes similar to those described in Droka_Basic and listed in Table 8. Frequency data is omitted here due to the limited number of paths per source point, so energy count and percent fields are not included.

5.3. Validation of the Modules

To validate the results produced by the algorithms implemented in Droka_Basic and Droka_Flow modules, an initial set of analyses was carried out on a synthetic slope DTM with 5 m cell size. The slope was generated by joining two planar surfaces: a runout plane inclined at 60° to the horizontal and a nearly horizontal plane with an inclination of approximately 2°. A single rockfall source point was located at an elevation of approximately 500 m, with a mass of 2600 kg and a volume of 1 m³. Parametric analyses were performed on this synthetic slope to assess the influence of key parameters on the results, using both Droka_Basic and Droka_Flow. The best parameters fit for both Droka modules was found to be lateral spreading angle α equal to 10°, and standard deviation equal to 0.1 m. The results are consistent for the modules, as reported in Figure 5, and align with the simplified analytical model in the calculation of block velocity and kinetic energy, following Equations 11 and 12.

However, for natural slopes, where both aspect and inclination vary widely, rockfall trajectories may spread much more. In addition, natural slopes feature niches, incisions, or debris accumulations due to erosion and instability phenomena over time, which significantly influence the rockfall paths. Therefore, it is suggested to calibrate the parameters for susceptibility analyses to the specific site conditions. A proper calibration can be achieved through the back analysis of actual events, using data gathered from in situ observations, as described in Section 6.

6. Application of Rockfall – Droka Modules: Bobbio Pellice Case Study

This section presents an example of the application of the Rockfall - Droka modules for a preliminary assessment of rockfall susceptibility and relative spatial hazard, illustrating the recommended procedure and discussing the results. The example focuses on a steep slope near the village of Bobbio Pellice, Piemonte (Western Alps), where recurrent rockfall events threaten a cluster of alpine buildings. On 24 June 2024, a large rockfall event destroyed one of these buildings, creating high-risk conditions.

6.1. Parameter Calibration

The Bobbio Pellice rockfall event involved a slope sector impacted by similar events in the past and was triggered by intense meteorological conditions. Observations from the survey conducted by GeoAlpi Consulting (2024) [34] include the following:

• the detachment occurred at approximately 1,940 m a.s.l;
• a block of approximately 30 m³ struck and destroyed a building located at about 1,750 m a.s.l., with visible traces along its trajectory (Figures 6a,c);
• this path was influenced by other blocks and debris on the slope, displaced during all the collapses;
• additional smaller blocks of approximately 1-2 m³ reached the bottom of the valley, with clearly defined trajectories. The stopping position of some of these blocks is shown in Figure 6b.

Based on the available information a back-analysis was performed; three source points were selected within the steepest area of the involved zone (Figure 6b). A block with a volume of 30 m³, mass of 78,000 kg and unit weight of 2600 kN/m³ was used to estimate kinetic energy. The 5 m DTM provided by Regione Piemonte (2010) was employed for an initial set of analyses to calibrate parameters and compare the results from Droka_Basic and Droka_Flow.

For Droka_Basic, the optimal results in terms of invasion area were achieved with a lateral spreading angle α of 22° and an energy angle φ_p of 37°, as shown in Figure 7, where the distribution of block kinetic energy within the invasion area is also included. The computed kinetic energy at the impact point with the building is approximately 5,000 kJ, consistent with the total destruction of the building.

To calibrate the Droka_Flow model, the influence of the assumed standard deviation was analyzed in this study through a parametric analysis, with standard deviation values ranging from 0.1m to 5m (2% to 100% of DTM cell size). Each simulation was repeated 5 times to generate about 20 paths (Figure 8). The results suggest that the most conservative simulations occur with standard deviations of 0.5 m and 1 m, corresponding to 10% and 20% of the 5 m DTM cell size, respectively. In these cases, at least one trajectory reaches the affected building, with kinetic energy around 5,000 kJ, consistent with the results from Droka_Basic. As shown in Figure 8, the simulated paths do not encompass the entire invasion area (in white), as only trajectories in the direction of maximum slope are simulated by the r.drain function. Nonetheless, this analysis demonstrates that the location of the building is the most prone to being impacted by a block detached from the source area, indicating it as the most critical area. The combination of Droka_Flow and Droka_Basic provides dual insights: first, it identifies the entire area potentially impacted by the rockfall, accounting for all possible trajectories, including the less probable ones; second, it highlights the most critical paths within this area.

The back-analyses with Droka_Basic and Droka_Flow modules were repeated with 1 m, 5 m and 10 m DTM resolutions for evaluating the influence of DTM resolution on the results. Figure 9 shows the results of Droka_Basic for the case with φ_p = 37° and α = 22°. It is worth noting that changes in DTM discretization slightly affect the orientation of source points. Specifically, the 1 m DTM results in greater variability in dip direction at source points compared to the 5 m and 10 m DTMs, leading to a wider invasion area. Despite this variability, a lateral spreading angle α of 22° remains the best calibration to cover the invasion area conservatively. In contrast, with the 10 m DTM, a conservative outcome requires increasing the α value. The optimal φ_p value for delineating the rockfall-affected area (37°) remains consistent across DTMs, suggesting it is a characteristic parameter of the phenomenon, dependent on the slope geometry and block properties rather than on DTM resolution.

Figure 10 presents the results of Droka_Flow for the case with φ_p = 37° and a standard deviation of 10% of the DTM size (i.e., 0.1, 0.5, and 1 m for DTM sizes of 1, 5, and 10 m, respectively), which yielded the most conservative results. This comparison shows that the 1 m DTM results in a very limited trajectory spread around the maximum slope path (this result is obtained regardless of the chosen standard deviation). The widest trajectory spread occurs with the 10 m DTM; however, many trajectories fall outside the designated rockfall invasion zone (shown in white in Figure 10). The analysis that better simulates the actual invasion area seems to be with the 5 m DTM, for both Rockall - Droka modules.

These observations suggest that effective parameter calibration for susceptibility and hazard analysis can be achieved through back-analysis conducted within the same geomorphological context and using the same DTM resolution, ensuring consistent results. As discussed, Droka_Basic and Droka_Flow modules complement each other, with their comparison highlighting both the most significant trajectories and the full invasion area. This combined approach offers valuable insights for defining rockfall susceptibility.

6.2. Susceptibility and Relative (Spatial) Hazard Analysis

To exemplify the suggested procedure for susceptibility analysis, this study was extended to the slope surrounding the area where the buildings are located. Since no information was available to quickly identify the possible source areas, all DTM cells with an inclination greater than a specific threshold were selected on the slope above the buildings. Following [35], the resulting threshold is 49° for a 5 m DTM. A total of 80 source points was randomly selected within cells exceeding this threshold, covering a slope sector approximately 200 m in length, as shown in Figure 11a. Finally, a block with a volume of 1m³ was used to estimate the kinetic energy.

Figure 11b compares the results of the analyses conducted with Rockfall - Droka modules in terms of runout area, while Figures 11c and 11d show the mean kinetic energy obtained by the two modules. To assess rockfall susceptibility, the Droka_Basic spatial frequency results, defined as the number of overlapping cones out of the total source points (Figure 12a), were multiplied by the mean kinetic energy values (Figure 11c) at each point. The final output, shown in Figure 12b, highlights the most critical zones by combining both spatial frequency and phenomenon intensity and can be interpreted as a relative (spatial) hazard.

The comparison between the results from Droka_Basic and Droka_Flow again confirms the consistency of the two algorithms. Droka_Basic is more conservative regarding total invasion area (Figure 11b) and mean kinetic energy (Figures 11c and 11d), as the area covered by rockfall paths simulated by Droka_Flow (Figure 11b) is concentrated in the northern part of the invasion zone and shows lower mean kinetic energy values in the central region. However, Droka_Basic indicates the highest frequency in this central zone (Figure 12a), suggesting that blocks from multiple source points can reach it. Differences in the mean kinetic energy calculated by the two algorithms arise from statistical comparisons based on different sample populations, although the maximum kinetic energy values computed are identical across the entire area. Regarding slope susceptibility and relative (spatial) hazard, the results from both algorithms are again consistent, highlighting the critical situation for certain buildings. This finding underscores the need for a more detailed analysis to accurately assess the risk to these structures and to establish appropriate mitigation measures.

7. Conclusions

In this paper, the Geohazard plugin is presented. The plugin is developed as an open-source tool within the QGIS environment and it offers a useful and user-friendly framework for assessing geohazards, related to landslides and rockfalls. By integrating several algorithms, it addresses challenges associated with complex geomorphological conditions, providing actionable insights for hazard management and mitigation at medium-small scale.

The Groundmotion – C Index module [20] demonstrated its reliability during validation on a training dataset with a 25 m cell size DTM, where its results closely matched external calculations. The analysis improved the interpretation of InSAR data for slow-moving landslides.

The Landslide - Shalstab module effectively evaluated the susceptibility to shallow landslides by calculating the critical rainfall infiltration required to trigger instability. Validation on the Dego area (located in the hilly Langhe region, north-western Italy), using a 10 m cell size DTM and consistent input parameters confirmed the ability of the model to delineate areas prone to shallow landslides, with results aligning well with hazard maps developed in previous studies [30]. The critical rainfall map (q_cr) and derived susceptibility classifications provided valuable insights into the spatial distribution of instability conditions, which can be further refined into hazard assessments through rainfall and infiltration modeling.

The Rockfall - Droka_Basic and Droka_Flow modules for rockfall analysis were applied to the Bobbio Pellice case study, an area in the Italian Western Alps [34], to assess rockfall affected zones and relative (spatial) hazard. Parameter calibration based on a well-documented rockfall event in the Bobbio Pellice area enabled accurate delineation of invasion zones and identification of critical trajectories. Droka_Basic captured the overall invasion zone using energy and lateral spreading angles, while Droka_Flow simulated individual trajectories, both providing preliminary estimates of block velocities and kinetic energy. Results showed that the destroyed building in the area was located in a critical zone, with estimated block kinetic energy exceeding 5,000 kJ.

These findings illustrate the reliability of the Geohazard plugin in addressing diverse geohazard challenges. The combination of algorithms provides a toolkit for geohazard assessment, enhancing the understanding of slope instability phenomena and supporting decision-makers in planning and implementing effective risk mitigation strategies.

In a time of increasing hydrogeological instability driven by climate change, the need for effective prevention strategies has become a critical priority for protecting both people and infrastructure from hazardous events. Prevention, as a key aspect of climate change adaptation, requires proactive measures to reduce the impacts of geohazards. This work aims to support and enhance prevention strategies, which are essential in the context of vulnerable territories such as Italy, where almost 94% of the Italian municipalities are at risk of landslides, floods, and coastal erosion, with over 8 million people living in high-hazard areas [36].