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

Abstract

Landslides are a major global threat, endangering lives, infrastructure, and economies. This paper introduces the Geohazard plugin, an open-source tool for QGIS, designed to support medium–small-scale landslide analysis and management. The plugin integrates several algorithms, including the Groundmotion–C index for evaluating SAR data reliability, Landslide–Shalstab for assessing shallow landslide susceptibility, and Rockfall–Droka for estimating rockfall invasion areas and the rockfall relative (spatial) hazard. An application example is provided for each module to facilitate validation and discussion. A case study from the Western Italian Alps highlights the practical application of the Rockfall–Droka modules, showcasing their potential to identify critical zones by integrating the results on affected areas, process intensity, and preferential paths. Emphasis is given to the calibration of model parameters, a critical aspect of the analysis, achieved through a back-analysis of a rockfall event that occurred in June 2024. The Geohazard plugin streamlines geohazard assessments, providing land managers with actionable insights for decision-making and risk mitigation strategies. This user-friendly GIS tool contributes to enhancing resilience in landslide-prone regions.

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 an 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].

The effective management and mitigation of landslide risk require robust tools for danger identification, susceptibility and hazard analysis, and landslide monitoring [9,10]. Modern advancements in remote sensing technologies and Geographic Information Systems (GISs) 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 the continuous and detailed observation of slope movements and targets exposed at risk. In particular, InSAR and its variant, Permanent Scatterer InSAR (PSInSAR), have proven invaluable for monitoring slow landslides over extended periods, providing millimeter-level precision in displacement measurements [11,12,13,14].

The rate of movement is a key parameter in landslide management, significantly influencing monitoring strategies and risk assessment outcomes. Large, slow movements can be effectively monitored using InSAR techniques, which provide valuable datasets for defining potential landslide scenarios and analyzing landslide behavior. These insights are essential for designing appropriate risk mitigation measures [15]. Conversely, rapid landslides, characterized by sudden instabilities and velocities often exceeding 5 m/s [16], present significant challenges for real-time monitoring. These events are typically widespread, recurrent across extensive areas, and difficult to predict, further complicating data collection and hazard assessment [17]. Thus, rapid and unpredictable events such as rockfalls and earth flows demand preliminary spatial assessments to identify critical zones that warrant detailed studies [18,19,20,21]. In this context, GIS-based models have become indispensable tools for territorial analysis, enabling large-scale evaluations based on equivalent parameters. Numerous GIS tools and plugins have been developed for susceptibility and hazard assessments of various landslide types, addressing different aspects.

Regarding slope susceptibility to landslides, models like Stability Index Mapping (SINMAP) [22] evaluate the slope stability by integrating hydrological and geotechnical factors. Similarly, Shallow Landslide Stability (SHALSTAB) [23] integrates topographic attributes with simplified mechanical and hydrological parameters to estimate the shallow landslide potential in relation to rainfall infiltration. To address uncertainties in physically based landslide susceptibility analysis, especially in seismic regions, First-Order Reliability Method (GIS-FORM) provides a user-friendly GIS extension [24]. Meanwhile, Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability (TRIGRS) [25] incorporates transient hydrological inputs in the stability analysis of slopes. The QPROTO plugin [10,21], designed for the QGIS environment, allows preliminary rockfall hazard assessments by estimating the potential invasion area and process intensity based on a set of source points and some equivalent parameters. Finally, open-source GIS tools such as QGIS and GRASS GIS support modular geohazard assessments. For instance, the GRASS GIS module r.slope.stability [26] enables physically-based slope stability analyses over extensive areas, integrating datasets, geotechnical parameters, and elevation models. This module provides a flexible framework for assessing the slope failure probability or factor of safety using limit equilibrium methods. Additionally, r.slopeunits facilitates reproducible terrain subdivision and classification, enhancing the landslide modelling reliability [27].

However, despite the widespread availability of these methods, their applicability is often restricted to a partial evaluation of the problem. Additionally, their computational intensity and reliance on extensive input data frequently limit their usability, particularly for non-expert users and large-scale applications. Integrating specialized algorithms into user-friendly GIS platforms thus remains challenging, and there is the need to extend accessibility for non-expert users and hinder 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 an SAR data reliability evaluation, a susceptibility analysis for shallow landslides, and a 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. Materials and Methods

The Geohazard plugin [28] 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 [29], 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 [23]. 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 [30]; 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.

2.1. 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 monitoring ground deformations and studying landslide phenomena with millimetric precision. Within the Copernicus program [31], data acquired by the Sentinel-1 satellites are freely available through the European Ground Motion Service (EGMS), 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 [32].

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. Input and output for the C index module.

Input	Format	Description
DTM	Raster	DTM of the area
Cell size	Number	DTM cell size (m)
Track angle φ	Number	Track angle of each measurement point; range (0–360°)
Incidence angle θ	Number	Incidence angle of each measurement point; range (0–90°)
Output	Format	Description
C index	Raster	C index computed in each DTM cell (−1; 1)
Percentage of visibility	Raster	VREAL percentage (0; 100)
Class of visibility	Raster	Classification of VREAL: Class 1: 0–20%
		Class 2: 20–40%
		Class 3: 40–60%
		Class 4: 60–80%
		Class 5: 80–100%

). 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 Line Of Sight (LOS), 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, and 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 incorporates 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. It also accounts for topographic parameters, including ground slope α and aspect β [29,33].

Figure 1 illustrates a geometric overview of the algorithm fundamentals. In Figure 1a, the orientation of the satellite’s orbit and the ground are depicted relative to the geographic north [33]. The angles φ and β represent the azimuth directions (aspect) of the satellite orbit and the ground, respectively. Figure 1b describes the relationship between the satellite’s orbit, the LOS, and the ground, where α is the inclination of the ground, θ is the incidence angle of the LOS relative to the ground, and γ is the angle between the ground and the LOS.

The track and incidence angles 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 [34]. 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}(S)·\mathit{sin}\left(A-{\displaystyle \frac{\pi }{2}}\right)\right]+\left[E·(-1)·\mathrm{c}\mathrm{o}\mathrm{s}(S)·\mathit{cos}\left(A-{\displaystyle \frac{\pi }{2}}\right)\right]+\left[H·\mathrm{s}\mathrm{i}\mathrm{n}(S)\right]$$

(1)

where S is the slope of the surface; A is the aspect of the surface with respect to north; and N, E, and H are direction cosines [29]. 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 [35,36]:

$$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. 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 [29].

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, and input data and results are shown in Table 1.

Example of Application and Validation of the Model

For each module, a training dataset is provided within the Geohazard plugin. Validation and demonstration of results were carried out by applying the C index algorithm to the training dataset, related to a high mountain area of 154 × 172 m² in Piemonte (Italy). The DTM used for this analysis has a cell size of 25 m [37] to balance accuracy and computational efficiency. SAR data acquired by the Sentinel-1 satellites [31] were employed with the following orbital parameters: for the ascending orbit, the incidence angle θ was 39° and the track angle (φ + 90°) was 350°; and, for the descending orbit, θ remained 39° with (φ + 90°) set to 190°. Figure 2 illustrates the resulting maps of the C index and visibility classes for both ascending and descending orbits. Slopes with higher visibility in the ascending orbit generally exhibit lower visibility in the descending orbit. Areas where the C index approaches 1 or −1 (indicated as darker green zones in Figure 2a,b) correspond to regions where V_LOS closely approximates V_REAL, aligning with high-visibility classes (Figure 2c,d). Conversely, areas where the C index is near 0 (white zones in Figure 2a,b) correspond to low-visibility classes (Figure 2c,d). Class 5 denotes regions with nearly 100% visibility, where satellite-detected displacement V_LOS closely matches actual ground movement V_REAL. In contrast, class 1 identifies areas with very low visibility, where interpretation of displacement and velocity measurements require greater caution. These results indicate the reliability of satellite data for detecting ground displacements and highlight their utility in small-to-medium-scale land-monitoring activities, enabling the comparison of displacements across areas and the identification of slow slope movements.

To validate the calculations made by the algorithm, an Excel spreadsheet was created, incorporating Equation (1) to derive the C index based on ground aspect β and slope α from the DTM, and incidence angle θ and track angle (φ + 90°) from the orbital parameters. A comparison between the external Excel calculations and the results of the algorithm was conducted on nine cells, with three cells selected for each pair of incidence and track angles. The results confirmed the validity and accuracy of the model, as summarized in

Table 2. Comparison of results for C index validation.

Aspect β (°)	Slope α (°)	Incidence Angle θ (°)	Track Angle φ + 90 (°)	C index (Excel)	C index (Model)
79	15	39	350	0.81	0.81
291	39	39	350	0.04	0.03
252	11	39	350	−0.46	−0.47
159	22	35	310	0.05	0.04
43	48	35	310	0.99	0.99
118	41	35	310	0.63	0.62
133	35	37	330	0.60	0.60
243	12	37	330	−0.42	−0.42
260	24	37	330	−0.19	−0.19

.

2.2. 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 [23,38]. The method (Figure 3) enables the analysis of topographic, soil, and rainfall data to evaluate slope stability and solves the relevant equations using a gridded terrain representation, combining a limit equilibrium stability model for infinite slopes with a steady-state hydrological model. The main assumptions underlying the formulation are as follows [39,40]:

• 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 3a and Equation (4)):

$$F_S={\displaystyle \frac{c^{\prime }+(\gamma ·z·{cos}^2\theta -u)·tan{\phi }^{\prime }}{\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 φ′ 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 considering 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 3b). For a detailed analytical discussion of the hydrological model employed in Shalstab, reference can be made to [38,39].

Figure 3. The hydro-mechanical model underlying the Shalstab module [39]: (a) infinite slope model, where z is the thickness of the shallow layer, h is the height of the saturated layer above the failure surface, and θ is the inclination of the slope; and (b) geometry of the catchment and flow path of water in the hydrological model, where a is the contributing upslope area draining across the contour length b.

Figure 3. The hydro-mechanical model underlying the Shalstab module [39]: (a) infinite slope model, where z is the thickness of the shallow layer, h is the height of the saturated layer above the failure surface, and θ is the inclination of the slope; and (b) geometry of the catchment and flow path of water in the hydrological model, where a is the contributing upslope area draining across the contour length b.

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^{\prime }+(-W·{\gamma }_w)·z·{cos}^2\theta ·tan{\phi }^{\prime }}{\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^{\prime }}{{\gamma }_w·z·sin\theta ·cos\theta }}+{\displaystyle \frac{\gamma }{{\gamma }_w}}·\left(1-{\displaystyle \frac{tan\theta }{tan{\phi }^{\prime }}}\right)\right]·{\displaystyle \frac{T}{q}}·sin\theta$$

(7)

where, according to Figure 3b:

• 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 the following:

$$tan\theta \le \left(1-{\displaystyle \frac{{\gamma }_w}{\gamma }}\right)·tan{\phi }^{\prime }(\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^{\prime }}{\gamma ·z·{cos}^2\theta }}+tan{\phi }^{\prime }(\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^{\prime }}{{\gamma }_w·z·sin\theta ·cos\theta }}+{\displaystyle \frac{\gamma }{{\gamma }_w}}·\left(1-{\displaystyle \frac{tan\theta }{tan{\phi }^{\prime }}}\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. Input for the Shalstab module.

Input	Format	Description
DTM	Raster	DTM of the area
Cell size	Number	DTM cell size (m)
Depth z	Raster	Thickness of the potentially unstable layer (m)
Unit weight γ	Raster	Soil unit weight (N/m3)
Friction angle φ′	Raster	Soil friction angle (°)
Permeability Ks	Raster	Soil permeability coefficient (m/h)
Soil cohesion csoil	Raster	Soil cohesion (N/m2)
Root cohesion croot	Raster	Root cohesion (N/m2)

. 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 the following: (1) absolute stable zones; (2) absolute unstable zones; (3) distribution of critical infiltrated rainfall q_cr (mm/day); and (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. Susceptibility classes.

Class	log(qcr/T) (1/m)
1	−inf	−3.4
2	−3.4	−3.1
3	−3.1	−2.8
4	−2.8	−2.5
5	−2.5	−2.2
6	−2.2	−1.9
7	−1.9	+inf

) as suggested by [38]. 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 [41].

Example of Application and Validation of the Model

To demonstrate the results on a small scale and validate the calculations, the Shalstab module was applied to the Dego area, a region of approximately 400 km² located in the hilly Langhe region, northwestern Italy. This area corresponds to sheet 211 of the Italian 1:50,000 map series produced by the Italian Cartographic Authority (IGM). The Dego area was previously analyzed as part of the Carg Project [42], which aimed to produce thematic geological maps of the Italian territory. In the Carg Project, absolute hazard was defined by incorporating the temporal component, represented by the return period of critical rainfall, and the assessment of the destructive potential of landslide phenomena (intensity). For validation purposes, the Shalstab module utilized input data derived from the Carg Project (

Table 5. Input parameters used for the validation on the Dego area.

Input	Value
DTM	[37]
Cell size	10 m
Depth z	0.52–1.09 m
Unit weigth γ	22,000 N/m3
Friction angle φ′	24–32°
Permeability Ks	0.144 m/h
Soil and root cohesion c	20,000–35,000 N/m2

), obtained through a detailed analysis of parameter variability within the Dego area. The results obtained in terms of absolute stable, absolute unstable cells, and distribution of q_cr were compared.

In Figure 4a, 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. Validation of the Shalstab module via comparison with the Carg project.

ID	qcr Shalstab (mm/day)	qcr Carg (mm/day)	Difference (mm/day)
1	39.04	37.04	−2.00
2	56.13	57.87	−1.74
3	87.76	88.45	−0.69
4	62.68	63.34	−0.66
5	117.52	116.64	0.88
6	37.49	40.90	−3.41
7	87.33	90.57	−3.24
8	61.82	63.12	−1.30
9	19.33	24.70	−5.37
10	216.83	211.98	4.85
11	167.18	164.83	2.35
12	66.49	68.53	−2.04
13	No data	No data	No data
14	77.45	81.81	−4.36
15	146.66	143.79	2.87
16	198.60	190.48	8.12
17	40.20	42.33	−2.13
18	18.89	18.87	0.02
19	3.98	1.76	2.22
20	40.67	40.08	0.59

): 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 Figure 4b, the results are presented for a specific portion of the study area (highlighted as a zoomed-in region in Figure 4a), and compared against the inventory of shallow landslides triggered during the November 1994 rainfall event. This comparison provides additional qualitative validation of the model, particularly given that the cumulative daily rainfall in this region exceeded 200 mm/day during the event, corresponding to a return period of approximately 100 years.

To evaluate the influence of DTM cell size on the analysis results, the calculation was repeated for the same area using the DTM 5 m provided by Regione Piemonte [37]. The parameters listed in Table 5 were used for consistency.

Table 7. Validation of the Shalstab module on different DTMs.

	DTM 10 m		DTM 5 m		Difference
ID	qcr,10 (mm/day)	Class	qcr,5 (mm/day)	Class	qcr,5–qcr,10 (mm/day)
1	39.04	1	112.02	1	72.98
2	56.13	1	52.74	1	−3.39
3	87.76	2	96.64	3	8.88
4	62.68	1	54.19	1	−8.49
5	117.52	3	334.85	7	217.33
6	37.49	1	No data	-	No data
7	87.33	1	352.74	3	265.41
8	61.82	1	84.97	2	23.15
9	19.33	1	0.55	1	−18.78
10	216.83	6	140.22	4	−76.61
11	167.18	3	103.28	2	−63.9
12	66.49	1	124.89	1	58.4
13	No data	-	No data	-	No data
14	77.45	2	80.96	2	3.51
15	146.66	3	173.97	3	27.31
16	198.60	5	285.56	7	86.96
17	40.20	1	77.27	2	37.07
18	18.89	1	13.31	1	−5.58
19	3.98	1	2.43	1	−1.55
20	40.67	1	27.41	1	−13.26

reports the critical rainfall q_cr values for the 20 representative points derived from the analyses using both 10 m and 5 m DTMs.

For comparison, the difference between the q_cr values and the corresponding susceptibility classes are included. An additional comparison is given through the profiles illustrated in Figure 4b, which represent two contour lines located at 450 m (profile 1) and 550 m (profile 2). The q_cr values along these profiles, obtained from their intersection with the respective maps (based on DTM 10 m and DTM 5 m), are shown in Figure 5. The results indicate that the difference in q_cr calculation between DTMs can be significant, particularly in areas where the slope is not uniform. This discrepancy is probably due to the influence of DTM resolution on both the mechanical stability factor Fs and the hydrological component of the model. Variations in slope inclination and catchment area calculations across cells are directly affected by the DTM resolution. These differences affect the susceptibility classification of the area, highlighting the importance of DTM resolution in such analyses. This aspect, which has implications for model accuracy and reliability, will be further investigated in future studies.

2.3. 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 φ_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.

2.3.1. Rockall–Droka_Basic

Droka_Basic module is based on the Cone method [10,21,30], 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 6a): the dip direction angle θ of the slope at the source point defines the orientation of the cone; and the energy angle φ_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 [21]:

$$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), and φ_p is the energy angle.

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

$$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 6a). 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. Input and output for the Droka_Basic module.

Input	Format	Description
DTM	Raster	DTM of the area
Cell size	Number	DTM cell size (m)
Source points	Point vector layer	Selected source points
Mass m	Number	Block mass (kg)
Energy angle φp	Number	Energy angle (°)
Lateral spreading angle α	Number	Lateral spreading angle (°)
Output	Format	Description
Fid, ID	Number	Points identification
Count	Number	Number of overlapped cones
Energy_min, max, mean	Number	Minimum, maximum, mean kinetic energy (kJ)
Velocity_min, max, mean	Number	Minimum, maximum, mean velocity (m/s)
Percent	Number	Number of overlapping cones at the point as a proportion of the total number of generated cones

. 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 in order 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.

2.3.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).

Assuming that the rockfall trajectory remains in constant contact with the ground and follows the steepest downhill gradient, the Droka_Flow module employs the GRASS GIS hydrological function r.drain to compute five paths from each source point. Starting from a user-defined input point or set of points, the algorithm calculates path directions by evaluating the elevation differences between adjacent cells. The output is a raster map delineating the calculated flow paths. Further details on the methodology and implementation are provided in [43]. As with Droka_Basic, the analysis result is a vector point layer, which is then generated, containing one point for each DTM cell covered by the simulated paths. The height difference ΔH between the energy line and the ground is calculated for each point, leading to block velocity and kinetic energy calculations when ΔH is positive (Equations (11) and (12); Figure 6b). 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. Input for the Droka_Flow module.

Input	Format	Description
DTM	Raster	DTM of the area
Cell size	Number	DTM cell size (m)
Source points	Point vector layer	Selected source points
Mass m	Number	Block mass (kg)
Energy angle φp	Number	Energy angle (°)
Mean of the normal distribution	Number	Height difference from DTM (m)
Standard deviation	Number	Standard deviation of the normal distribution of DTM elevation (m)

. Once more, 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 3. The output is a vector layer with point attributes similar to those described in Droka_Basic and listed in Table 8. Frequency data are omitted here due to the limited number of paths per source point, so energy count and percent fields are not included.

2.3.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 parameter 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 7, 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 analysis 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 3.

3. Results

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.

3.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) [44] include the following:

• The detachment occurred at approximately 1940 m a.s.l;
• A block of approximately 30 m³ struck and destroyed a building located at about 1750 m a.s.l., with visible traces along its trajectory (Figure 8a,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 8b.

Based on the available information, a back-analysis was conducted by selecting three source points located within the steepest area of the involved zone (Figure 8b). A block with a volume of 30 m³, mass of 78,000 kg, and unit weight of 2600 kN/m³ was used to estimate the kinetic energy. The analysis utilized a 5 m resolution DTM provided by Regione Piemonte [37] for an initial set of simulations aimed at parameter calibration and a comparison of results between the Droka_Basic and Droka_Flow modules.

For Droka_Basic, the optimal results in terms of the invasion area were achieved with a lateral spreading angle α of 22° and an energy angle φ_p of 37°, as shown in Figure 9, 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 5000 kJ, consistent with the 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.1 m to 5 m (2% to 100% of the DTM cell size). Each simulation was repeated five times to generate about 20 paths (Figure 10). 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 5000 kJ, consistent with the results from Droka_Basic. As shown in Figure 10, 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; and, 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, (available at the Geoportal of the Italian Ministry of the Environment [45]), 5 m, and 10 m DTM resolutions [37] for evaluating the influence of DTM resolution on the results. Figure 11 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 a greater variability in the dip direction at the 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 the DTM resolution. Figure 12 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 12). 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 the susceptibility and hazard analysis can be achieved through the back-analysis conducted within the same geomorphological context and using the same DTM resolution, ensuring consistent results. As discussed, the 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.

3.2. Susceptibility and Relative (Spatial) Hazard Analysis

To exemplify the suggested procedure for the 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 [46], 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 13a. Finally, a block with a volume of 1 m³ was used to estimate the kinetic energy.

Figure 13b compares the results of the analyses conducted with Rockfall–Droka modules in terms of runout area, while Figure 13c,d show the mean kinetic energy obtained by the two modules. To assess the rockfall susceptibility, the Droka_Basic spatial frequency results, defined as the number of overlapping cones out of the total source points (Figure 14a), were multiplied by the mean kinetic energy values (Figure 13c) at each point. The final output, shown in Figure 14b, highlights the most critical zones by combining both the 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 the invasion area (Figure 13b) and mean kinetic energy (Figure 13c,d), as the area covered by rockfall paths simulated by Droka_Flow (Figure 13b) 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 14a), 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 the 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.

4. Discussion

The assessment of the landslide susceptibility and hazard is crucial for territorial management and risk mitigation. At medium to small scales, such as roads, municipalities, or provincial areas, empirical or simplified analytical methods are commonly used to identify critical areas requiring detailed investigation [10]. Numerous methodologies exist for small-scale analyses using territorial data; most have been implemented within GIS environments due to their analytical capabilities and the growing availability of public spatial datasets.

However, there is still a lack of a comprehensive GIS plugin capable of analyzing various landslide types, including slope susceptibility, hazard assessment, affected areas, and associated energy. In this context, the Geohazard plugin was developed as an open-source, user-friendly tool for preliminary landslide susceptibility and hazard assessments at medium to small scales within the QGIS environment. The current version includes three modules tailored to different landslide hazard aspects, described in Section 2, with application examples provided to facilitate the discussion.

The Groundmotion–C Index (Section 2.1) module reliably classified slopes based on their visibility from satellite orbits, effectively managing uncertainties in SAR data interpretation for slow slope movements. Validation using a training dataset confirmed the accuracy of the module, highlighting its utility for the long-term monitoring of gravitational processes.

The Landslide–Shalstab module (Section 2.2) assessed shallow landslide susceptibility by calculating critical rainfall infiltration q_cr across large areas. The analysis produced susceptibility maps that aligned with prior hazard studies, illustrating its reliability. These maps provide a foundation for further analyses, such as evaluating the landslide runout and kinetic energy, which are essential for defining hazard scenarios. However, the application of the Shalstab model to DTMs with different resolutions has revealed its sensitivity to this parameter. High-resolution DTMs can capture fine-scale terrain features, enabling precise calculations of slope gradients and contributing areas. In contrast, lower-resolution DTMs may smooth out critical topographic details, potentially leading to less reliable predictions and the potential overestimation or underestimation of hazard zones [47]. Accurate predictions also depend on the precision and accuracy in estimating hydrological and mechanical parameters and their spatial distribution, a key point when analyzing extensive areas. Therefore, for Shalstab applications, a balanced approach that combines an appropriate DTM resolution with accurate parameter estimations is crucial. This ensures the reliable identification of shallow landslide-prone areas and supports informed decision-making in land-use planning and hazard mitigation. The Rockfall–Droka_Basic and Droka_Flow modules (Section 2.3) proved effective in evaluating rockfall hazards. By calibrating parameters, the modules delineated invasion zones and estimated the kinetic energy, offering complementary perspectives. Droka_Basic provided an overall depiction of the hazard area, while Droka_Flow simulated specific trajectories, enhancing our understanding of the relative hazard dynamics. The influence of the DTM resolution was analyzed in the paper, emphasizing the importance of calibrating model parameters according to DTM cell resolutions. When detailed information about past events is available for the study area, as shown in the Bobbio Pellice case study, performing a back-analysis on the same DTM used for the predictive analysis of the entire region enhances the robustness of the calibration process.

The proposed Geohazard plugin offers an innovative approach to preliminary landslide susceptibility and hazard assessment by integrating multiple analytical modules within the QGIS environment. Unlike previous methodologies that typically address isolated aspects of landslide analysis, this plugin combines slope susceptibility, critical rainfall infiltration, and rockfall impact zones’ evaluation in a unified platform. This integrated framework facilitates a more comprehensive hazard assessment while simplifying the workflow for geohazard studies.

5. Conclusions

The Geohazard plugin, presented in this paper, integrates multiple modules to address different aspects of geohazard analysis, offering a practical and versatile framework for hazard management at medium to small scales.

The Groundmotion–C Index module demonstrated its reliability during validation on a training dataset with a 25 m cell size DTM, where its results closely matched manual calculations. The application results further illustrated the ability to classify slopes based on their visibility from the satellite orbit, providing a valuable tool for reducing uncertainties in SAR data interpretation for the slow slope movement estimation.

The Landslide–Shalstab module effectively assessed the susceptibility to shallow landslides by calculating the critical rainfall infiltration required to trigger slope instability. Its application to the Dego area (located in the hilly Langhe region, northwestern Italy) demonstrated the results of the analysis at a small scale and the validation, based on hazard maps developed in previous studies [39]. The critical rainfall map q_cr and the resulting susceptibility classifications provided valuable insights into the spatial distribution of instability conditions. Further improvements to this method could include the implementation of a simplified infiltration model to express the results in terms of gross rainfall. This information is essential for defining potential landslide scenarios and forms the basis for subsequent analyses of the affected areas and process intensity, contributing to the overall calculation of the landslide hazard.

The Rockfall–Droka_Basic and Droka_Flow modules, designed for rockfall analysis, were tested in the Bobbio Pellice case study, an area in the Italian Western Alps [44]. Parameter calibration based on a well-documented rockfall event in the Bobbio Pellice area enabled the accurate delineation of invasion zones and the 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. The results showed that the destroyed building in the area was in a critical zone, where the estimated block kinetic energy exceeded 5000 kJ. The case study showed how the modules perform under real case studies. The calibration process, grounded in field observations, emphasized the critical importance of systematically aligning model inputs with observed data to enhance their accuracy. Additionally, the results showed that the modules can replicate observed events, offering important insights into trajectory precision and hazard mapping. The findings also stressed the importance of maintaining consistency between the DTM used for calibration and the one applied in the analysis to ensure consistent and reliable results, supporting better hazard prediction and land-use planning. Future research will focus on addressing these sensitivities to further improve the applicability and robustness of the plugin. The two modules were complementary, with Droka_Basic delineating the overall invasion area and Droka_Flow providing detailed insights into critical trajectories.

These findings validate the Geohazard plugin as a reliable and versatile tool for geohazard assessment. By integrating multiple analytical approaches, the plugin enhances our understanding of slope instability and aids in developing effective risk mitigation strategies. Looking ahead, the Geohazard plugin will be expanded to include additional modules specifically tailored for territorial landslide analysis, with a particular focus on diffuse and recurrent phenomena. Planned improvements include tools for slope classification related to rockfall susceptibility, defined as the propensity of a slope to trigger rockfall events. Additionally, the integration of triggering and run-out modules will provide a comprehensive analytical framework suitable for various types of landslides. A simplified tool for assessing the vulnerability of different types of exposed elements will also be introduced, enabling a preliminary estimation of landslide risk at medium to small scales.

Finally, the possibility of the integration of machine-learning modules for landslide susceptibility mapping [48] will be explored to enhance predictive accuracy and facilitate the incorporation of diverse environmental and geospatial datasets into the analysis.

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 [49].